PRO eis_orbit_correction, data, wvl, fitdata=fitdata, $
                          yrange=yrange, plot=plot, iwin=iwin, $
                          noprep=noprep

;+
; NAME
;
;    EIS_ORBIT_CORRECTION
;
; PROJECT
;
;    Hinode/EIS
;
; EXPLANATION
;
;    Calculates the wavelength vector for each exposure in the EIS data-set,
;    corrected for the orbital variation. The method is as follows.
;
;    The data are averaged in the Y-direction to yield a single line profile
;    for each exposure. The centroid of this profile is determined with the
;    routine LINE_CENTROID. We then have the centroid defined as a function
;    of time. This curve is then smoothed over 5 positions with the IDL SMOOTH
;    function.
;
;    A spline function with nodes spaced at 5 minute intervals is then fit
;    to this function. This spline is then used to re-derive the centroid
;    for each exposure. The difference between this centroid and the average
;    centroid is then used to correct the wavelength vector for the window
;    at each exposure. The resulting 2D array of wavelengths is output.
;
;    The offsets calculated for one window can be applied to another window
;    using the FITDATA optional input/output (see below). Thus a strong line
;    (e.g., Fe XII 195) can be used to correct the wavelengths for a weaker
;    line. Note that no attempt is made to scale the offset by the ratio of
;    the wavelengths, which seems to be consistent with observations.
;
; INPUTS
;
;    DATA    An EIS data object.
;
; OPTIONAL INPUTS
;
;    IWIN    The index of the wavelength window to be analysed. If not set,
;            then the user is asked to specify a window number on the command
;            line.
;
;    FITDATA A structure (see below for definition) that contains orbital
;            wavelength offsets determined from a previous run of
;            EIS_ORBIT_CORRECTION.
;
;    YRANGE  By default the routine averages all pixels in the Y direction
;            to determine the average profile for an exposure. By setting
;            YRANGE, a subset of the Y-pixels can be specified. E.g.,
;            YRANGE=[50,100] means pixels 50 to 100 will be averaged.
;
; OPTIONAL OUTPUTS
;
;    FITDATA A structure that contains the derived orbital wavelength offsets.
;            The tags are:
;              .t  The time in minutes for each exposure, specified relative
;                  to the start of the first exposure.
;              .offset  The offset in angstroms of the exposure wavelengths
;                       from the reference wavelength.
;              .refwvl  The reference wavelength in angstroms, defined to be
;                       the average centroid position over the raster.
;              .iwin    The index of the window used to derive the offsets.
;
; KEYWORDS
;
;    PLOT   If set, then the smoothed centroid position is plotted as a
;           function of time, with the spline fit overplotted with a thicker
;           line. If FITDATA has been set, then the plot compares the
;           calculated centroid position with the fit from the reference
;           window.
;
;    NOPREP In the standard mode of operation, this routine assumes the EIS
;           data have been processed with EIS_PREP. If you have not run
;           eis_prep and you want a "quick-and-dirty" solution, then set
;           this keyword. Note that the data array will have the pedestal
;           removed by the routine EIS_DEBIAS, since otherwise the centroid
;           is not accurately derived.
;
; OUTPUTS
;
;    WVL     A 2-dimensional array of size (NL,NX) where NL is the number
;            of wavelength pixels for the wavelength window, and NX is the
;            number of exposures. Each column WVL[*,i] contains the wavelength
;            vector for exposure i, corrected for the orbital variation.
;
; CALLS
;
;    LINE_CENTROID, EIS_EXPOSURE_TIMES, EIS_PIXEL_TO_WAVELENGTH
;    EIS_DEBIAS
;
; HISTORY
;
;    Ver.1, 25-Apr-2007, Peter Young
;-

id=data->getline_id()
IF n_elements(iwin) EQ 0 THEN BEGIN
  i=where(id NE '')
  nwind=n_elements(id[i])

  ; -----------------------------------<
  ; Choose which wavelength window to process
  ;
  PRINT,''
  PRINT,'Choose a wavelength window to process:'
  PRINT,''
 ;
  FOR i=0,nwind-1 DO PRINT,i+1,'.  ',trim(id[i])
 ;
  READ,choice
  iwin=FIX(choice-1)
 ; -----------------------------------<
ENDIF


wd=data->getvar(iwin)
xs=data->getxs()   ; detX start positions of windows
xs=xs[iwin]
xw=data->getxw()   ; width of each window
xw=xw[iwin]
slit_ind=data->getslit_ind(0)   ; get slit index


;
; Find missing data values, and set the values in wd to the value of
; MISSING.
;
IF NOT keyword_set(noprep) THEN BEGIN
  wde=data->geterr(iwin)
  missing=-100.
  i=where(wde EQ 0.,ni)         ; this is definition of missing data
  IF ni NE 0 THEN wd[i]=missing
ENDIF


;
; TT are the start times of each exposure relative to the start time of the
; first exposure
;
eis_exposure_times,data,tt

IF max(tt)-min(tt) LT 600. THEN BEGIN
  print,'% EIS_ORBIT_CORRECTION: the raster duration is too short to derive an accurate'
  print,'                        orbital correction. Returning...'
ENDIF


;
; Collapse the 3D array along the slit (2nd dimension)
;
siz=size(wd)
ny=siz[2]
IF n_elements(yrange) EQ 0 THEN yrange=[0,ny-1]
IF keyword_set(noprep) THEN BEGIN
  eis_debias,wd
  wdx=average(wd[*,yrange[0]:yrange[1],*],2)
ENDIF ELSE BEGIN
  wdx=average(wd[*,yrange[0]:yrange[1],*],2,missing=missing)
ENDELSE

siz=size(wdx)
nx=siz[2]
nl=siz[1]

x=findgen(xw)+xs    ; 

cen=fltarr(nx)

;
; Derive line centroid for each exposure, and store in CEN. The centroid is
; given in pixels, not wavelength.
;
FOR i=0,nx-1 DO BEGIN
  y=wdx[*,i]
  line_centroid,x,y,cent,/silent
  cen[i]=cent
ENDFOR 

;
; Convert pixel centroid to wavelength
;
IF xs GE 2198 THEN BEGIN
  xoff=2198 
  ccd=2
ENDIF ELSE BEGIN
  xoff=50
  ccd=1
ENDELSE
;
wcen=fltarr(nx)
FOR i=0,nx-1 DO BEGIN
  pix=cen[i]-xoff
  wcen[i]=eis_pixel_to_wavelength(pix,ccd,slit_ind)
ENDFOR


;
; Smooth the wavelength centroid
;
scen=smooth(wcen,5)


t=tt/60.    ; convert times to minutes



IF keyword_set(plot) THEN BEGIN
  plot,t,scen,/xsty,/ynozero, $
       xtit='Time / minutes', $
       ytit='Centroid / '+string(197b)
ENDIF

;
; The average centroid position is considered to be the reference wavelength
; (REFWVL). SCEN is re-defined by subtracting this value.
;
refwvl=average(scen)
scen=scen-refwvl


;
; The following code calculates the orbital offset in angstroms between the
; centroid for each exposure and the reference wavelength (REFWVL).
; If FITDATA has been specified, then the offsets contained in it are used.
; Otherwise the offsets are calculated by fitting a spline.
;
IF n_tags(fitdata) EQ 0 THEN BEGIN
 ;
 ; define a spline (ti,yi) with nodes at every 5 minutes
 ;
  t_max=fix(max(t)/5)*5
  ti=findgen(t_max/5+1)*5
 ;
  y2=spl_init(t,scen)
  yi=spl_interp(t,scen,y2,ti)
 ;
  y2=spl_init(ti,yi)  
 ;
  offset=spl_interp(ti,yi,y2,t)
 ;
  fitdata={ t: t, $
            offset: offset, $
            iwin: iwin, $
            refwvl: refwvl}
ENDIF ELSE BEGIN
  print,'% EIS_ORBIT_CORRECTION: using correction contained in FITDATA input.'
  offset=fitdata.offset
ENDELSE


IF keyword_set(plot) THEN oplot,t,offset+refwvl,th=2

wvl=fltarr(nl,nx)
wvl0=eis_pixel_to_wavelength(x-xoff,ccd,slit_ind)
;
FOR i=0,nx-1 DO wvl[*,i]=wvl0-offset[i]


END