Pro tell_corr,flux_file,corrflux_file,error_file=error_file,$
              correrror_file=correrror_file,wave_file=wave_file,$
              fwhm=fwhm,newextinct=newextinct
  
;+---------------------------------------------------------------------
;
; TELL_CORR
;
; Interactively perform telluric absorption correction using model
; transmission spectra.
;
; Telluric transmission models for Paranal calculated using SKYCALC,
; Cerro Paranal Advanced Sky Model, version 1.3.5.  Resolution R =
; 100,000.  Website: http://www.eso.org/observing/etc/bin/gen/
; form?INS.MODE=swspectr+INS.NAME=SKYCALC
;
; Models are a function of airmass and precipitable water vapor.  The
; user interactively adjusts these features, along with a velocity
; offset and smoothing FWHM, in order to best match the telluric
; absorption in the observed spectrum.
;
; A list of commands can be obtained by typiing '?' at the prompt.
;
; INPUTS
;         flux_file   - String, FITS filename of 1D uncorrected flux
;                       spectrum 
;         
; KEYWORDS
;         error_file  - String, FITS filename of 1D uncorrected error
;                       spectrum.  If set along with the
;                       correrror_file keyword, the telluric
;                       correction will also be applied to the error
;                       spectrum.
;         correrror_file - String, name of FITS file to which to write
;                          corrected error spectrum.
;         wave_file   - String, FITS filename of 1D wavelength array
;                       for uncorrected spectrum.  Default is to
;                       determine wavelengths from header keywords.
;                       NOTE: Wavelength bins must have a fixed
;                       velocity width.
;         fwhm        - Scalar, initial FWHM (km/s) of Gaussian kernel
;                       by which to smooth model transmission
;                       spectrum.  Default is to set this
;                       interactively.
;         newextinct  - FOR X-SHOOTER LP DATA.  If set, the extinction
;                       correction applied in xshooter_multiextract is
;                       undone before the new telluric correction is
;                       applied.
;
; OUTPUTS
;         corrflux_file - String, name of FITS file to which to write
;                         corrected flux spectrum.
;
; HISTORY
;         Written 12/19/2013 GDB
;----------------------------------------------------------------------


if (n_params() lt 2) then begin
   print,'CALLING SEQUENCE:'
   print,'    tell_corr,flux_file,corrflux_file,error_file=error_file,'
   print,'              correrror_file=correrror_file,wave_file=wave_file,'
   print,'              fwhm=fwhm'
   return
end

;;;
;;; Read in object 
;;;

; Read in the uncorrected spectrum
if keyword_set(wave_file) then begin
   rawflux = readfits(flux_file,hflux)
   wave    = readfits(wave_file)
endif else begin
   rawflux = rflam(flux_file,hflux,wave)
endelse

; Optional - Read in uncorrected error array
if keyword_set(error_file) then begin
   rawerror = readfits(error_file,herror)
endif else begin
   rawerror = 0
   herror   = ''
endelse

; Inidices of object spectrum
obj_npix = n_elements(rawflux)
obj_ind  = dindgen(obj_npix)

; Check that wavelengths have costant velocity widths
dloglam_lo = alog10(wave(1))-alog10(wave(0))
dloglam_hi = alog10(wave(obj_npix-1))-alog10(wave(obj_npix-2))
if (abs((dloglam_lo-dloglam_hi)/dloglam_lo) gt 1d-5) then begin
   print,'Object pixels must have constant velocity width.'
   return
endif

; Object pixel size
c_kms      = 2.99792458d5
obj_pix_dv = c_kms * (10d^dloglam_lo - 1)

; Wavelengths at the pixel edges
wave_lo = wave * (1 - 0.5*obj_pix_dv/c_kms)
wave_hi = wave * (1 + 0.5*obj_pix_dv/c_kms)

; Get airmass from header
object_airmass = gxpar(hflux,'AIRMASS',count=n_airmass)
if (n_airmass eq 0) then $
   object_airmass = gxpar(hflux,'HIERARCH ESO TEL AIRM END',count=n_airmass)
if (n_airmass eq 0) then begin
   print,'Airmass unknown.  Assuming 1.0'
   object_airmass = 1.0
endif

; Get heliocentric velocity information from header
vhelio = gxpar(hflux,'HELIOVEL',count=n_vhelio)
if (n_vhelio eq 0) then begin
   print,'Heliocentric velocity unknown.  Assuming 0.0'
   vhelio = 0.0
endif

; Optional - Undo previous atmospheric extinction correction
; This may be necessary for X-Shooter LP data.
if keyword_set(newextinct) then begin
   print,'Undoing previous extinction correction...'
   ; Extinction model as implemented in xshooter_multiextract
   ext_dir   = '/Users/gdb/eso/xshooter/flux_cal/calibrations/'
   ext_file  = ext_dir+'xsh_paranal_extinct_model_uvb.fits'
   ext_model = mrdfits(ext_file,1,ext_hdr)
   ext_w     = 10.*ext_model.lambda   ; Angstroms
   ext_k     = ext_model.extinction   ; mag/airmass
   ; Assume a constant value at wavelength higher than those modeled
   nbin_ext  = n_elements(ext_w)
   ext_dw    = ext_w(nbin_ext-1) - ext_w(nbin_ext-2)
   ext_w     = [ext_w,ext_w(nbin_ext-1)+ext_dw]
   ext_k     = [ext_k,ext_k(nbin_ext-1)]
   ; Atmospheric extinction
   corr_ls    = where(wave gt 0)
   k          = fltarr(obj_npix)
   k(corr_ls) = interpol(ext_k,ext_w,wave(corr_ls))
   frac_trans = 1. + fltarr(obj_npix)
   frac_trans(corr_ls) = 10.^(-0.4*object_airmass*k(corr_ls))
   ; Undo extinction correction
   rawflux    = rawflux  * frac_trans
   rawerror   = rawerror * frac_trans
endif

;;;
;;; Read in telluric absorption models
;;;

; Directory in which telluric models are stored
model_dir = '/Users/gdb/eso/xshooter/flux_cal/telluric_models/'

; Model air masses
airmass_list = [1.0, 1.5, 2.0, 2.5, 3.0]
airmass_str  = '_air'+['1.0', '1.5', '2.0', '2.5', '3.0']
n_airmass    = n_elements(airmass_list)

; Model PWV
pwv_list = [0.5, 1.0, 1.5, 2.5, 3.5, 5.0, 7.5, 10.0, 20.0]
pwv_str  = '_pwv'+['0.5', '1.0', '1.5', '2.5', '3.5', '5.0', '7.5', '10.0', '20.0']
n_pwv    = n_elements(pwv_list)

; Read in transmission models
print,'Reading in telluric models...'
for i=0,n_airmass-1 do begin
   for j=0,n_pwv-1 do begin
      model_file = model_dir+'skytable'+airmass_str(i)+pwv_str(j)+'.fits'
      model_data = mrdfits(model_file,1,hdr,/silent)
      if (i eq 0 and j eq 0) then begin
         model_wave      = 10000 * model_data.lam
         model_npix      = n_elements(model_wave)
         all_model_trans = dblarr(model_npix,n_pwv,n_airmass)
      endif
      all_model_trans(*,j,i) = model_data.trans
   endfor
endfor

; Use only wavelengths that cover the object wavelengths
keep_ls         = where(model_wave ge min(wave)-50 and $
                        model_wave le max(wave)+50,n_keep)
all_model_trans = all_model_trans(keep_ls,*,*)
model_wave      = model_wave(keep_ls)
model_npix      = n_keep

; Indices of model pixels
model_ind = dindgen(model_npix)

; Pixel size for transmission models.  For Paranal models these will
; have a constant size in km/s.  Interpolation routine below currently
; requires constant velocity width.
model_dloglam = alog10(model_wave(1)) - alog10(model_wave(0))
model_pix_dv  = c_kms * (10d^model_dloglam - 1)

; Number of model pixels per object pixels
pixsz_ratio = obj_pix_dv / model_pix_dv

; Substitute small values for zeros (there should be no negative
; values)
zero_ls = where(all_model_trans eq 0,n_zero)
if (n_zero ne 0) then all_model_trans(zero_ls) = 10d^(-100)

; Interpolate over logarithmic values
all_model_logtrans = alog10(all_model_trans)


;;;
;;; Interactively find best-fitting telluric model
;;;

; Set initial airmass, PWV, velocity offset, and instrumental FWHM
airmass = object_airmass
pwv     = 2.5
dv      = vhelio  ; km/s
if keyword_set(fwhm) then ker_fwhm = fwhm else ker_fwhm = 30.0  ; km/s
   
; Set initial display parameters
!p.multi    = [0,1,3]
!p.charsize = 2.5
!x.style    = 1
!y.style    = 1
xmin        = min(wave)
xmax        = max(wave)
ymin        = min(rawflux)
ymax        = max(rawflux)
obj_plot_ls = where(wave ge xmin and wave le xmax)
device,decomp=0
loadct,39

; Transmission model integrated over object wavelength bins
binned_model_trans = fltarr(obj_npix)

option = ''

while (option ne 'q') do begin

   ; Get indices of closest model airmasses
   if (airmass ge min(airmass_list)) then begin
      airmass_ind_lo = max(where(airmass_list le airmass))
   endif else begin
      airmass_ind_lo = 0
   endelse
   airmass_ind_lo <= (n_airmass - 2)
   airmass_ind_hi = airmass_ind_lo + 1
   
   ; Interpolate models over airmass
   airmass_lo        = airmass_list(airmass_ind_lo)
   airmass_hi        = airmass_list(airmass_ind_hi)
   airmass_weight_lo = (airmass_hi - airmass) / (airmass_hi - airmass_lo)
   airmass_weight_hi = 1 - airmass_weight_lo
   airmass_model_logtrans = dblarr(model_npix,n_pwv)
   for j=0,n_pwv-1 do airmass_model_logtrans(*,j) = $
                   airmass_weight_lo*all_model_logtrans(*,j,airmass_ind_lo) + $
                   airmass_weight_hi*all_model_logtrans(*,j,airmass_ind_hi)
   
   ; Get indices of closest PWV
   if (pwv ge min(pwv_list)) then begin
      pwv_ind_lo = max(where(pwv_list le pwv))
   endif else begin
      pwv_ind_lo = 0
   endelse
   pwv_ind_lo <= (n_pwv - 2)
   pwv_ind_hi = pwv_ind_lo + 1
   
   ; Interpolate models over pwv
   pwv_lo         = pwv_list(pwv_ind_lo)
   pwv_hi         = pwv_list(pwv_ind_hi)
   pwv_weight_lo  = (pwv_hi - pwv) / (pwv_hi - pwv_lo)
   pwv_weight_hi  = 1 - pwv_weight_lo
   model_logtrans = pwv_weight_lo*airmass_model_logtrans(*,pwv_ind_lo) + $
                    pwv_weight_hi*airmass_model_logtrans(*,pwv_ind_hi)
   model_trans    = 10d^model_logtrans
   
   ; Smooth to instrumental resolution
   ker                = gaussker(fwhm=ker_fwhm/model_pix_dv)
   smooth_model_trans = convol(model_trans,ker,/edge_truncate)
   
   ; Shift model wavelengths
   shifted_model_wave    = model_wave * (1 + dv/c_kms)

   for j=0d,obj_npix-1 do begin
      ; Occasionally calculate a zero point for model indices
      ; corresponding to the low-wavelength edge of a relatively
      ; nearby object pixel.  This needs to be done more than once to
      ; avoid drift due to floating point errors.
      if (j mod 1000 eq 0) then begin
         j0   = j
         ind0 = interpol(model_ind,shifted_model_wave,wave_lo(j))
      endif
      ; Model indices corresponding to the lower and upper edges
      ; of the current bin
      ind_loedge = ind0 + (j-j0)*pixsz_ratio
      ind_hiedge = ind0 + (j-j0+1)*pixsz_ratio
      ; Index range of input pixels that fall at least partially in
      ; the current bin
      ind_lo   = ceil(ind_loedge - 0.5)
      ind_hi   = floor(ind_hiedge + 0.5)
      ind_list = model_ind(ind_lo:ind_hi)
      ; Fraction of each model pixel that falls in the current bin
      pix_loedge = (ind_list - 0.5) > ind_loedge
      pix_hiedge = (ind_list + 0.5) < ind_hiedge
      frac       = pix_hiedge - pix_loedge
      ; Sum counts
      binned_model_trans(j) = total(smooth_model_trans(ind_list)*frac) / $
                              total(frac)
   endfor

   ; Apply telluric correction
   corr_flux = rawflux / binned_model_trans

   ; Rough scaling for overplotting model transmission
   model_scale   = median(rawflux(obj_plot_ls)) / $
                   median(binned_model_trans(obj_plot_ls))

   ; Plot unsmoothed telluric model
   plot,shifted_model_wave,model_trans,xr=[xmin,xmax],yr=[-0.1,1.1],$
        psym=10,title='Model Transmission'
   oplot,[-1e6,1e6],[0,0],linestyle=1
   oplot,[-1e6,1e6],[1,1],linestyle=1

   ; Plot uncorrected flux with smoothed telluric model overlaid
   plot,wave(obj_plot_ls),rawflux(obj_plot_ls),xr=[xmin,xmax],yr=[ymin,ymax],$
        psym=10,title='Uncorrected Flux + Smoothed/Binned Model (Scaled)'
   oplot,[-1e6,1e6],[0,0],linestyle=1
   oplot,wave(obj_plot_ls),model_scale*binned_model_trans(obj_plot_ls),$
         psym=10,color=250

   ; Plot corrected flux
   plot,wave(obj_plot_ls),corr_flux(obj_plot_ls),xr=[xmin,xmax],$
        yr=[ymin,ymax],psym=10,title='Corrected Flux'
   oplot,[-1e6,1e6],[0,0],linestyle=1

   ; Print current parameters
   print,'Airmas   = ',airmass
   print,'PWV      = ',pwv
   print,'Delta_v  = ',dv
   print,'FWHM     = ',ker_fwhm

   ; Get command
      print,'Enter option (? for list):'
      option = get_kbrd(1)
   case option of
      '6': airmass -= 0.1 
      '7': airmass -= 0.01
      '8': airmass += 0.01
      '9': airmass += 0.1 
      'y': pwv -= 0.1 
      'u': pwv -= 0.01
      'i': pwv += 0.01
      'o': pwv += 0.1 
      'h': dv -= 1. 
      'j': dv -= 0.1
      'k': dv += 0.1
      'l': dv += 1. 
      ',': ker_fwhm -= 0.5
      '.': ker_fwhm += 0.5
      'b': begin
             read,airmass,prompt='Enter airmass: '
             read,pwv,prompt='Enter PWV: '
             read,dv,prompt='Enter Delta_v: '
             read,ker_fwhm,prompt='Enter FWHM: '
           end
      'z': begin
             read,ymin,prompt='Enter min y-value: '
             read,ymax,prompt='Enter max y-value: '
          end
      'x': begin
             read,xmin,prompt='Enter min x-value: '
             read,xmax,prompt='Enter max x-value: '
             obj_plot_ls = where(wave ge xmin and wave le xmax)
          end
      'v': begin
             print,'Click on plot twice for new x bounds.'
             getpoints,temp_x,temp_y
             if (n_elements(temp_x) ge 2 and $
                 min(temp_x) ne max(temp_x)) then begin
                xmin = min(temp_x)
                xmax = max(temp_x)
                obj_plot_ls = where(wave ge xmin and wave le xmax)
             endif
          end
      'c': begin
             print,'Click on plot twice for new y bounds.'
             getpoints,temp_x,temp_y
             if (n_elements(temp_y) ge 2 and $
                 min(temp_y) ne max(temp_y)) then begin
                ymin = min(temp_y)
                ymax = max(temp_y)
             endif
          end
      'a': begin
              xmin = min(wave)
              xmax = max(wave)
              ymin = min(rawflux)
              ymax = max(rawflux)
              obj_plot_ls = where(wave ge xmin and wave le xmax)
           end
      '[': begin
             current_xrange = xmax - xmin
             xmin = (xmin - 0.75*current_xrange)(0) > min(wave)
             xmax = xmin + current_xrange
             obj_plot_ls = where(wave ge xmin and wave le xmax)
          end
      ']': begin
             current_xrange = xmax - xmin
             xmax = (xmax + 0.75*current_xrange)(0) < max(wave)
             xmin = xmax - current_xrange
             obj_plot_ls = where(wave ge xmin and wave le xmax)
           end
      '?': begin
             print,'   6: airmass - 0.1' 
             print,'   7: airmass - 0.01'
             print,'   8: airmass + 0.01'
             print,'   9: airmass - 0.1' 
             print,'   y: PWV - 0.1' 
             print,'   u: PWV - 0.01'
             print,'   i: PWV + 0.01'
             print,'   o: PWV - 0.1' 
             print,'   h: Delta_v - 1'
             print,'   j: Delta_v - 0.1'
             print,'   k: Delta_v + 0.1'
             print,'   l: Delta_v + 1'
             print,'   ,: FWHM - 0.5'
             print,'   .: FWHM + 0.5'
             print,'   b: Enter values manually'
             print,'   v: Click on plot for new x range'
             print,'   c: Click on plot for new y range'
             print,'   x: Enter new x range'
             print,'   z: Enter new y range'
             print,'   a: Reset display'
             print,'   [: Pan display left'
             print,'   ]: Pan display right' 
             print,'   q: Quit'
          end
      'q' : print,'Quitting'
      else: print,'Unknown key'
   endcase

endwhile

;;; Write parameters to FITS header
sxaddpar,hflux,'TELL_AIR',airmass,' Telluric Model Airmass'
sxaddpar,hflux,'TELL_PWV',pwv,' Telluric Model PWV'
sxaddpar,hflux,'TELL_DV',dv,' Telluric Model Velocity Offset (km/s)'
sxaddpar,hflux,'TELL_FWHM',ker_fwhm,' Telluric Model FWHM (km/s)'

;;; Save results
print,'Writing results...'
writefits,corrflux_file,corr_flux,hflux
if keyword_set(correrror_file) then begin
   sxaddpar,herror,'TELL_AIR',airmass,' Telluric Model Airmass'
   sxaddpar,herror,'TELL_PWV',pwv,' Telluric Model PWV'
   sxaddpar,herror,'TELL_DV',dv,' Telluric Model Velocity Offset (km/s)'
   sxaddpar,herror,'TELL_FWHM',ker_fwhm,' Telluric Model FWHM (km/s)'
   corr_error = rawerror / binned_model_trans
   writefits,correrror_file,corr_error,herror
endif

return
end