;+
; 
;             NAME : eis_wave_corr_old
; 
;          PURPOSE : PRELIMINARY EIS wavelength corrections for orbital variations and
;                    slit tilt using moments. *** THIS IS NOW THE LEGACY VERSION ***
;
; CALLING SEQUENCE : eis_wave_corr_old,file,wvl_corr,dw_tilt,dw_t
; 
;           INPUTS : EIS level 1 file name
; 
;  OPTIONAL INPUTS : ymin, ymax - select a subregion along the slit
;                    for averaging.
; 
;         KEYWORDS : show - display the results on the screen, the
;                    default is to write the plots to the zbuffer. In
;                    either case jpeg images are written to the
;                    current directory.
;
;                  : nospline - do not do smoothing and spline fitting
;                    to the wavelength variation with time. The spline
;                    can only be done if the duration of the raster is
;                    longer than 300s.
;
;                  : do_fit - use a Gaussian fit to the line profile
;                    instead of a moment calculation. This takes MUCH
;                    longer but any residual wavelength correction is
;                    likely to be smaller.
; 
;          OUTPUTS : wvl_corr - the wavelength correction for each pixel
;
;                  : dw_tilt - the slit tilt component
;
;                  : dw_t - the temporal component
;
;                  : also saves these arrays to a genx file for future
;                    use.
;
; OPTIONAL OUTPUTS : none
;
;    COMMON BLOCKS : none
; 
;          EXAMPLE : IDL> eis_wave_corr_old,file,wvl_corr,dw_tilt,dw_t,/show
;                    IDL> d.wvl = d.wvl - wvl_corr[i,j]
;
;    SPECIAL CALLS : eis_slit_tilt, eis_setup_xwindow
;
;          WRITTEN : HPW : 17-MAY-2007
;                    HPW : 10-OCT-2007 : Added 2" slit, changed layout
;                    HPW : 27-DEC-2007 : Added spline and show keyword
;                    HPW : 09-JAN-2008 : Added mpfit option
;                    HPW : 23-JAN-2008 : Fixed SNS bug
;                    HPW : 16-SEP-2008 : Fixed to work with missing data
;                    HPW : 25-SEP-2008 : Modified handling of bad data 
;                    HPW : 07-NOV-2008 : Fixed to work with NEXP_PRP eq 2
;                    HPW : 29-DEC-2008 : Fixed file name bug
;
;-

;; **********************************************************************

function eis_wave_corr_moment,x,y,ints=ints

  idx  = sort(y)
  back = median(y[idx[0:4]])
  prof = y - back
  dx   = x[1]-x[0]

  ints  = total(prof)*dx
  peak  = max(prof)
  cent  = total(x*prof)*dx/ints
  width = ints/(peak*sqrt(2*!pi))

return,[peak,cent,width,back]
end

;; **********************************************************************

function eis_wave_corr_fit,x,y,e,ints=ints

  yfit  = mpfitpeak(x,y,fit,error=e)

  peak  = fit[0]
  cent  = fit[1]
  width = fit[2]
  back  = fit[3]
  ints  = sqrt(2*!pi)*fit[0]*fit[2]

return,[peak,cent,width,back]
end

;; **********************************************************************

pro eis_wave_corr_old,file,wvl_corr,dw_tilt,dw_t,ymin=ymin,ymax=ymax,$
                      show=show,$
                      nospline=nospline,$
                      do_fits=do_fits

  if keyword_set(do_fits) then begin
    message,'USING GAUSSIAN FITS . . . THIS MAY TAKE A VERY LONG TIME',/info
  endif

  ;; --- save some of the existing display parameters
  p_old = !p
  device_name_old = !d.name
  tvlct,r_old,g_old,b_old,/get

  ;; --- assumed rest wavelegnth for Fe XII 195
  wave0 = 195.12

  ;; --- read the eis data file 
  break_file,file,disk, dir, filnam, ext, ver, node,/last

  if ext eq '.genx' then begin
    restgen,file=file,struct=d
    dwave=str2number(d.line_id)
    if dwave ne wave0 then message,'WARNING. LINE_ID: '+d.line_id+'. EXPECTING 195.12',/info
  endif else begin
    iwin = eis_get_iwin(file,wave0)
    if iwin lt 0 then begin
    	message,'Fe XII 195.12 not found, returning . . .',/info
    	return
    endif 
    d = eis_getwindata(file,wave0)
  endelse

  ;; --- deal with multiple exposures at each position
  if d.hdr[0].nexp_prp eq 2 then begin

    ;;
    ;; To deal with this type of data we need to rearrange things so that
    ;; the exposures are in time order. The wvl_corr array will be unpacked 
    ;; at the end.
    ;;

    nx_new = d.hdr[0].nexp_prp*d.nx

    if d.hdr.nraster eq 1 then begin
      i0 = 1
      i1 = 0
    endif else begin
      i0 = 0
      i1 = 1
    endelse

    ;; -- spectra and error
    int_new = fltarr(d.nl,nx_new,d.ny)
    err_new = fltarr(d.nl,nx_new,d.ny)
    idx1 = 2*findgen(d.nx)+1
    idx2 = 2*findgen(d.nx)
    int_new[*,idx1,*] = d.int[*,*,*,i0]
    int_new[*,idx2,*] = d.int[*,*,*,i1]
    err_new[*,idx1,*] = d.err[*,*,*,i0]
    err_new[*,idx2,*] = d.err[*,*,*,i1]

    ;; --- time
    time_ccsds_new = strarr(nx_new)
    time_new = fltarr(nx_new)
    time_new[idx1] = d.time[*,i0]
    time_new[idx2] = d.time[*,i1]

    time_ccsds_new[idx1] = d.time_ccsds[*,i0]
    time_ccsds_new[idx2] = d.time_ccsds[*,i1]

    ;; --- data quality
    if tag_exist(d,'data_quality') then begin
      data_quality_new = bytarr(nx_new)
      data_quality_new[idx1] = d.data_quality[*,i0]
      data_quality_new[idx2] = d.data_quality[*,i1]
    endif
    
    d = rem_tag(d,['nx','int','err','time','time_ccsds'])
    
    d = add_tag(d,nx_new,'nx')
    d = add_tag(d,int_new,'int')
    d = add_tag(d,err_new,'err')
    d = add_tag(d,time_new,'time')
    d = add_tag(d,time_ccsds_new,'time_ccsds')
    if tag_exist(d,'data_quality') then begin
      d = rem_tag(d,'data_quality')
      d = add_tag(d,data_quality_new,'data_quality')
    endif

  endif

  ;; --- create output file name
  break_file,file,disk,dir,name,ext,/last
  opf = strmid(name,0,22)+'.wave'

  ;; --- set the yrange for averaging
  if not(keyword_set(ymin)) then ymin = 0 else ymin = ymin > 0
  if not(keyword_set(ymax)) then ymax = d.ny-1 else ymax = ymax < (d.ny-1)

  ;; ----------------------------------------------------------------------
  ;; --- compute slit tilt

  dw_tilt = eis_slit_tilt(d)

  ;; ----------------------------------------------------------------------
  ;; --- velocity oscillation

  ;; --- wavelegnth range for moments
  w1 = 194.9
  w2 = 195.3
  mm = where(d.wvl ge w1 and d.wvl le w2)

  ;; --- compute moments
  cent = fltarr(d.nx,d.ny)
  ints = fltarr(d.nx,d.ny)
  nfit = float(n_elements(cent))
  cntr = 0L
  for ii=0,d.nx-1 do begin
    for jj=0,d.ny-1 do begin

      if keyword_set(do_fits) then begin
        fit = eis_wave_corr_fit(d.wvl[mm],d.int[mm,ii,jj],$
                                d.err[mm,ii,jj],ints=intensity)
        percent = 100.*cntr/nfit
        if cntr eq 0 then reset = 1 else reset = 0
        progress,percent,frequency=2.,reset=reset
        cntr = cntr + 1
      endif else begin
        fit = eis_wave_corr_moment(d.wvl[mm],d.int[mm,ii,jj],$
                                   ints=intensity)
      endelse

      cent[ii,jj] = fit[1]
      ints[ii,jj] = intensity
    endfor
  endfor
  if keyword_set(do_fit) then progress,100.,/last

  ;; --- average centroids along Y
  dw_t = fltarr(d.nx)
  for ii=0,d.nx-1 do begin
    dw_t[ii] = average(cent[ii,ymin:ymax]-dw_tilt[ymin:ymax]) - wave0
  endfor

  ;; --- if possible, do a spline to the velocities 

  delta_t  = 300. ;; spacing for the knots in the spline
  n_smooth = 3    ;; smoothing parameter

  dt = max(d.time)-min(d.time)
  if (dt le delta_t) and not(keyword_set(nospline)) then begin
    message,'NOT ENOUGH DATA FOR SPLINE FIT',/info
    message,'TURING SPLINE FIT OFF',/info
    nospline = 1
  endif

  ;; --- locate missing data
  if tag_exist(d,'data_quality') then begin
    good = where(d.data_quality eq 1,ngood)
    bad  = where(d.data_quality eq 0,nbad)
  endif else begin
    ;; --- data_quality is not defined, improvise!
    ngood = n_elements(d.time) 
    good  = indgen(ngood)
    nbad  = 0
  endelse

  MAX_BAD = 0.25 ;; 
  if float(nbad)/float(ngood) gt MAX_BAD then begin
    message,'NOT ENOUGH GOOD DATA FOR SPLINE FIT',/info
    message,'TURING SPLINE FIT OFF',/info
    nospline = 1
  endif

  ;; --- default is to use the Z buffer for the display
  if not(keyword_set(show)) then zbuff = 1 else zbuff = 0

  if not(keyword_set(nospline)) then begin

    time = d.time

    time_good = time[good]
    dw_t_good = dw_t[good]
      
    ;; --- spline fit the original data "delta_t" intervals 
    ;;
    ;; Note that rasters and SNS are ordered differently. Rasters are
    ;; stepped from east to west and so the last exposure is in the
    ;; first column. Thus they need to be reversed. SNS don't need to
    ;; be reversed. HPW: 23-JAN-2008.
    ;;
    time_diff = time_good[1] - time_good[0]
    if time_diff lt 0 then begin
      x = reverse(time_good-time_good[ngood-1])
      y = reverse(dw_t_good)
    endif else begin
      x = time_good - time_good[0]
      y = dw_t_good
    endelse

    y  = smooth(y,n_smooth)
    ns = (max(time_good)-min(time_good))/delta_t+1
    xs = findgen(ns)*delta_t
    ys = spline(x,y,xs)

    ;; --- spline back to the orginal time grid
    dw_t_new = spline(xs,ys,x)
    if time_diff lt 0 then dw_t_new = reverse(dw_t_new)
    
    ;; --- display
    eis_setup_xwindow,900,450,0,zbuff=zbuff
    linecolors
    utplot,d.time_ccsds[good],dw_t_good,xstyle=1,thick=3,psym=1
    outplot,d.time_ccsds[good],dw_t_new,color=2,thick=2

    legend,['original','spline fit'],linestyle=[0,0],color=[255,2],$
        box=0,spacing=1.5,pspacing=1.5
    write_jpeg,opf+'.a.jpg',tvrd(true=1),true=1,quality=95
    
    dw_t[good] = dw_t_new
    if nbad gt 0 then dw_t[bad] = 0.0

  endif else begin

    if nbad gt 0 then dw_t[bad] = 0.0

    eis_setup_xwindow,900,450,0,zbuff=zbuff
    linecolors
    utplot,d.time_ccsds,dw_t,xstyle=1,thick=2
    legend,['original'],linestyle=[0],color=[255],$
        box=0,spacing=1.5,pspacing=1.5
    write_jpeg,opf+'.a.jpg',tvrd(true=1),true=1,quality=95

  endelse

  ;; ----------------------------------------------------------------------
  ;; --- correction array
  
  wvl_corr = fltarr(d.nx,d.ny)
  for i=0,d.nx-1 do wvl_corr[i,*] = dw_tilt
  for j=0,d.ny-1 do wvl_corr[*,j] = wvl_corr[*,j] + dw_t

  ;; ----------------------------------------------------------------------
  ;; --- display

  ;; --- aspect ratio of raster
  eis_raster_fov,d,scale,origin,fovx,fovy
  aspect = fovy/fovx
  
  ;; --- setup the display don't let the window get too large
  IMGnx = 300  ;; each window is (300)x(300*aspect) pixels
  IMGny = round(IMGnx*aspect)
  max_win_ny = 900 ;; maximum vertical extent of window
  if IMGny gt max_win_ny then eps = max_win_ny/float(IMGny) else eps = 1.0
  eis_setup_xwindow,round(eps*IMGnx*4),round(eps*IMGny),1,zbuff=zbuff
  !p.multi=[0,4,1,0,0]
  !p.charsize=2
    
  ;; --- display intensity
  eis_colors,/intensity
  temp  = sigrange(ints,f=0.999,range=range)
  range = range[1]*[1.0E-2,1]
  image = eis_scale_image(ints,'logarithmic',range)
  plot_image,image,title='Fe XII Intensity',scale=scale,origin=origin
  ;; --- 
  x1 = d.solar_x[0]
  x2 = x1 + d.nx*d.scale[0]
  plots,[x1,x2],d.solar_y[ymin],linestyle=1
  plots,[x1,x2],d.solar_y[ymax],linestyle=1
  
  ;; --- display velocity correction
  eis_colors,/velocity
  plot_image,wvl_corr,title='Wavelength Correction',scale=scale,origin=origin
  
  ;; --- display uncorrected velocities
  vel   = 2.99792e+5*(cent-wave0)/wave0
  image = bytscl(vel,-20,20,top=254)
  plot_image,image,title='Uncorrected Velocities',/noscale,$
      scale=scale,origin=origin
  
  ;; --- display corrected velocities  
  cent  = cent - wvl_corr
  vel   = 2.99792e+5*(cent-wave0)/wave0
  image = bytscl(vel,-20,20,top=254)
  plot_image,image,title='Corrected Velocities',/noscale,$
      scale=scale,origin=origin
  
  write_jpeg,opf+'.b.jpg',tvrd(true=1),true=1,quality=95
    
  ;; ----------------------------------------------------------------------
  ;; --- save

  if keyword_set(do_fit) then method = 'MPFIT' else method = 'MOMENT'

  ;; --- unpack the wvl_corr array if necessary
  if d.hdr[0].nexp_prp gt 1 then begin
    nx_new = d.nx/d.hdr[0].nexp_prp
    help,wvl_corr
    wvl_corr_new = fltarr(nx_new,d.ny,d.hdr[0].nexp_prp)
    wvl_corr_new[*,*,0] = wvl_corr[idx1,*]
    wvl_corr_new[*,*,1] = wvl_corr[idx2,*]
    wvl_corr = wvl_corr_new
  endif

  text = [" ---------- EIS_WAVE_CORR OUTPUT ---------- ",$
          " computed with "+method,$
          " ymin (pixels) = "+trim(ymin),$
          " ymax (pixels) = "+trim(ymax),$
          " file = "+opf+".genx",$
          " read with ",$
          "  IDL> restgen,wvl_corr,dw_tilt,dw_t,file='"+opf+".genx'",$
          "  wvl_corr = 2D wavelength correction array",$
          "  dw_t = 1D time correction",$
          "  dw_tilt = 1D slit tilt correction",$
          " use with ",$
          "  d.wvl = d.wvl - wvl_corr[i,j]",$
          " ----------------------------------------- "]
  if keyword_set(show) then hprint,text
  savegen,wvl_corr,dw_tilt,dw_t,file=opf+'.genx',text=text

  ;; --- clean up the display
  !p = p_old
  tvlct,r_old,g_old,b_old
  set_plot,device_name_old
  
return
end