pro inter,wi,wf,t,g,f,x,y,yerr,res,x1d,y1d,path=path,plot=plot

;+
;
;  3D spectra interpolation of the CIFIST grid (see Bertran de Lis et al. 2017)
;
;	IN:  
;            wi    -- float          Initial wavelegth (Angstroms)
;	         wf    -- float          Final wavelegth (Angstroms)
;            t     -- float          Effective temperature
;	         g     -- float          Log of the surface gravity
;            f     -- float          Metallicity [Fe/H]
;            res   -- float          If 1D is provided, this is the resolving power of the 1D
;                                    spectra, and is mandatory. Otherwise, is the desired output
;                                    resolution, and it is optional.
;            x1d   -- 1d float array Optional. 1D wavelength array. If provided, the quotient 3D/1D
;                                    is interpolated, and the the result multilied by y1d. Otherwise, 
;                                    straight 3D interpolation is performed.
;            y1d   -- 1d float array Optional. 1D normalized flux array
;            path  -- string         Optional. Path where "*.fits" files are located
;
;	OUT: 
;            x     -- float          Output wavelegth (Angstroms)
;	         y     -- float          Output 3D normalized flux
;            yerr  -- float          Uncertainty of output 3D flux		
;
;	KEYWORDS: 
;            plot  -- bool           When set, several plots are produced in "Interpolchecks_xxx.eps" in 
;                                    order to help the user to evaluate the quality and effectiveness 
;                                    of the interpolation:
;                                    a) Countour maps of the mean flux ratio 3D/1D averaged
;                                       at all requested wavelengths, for each [Fe/H].                        
;                                    b) Contour maps of the RMS of the difference between the 3D 
;                                       and 1D normalized fluxes at the wavelength range requested, 
;                                       for each [Fe/H].
;                                    c) Contour maps of the RMS between the 3D normalized flux and 
;                                       its interpolated 3D (removing the corresponding point of the 
;                                       grid), at the wavelength range requested and for each [Fe/H].
;                                    d) Contour maps of the ratio between plot b) and c). When 1D 
;                                       spectra aproximates better to real 3D than the 3D interpolation,
;                                       due to large errors, it appears in grey scale. It is not 
;                                       recommended to use the interpolation in this case.
;                                    e) Contour maps of the mean absolute uncertainty averaged at all 
;                                       requested wavelengths, for each [Fe/H].
;
;       USAGE: inter,8535,8548,5200,4.2,-0.5,x,y,yerr,50000,/plot
;
;
; S. Bertran de Lis, IAC, May, 2017
;
;-

;Check input parameters
if N_params() lt 8 then begin
  print,'% inter: use -- inter,wi,wf,t,g,f,x,y,yerr[,res,x1d,y1d,path],plot=plot'
  return
endif

;Check wavelength range
if (wi ge wf) then begin
   print,'% inter: error -- Initial wavelength should be smaller than final wavelegth'
   return
endif
if (wi lt 3000. or wi gt 30000. or wf lt 3000. or wf gt 30000.) then begin
   print,'% inter: error -- Wavelength range should be between 3000 and 30000A'
   return
endif

;Check range input stellar parameters
if (t lt 4000.  or t gt 6500.) then begin
   print,'% inter: error -- input Teff must be between 4000 and 6500 K'
   return
endif
if (g lt 1.  or g gt 4.9) then begin
   print,'% inter: error -- input log(g) must be between 1.7 and 4.5 dex'
   return
endif
if (f lt -3.0 or f gt 0.5) then begin
   print,'% inter: error -- input [Fe/H] must be between -3.0 and 0.0 dex'
   return
endif

;Check 1D input spectra
sx1d=size(x1d)
sy1d=size(y1d)
if (sx1d[0] gt 1 or sy1d[0] gt 1) then begin
  print,'% inter: error -- 1D wavelength and normalized flux should be 1D arrays'
  return
endif
if (sx1d[0] eq 1 and sy1d[0] eq 1) then begin
   if (x1d[0] gt wi or x1d[sx1d[1]-1] lt wf) then begin
      print,'% inter: error -- 1D wavelength array do not cover the requested wavelength array'
      return
   endif
   ;if(n_elements(res) ne 1) then begin
   ;   print,'% inter: error -- A resolution is required when 1D spectra is provided'
   ;   return
   ;endif
endif

;Interpolation mode
if (sx1d[0] eq 1 and sy1d[0] eq 1) then begin
   print,'-- 1D spectra found. 3D/1D ratio interpolation is performed' 
   intmode = '3d1d'
endif else begin
   print,'-- 1D spectra not found. Straight 3D interpolation is performed'
   intmode = '3d'
endelse

;Check path
if (n_elements(path) eq 0) then path='./'
if (intmode eq '3d1d') then tables = ['W3D1D.fits', 'W3D1Derr.fits']
if (intmode eq '3d') then tables = ['W1D.fits', 'W3D1Derr.fits']
for i=0,1 do begin
   if (file_test(path+tables[i]) eq 0) then begin
      print,'% inter: error -- File `'+tables[i]+'` not found in directory specified in path parameter'
      return
   endif
endfor
if keyword_set(plot) then begin
   for i=2,n_elements(tables)-1 do begin
      if (file_test(path+tables[i]) eq 0) then begin
         print,'% inter: error -- File `'+tables[i]+'` not found in directory specified in path parameter'
         return
      endif
   endfor
endif

;Open & read data table
if (intmode eq '3d1d') then begin 
   tabin  = path+'W3D1D.fits' 
   tabein = path+'W3D1Derr.fits'
endif else begin 
   tabin  = path+'W3D1D.fits'
   tabein = path+'W3D1Derr.fits'
   tab1din = path+'W1D.fits'
endelse
tab  = readfits(tabin,EXTEN_NO=1,/SILENT)
if intmode ne '3d1d' then tab1d=readfits(tab1din, EXTEN_NO=1,/SILENT)

w    = tab[0,3:*]
teff = tab[1:*,0]
logg = tab[1:*,1]
feh  = tab[1:*,2]

;Interpolation
wc = where(w ge wi and w le wf)
x  = w[wc]
nx = n_elements(x)
y  = fltarr(nx)
p  = [(transpose(teff)-min(teff))/(max(teff)-min(teff)), (transpose(logg)-min(logg))/(max(logg)-min(logg)), $
      (transpose(feh)-min(feh))/(max(feh)-min(feh)) ] ; Normalized stellar parameters
pp = [(t-min(teff))/(max(teff)-min(teff)), (g-min(logg))/(max(logg)-min(logg)), $
      (f-min(feh))/(max(feh)-min(feh)) ] ;Interpolation point
for i=0,nx-1 do begin
   rbf_eval,tab[1:*,3+wc[i]],p,pp,out
   if intmode eq '3d1d' then y[i] = out else begin
     rbf_eval,tab1d[1:*,3+wc[i]],p,pp,out1d
     y[i] = out*out1d
   endelse
endfor
tab = 0

;Uncertainties
taberr = readfits(tabein,EXTEN_NO=1,/SILENT)
yerr   = fltarr(nx)
for i=0,nx-1 do begin
   rbf_eval,taberr[1:*,3+wc[i]],p,pp,out
   yerr[i]=out/2. ;Divided by two: Reduction factor
endfor
taberr = 0
print,'-- Interpolation finished!'

;Smooth spectra
if(n_elements(res) eq 1) then begin
  print,'-- Smoothing'
  fwhm = 1./res/alog(10.)
  gconv,x,y,fwhm,xres,yres
  gconv,x,yerr,fwhm,xres,yerrres
  x = xres
  y = yres
  yerr = yerrres
endif
if (intmode eq '3d1d') then begin
  wc1d = where(x1d ge wi and x1d le wf)
  qint = interpol(y,x,x1d[wc1d],/SPLINE)
  eint = interpol(yerr,x,x1d[wc1d],/SPLINE)
  y = qint*y1d[wc1d]
  x = x1d[wc1d]
  yerr = eint
endif

;Plots
if keyword_set(plot) then begin

  print,'-- Starting to plot. This may take few minutes...'   

  ;Set output
  rand = SYSTIME(/SECONDS)
  output = 'Interpolchecks_'+string(rand,FORMAT='(F13.2)')+'.ps'
  SET_PLOT,'PS'
  DEVICE,filename=output,landscape=0,/COLOR,xsize=20,ysize=27,yoffset=0,xoffset=0
  !P.MULTI=[0,4,4,0,0]
  charsize = 1.5

  ; a) Contours mean(3D/1D)
  f3d1d   = readfits(path+'F3D1D.fits',EXTEN_NO=1,/SILENT)
  f3d1d   = f3d1d[1:*,3+wc]
  nlevels = 50.
  minlev = mean(f3d1d)-2*stddev(f3d1d)
  if (minlev lt 0) then minlev  = 0.
  maxlev  = mean(f3d1d)+2*stddev(f3d1d)
  levels  = IndGen(nlevels)*((maxlev-minlev))/nlevels+minlev
  cgLoadCT,33,ncolors=nlevels,bottom=3
  cgColorBar,ncolors=nlevels,bottom=3,divisions=10,range=[minlev,maxlev],format='(F4.2)',charsize=charsize,$
    tcharsize=0.8,annotatecolor='black',title='Mean(F3D/F1D)',tlocation='top',position=[.08,0.97,.47,0.98] 
  metal = [0.,-1., -2., -3.]
  for i=0,n_elements(metal)-1 do begin 
     filter = where(feh eq metal[i])
     data = mean(f3d1d[filter,*],dimension=2)
     cgcontour,data,teff[filter]/1000.,logg[filter],/irregular,/fill,charsize=charsize,/noerase,$
      c_colors=indgen(nlevels)+3,nlevels=nlevels,levels=levels,resolution=[28,16],label=0,$
      xtitle='T!Deff!N',ytitle='log(g)',title=strjoin(['[Fe/H] =',string(metal[i],format='(F4.1)')],' '),$
      xrange=[7.,3.7],yrange=[4.5,1.5],Position=[0.08+i*0.23,0.81,0.25+i*0.23,0.93]      
  endfor

  ; b) Contours RMS(3D-1D)
  f3d = readfits(path+'F3D.fits',EXTEN_NO=1,/SILENT)
  f3d = f3d[1:*,3+wc]
  f1d = f3d/f3d1d
  rms     = sqrt(total((f3d-f1d)^2.,2)/nx)
  minlev  = 0
  maxlev  = mean(rms)+2*stddev(rms)
  levels  = IndGen(nlevels)*((maxlev-minlev))/nlevels+minlev
  cgLoadCT,33,ncolors=nlevels,bottom=3
  cgColorBar,ncolors=nlevels,bottom=3,divisions=6,range=[minlev,maxlev],format='(F6.3)',charsize=charsize,$
    tcharsize=0.8,annotatecolor='black',title='RMS(F3D-F1D)',tlocation='top',position=[.08,0.73,.47,0.74]
  for i=0,3 do begin
    filter  = where(feh eq metal[i])
    data    = sqrt(total((f3d[filter,*] - f1d[filter,*])^2.,2)/nx)
    cgcontour,data,teff[filter]/1000.,logg[filter],/irregular,/fill,charsize=charsize,/noerase,$
      c_colors=indgen(nlevels)+3,nlevels=nlevels,levels=levels,resolution=[28,16],background='black',$
      xtitle='T!Deff!N',ytitle='log(g)',title=strjoin(['[Fe/H] =',string(metal[i],format='(F4.1)')],' '),$
      xrange=[7.,3.7],yrange=[4.5,1.5],Position=[0.08+i*0.23,0.57,0.25+i*0.23,0.69]
  endfor

  ; c) Contours RMS(3D-3Dint)
  if (intmode eq '3d1d') then f3derr  = readfits(path+'F3D1Derr.fits',EXTEN_NO=1,/SILENT)
  if (intmode eq '3d')   then f3derr  = readfits(path+'F3Derr.fits',EXTEN_NO=1,/SILENT)
  f3derr  = f3derr[1:*,3+wc]
  rms     = sqrt(total((f3derr)^2.,2)/nx)
  minlev  = 0
  maxlev  = mean(rms)+2*stddev(rms)
  levels  = IndGen(nlevels)*((maxlev-minlev))/nlevels+minlev
  cgLoadCT,33,ncolors=nlevels,bottom=3
  cgColorBar,ncolors=nlevels,bottom=3,divisions=6,range=[minlev,maxlev],format='(F6.3)',charsize=charsize,$
    tcharsize=0.8,annotatecolor='black',title='RMS(F3D-F3D!Dint!N)',tlocation='top',position=[.08,0.49,.47,0.50]
  for i=0,3 do begin
    filter  = where(feh eq metal[i])
    data    = sqrt(total((f3derr[filter,*])^2.,2)/nx)
    cgcontour,data,teff[filter]/1000.,logg[filter],/irregular,/fill,charsize=charsize,/noerase,$
      c_colors=indgen(nlevels)+3,nlevels=nlevels,levels=levels,resolution=[28,16],background='black',$
      xtitle='T!Deff!N',ytitle='log(g)',title=strjoin(['[Fe/H] =',string(metal[i],format='(F4.1)')],' '),$
      xrange=[7.,3.7],yrange=[4.5,1.5],Position=[0.08+i*0.23,0.33,0.25+i*0.23,0.45]
  endfor
  
  ; d) Contours RMS(3D-1D)/RMS(3D-3Dint)
  rms = sqrt(total((f3d-f1d)^2.,2)/nx) / sqrt(total((f3derr/2)^2.,2)/nx);Divided by two: Reduction factor
  minlev  = 0
  maxlev  = 10
  levels  = IndGen(nlevels)*((maxlev-minlev))/nlevels+minlev
  rcolor = [ 192+(0-192)*findgen(5)/(5-1), replicate(255,20), 255+(0-255)*findgen(30)/(30-1)] 
  gcolor = [ 192+(0-192)*findgen(5)/(5-1), 255+(0-255)*findgen(20)/(20-1), replicate(0,30)]
  bcolor = [ 192+(0-192)*findgen(5)/(5-1), replicate(0,20), 0+(255-0)*findgen(30)/(30-1)]
  TVLCT,rcolor,gcolor,bcolor
  cgColorBar,ncolors=nlevels,bottom=0,divisions=10,range=[minlev,maxlev],format='(F5.1)',charsize=charsize,$
    tcharsize=0.8,annotatecolor='black',title='RMS(F3D-F1D)/RMS(F3D-F3D!Dint!N)/2',tlocation='top',position=[.08,0.25,.47,0.26]
  for i=0,3 do begin
    filter  = where(feh eq metal[i])
    difreal = f3d[filter,*] - f1d[filter,*]
    difint  = f3derr[filter,*]/2;Divided by two: Reduction factor
    data    = sqrt(total((difreal)^2.,2)/nx) / sqrt(total((difint)^2.,2)/nx)
    cgcontour,data,teff[filter]/1000.,logg[filter],/irregular,/fill,charsize=charsize,/noerase,$
      c_colors=indgen(nlevels),nlevels=nlevels,levels=levels,resolution=[28,16],$
      xtitle='T!Deff!N',ytitle='log(g)',title=strjoin(['[Fe/H] =',string(metal[i],format='(F4.1)')],' '),$
      xrange=[7.,3.7],yrange=[4.5,1.5],Position=[0.08+i*0.23,0.09,0.25+i*0.23,0.21]
  endfor
  
  ;Change page to plot
  ERASE

  ; e) Contours uncertainties
  minlev  = 0.
  maxlev  = mean(f3derr/2.)+2*stddev(f3derr/2.) ;Divided by two: Reduction factor
  levels  = IndGen(nlevels)*((maxlev-minlev))/nlevels+minlev
  cgLoadCT,33,ncolors=nlevels,bottom=3
  cgColorBar,ncolors=nlevels,bottom=3,divisions=6,range=[minlev,maxlev],format='(F6.3)',charsize=charsize,$
    tcharsize=0.8,annotatecolor='black',title='Mean(Abs(F3D - F3D!Dint!N))/2',tlocation='top',position=[.08,0.95,.47,0.96]
  for i=0,3 do begin
    filter  = where(feh eq metal[i])
    data    = mean(abs(f3derr[filter,*]/2.),dimension=2)
    cgcontour,data,teff[filter]/1000.,logg[filter],/irregular,/fill,charsize=charsize,/noerase,$
      c_colors=indgen(nlevels)+3,nlevels=nlevels,levels=levels,resolution=[28,16],background='black',$
      xtitle='T!Deff!N',ytitle='log(g)',title=strjoin(['[Fe/H] =',string(metal[i],format='(F4.1)')],' '),$
      xrange=[7.,3.7],yrange=[4.5,1.5],Position=[0.08+i*0.23,0.79,0.25+i*0.23,0.91]
  endfor
  f3derr = 0
 
  ; f) Interpolated spectra
  err_high = y + abs(yerr)
  err_low  = y - abs(yerr)

  if (intmode eq '3d') then begin
     w1d   = readfits(path+'W1D.fits',EXTEN_NO=1,/SILENT)
     y1d   = fltarr(nx)
     for i=0,nx-1 do begin
       rbf_eval,w1d[1:*,3+wc[i]],p,pp,out
       y1d[i] = out
     endfor
     x1d = x
     w1d = 0
  endif
  
  cgplot,x1d,y1d,color='blue',charsize=charsize,xtitle='Wavelength ('+string(197B)+')', ytitle='Norm Flux',/noerase,$
      position=[0.08,0.58,0.71,0.70],xrange=[wi,wf],xstyle=1,yrange=[min(y)-0.15*min(y),max(y)+0.1*max(y)],ystyle=1,/normal
  cgcolorfill,[x,reverse(x),x[0]],[err_high,reverse(err_low),err_high[0]],color='gray'
  cgoplot,x,y,color='red'
  cgoplot,x1d,y1d,color='blue'
  cgLegend, title=[' '],LineStyle=[0], Colors=['gray'],thick=20,length=0.04,Location=[0.1,0.597]
  cgLegend, title=['1D','3D!9'+string("53B)+'!X3D!Derr!N'], LineStyle=[0,0], Colors=['blue','red'],vspace=1.,$
    length=0.04,charsize=.7,Location=[0.1,0.61]

    
  cgLoadCT,1
  !P.COLOR=0
  xyouts,0.77,0.68,strjoin(['T!Deff!N    =',string(t,format='(I4)')],' '),/normal
  xyouts,0.77,0.66,strjoin(['Log(g) =',string(g,format='(F3.1)')],' '),/normal
  xyouts,0.77,0.64,strjoin(['[Fe/H] =',string(f,format='(F4.1)')],' '),/normal
  xyouts,0.77,0.62,strjoin(['Wave_i =',string(wi,format='(F8.2)')],' '),/normal
  xyouts,0.77,0.60,strjoin(['Wave_f =',string(wf,format='(F8.2)')],' '),/normal

  DEVICE, /close_file
  set_plot, 'X'
  print,'Done!'

endif

end