function eis_orbit_spline, fitdata, setzero=setzero, plotvel=plotvel, $
                           spline=spline, error=error, yrange=yrange, $
                           psym=psym, quiet=quiet

;+
; NAME
;
;    EIS_ORBIT_SPLINE()
;
; PROJECT
;
;    Hinode/EIS
;
; EXPLANATION
;
;    Takes the centroid tag from the output of the EIS_AUTO_FIT
;    routine and computes a fit to the orbital variation. The velocity array
;    is averaged in the Y-direction (see the keyword YRANGE=), and a cubic
;    spline is fit to the averaged function (see the routine 
;    EIS_CUBIC_SPLINE). The spline fit is then subtracted from the centroid
;    array to give the corrected centroid array. The velocity array is also
;    corrected by converting the centroid to km/s.
;
;    The SPLINE= keyword can be used to save the spline fit so that it can
;    then be used for correcting other arrays. E.g., typically one would fit
;    the variation of the Fe XII 195 line, and then use this to correct all
;    other lines. Note that it seems most lines have the same centroid
;    variation in Angstroms (not km/s) over the orbit
;
;    The average orbit variation of the centroid is plotted to the screen,
;    with the fit overlaid as a thick line. If the spline= input is set,
;    then the fit from the spline structure will be overplotted for
;    comparison with the actual variation of the line.
;
; INPUTS
;
;    FITDATA  The Gaussian fit structure output by EIS_AUTO_FIT.
;
; OPTIONAL INPUTS
;
;    SPLINE   The fit to the orbit data from a previous run of 
;             EIS_ORBIT_SPLINE. For example, one may want to use the
;             Fe XII 195 fit for other, weaker emission lines.
;
;    YRANGE   By default the routine averages the centroid along the entire
;             field-of-view in Y before doing the spline fit. With this
;             keyword you can restrict the averaging to only the specified
;             pixel range. E.g., YRANGE=[0,19] means only Y-pixels 0 to 19
;             are used.
;
;    PSYM     Use PSYM= to specify a symbol for
;             plotting the averaged centroid data instead of the default
;             thin line.
;
; OPTIONAL OUTPUTS
;
;    SPLINE   If SPLINE does not exist when calling EIS_ORBIT_SPLINE, then
;             it will be created and contain the fit to the input structure
;             FITDATA. The tags of SPLINE are:
;
;              .ti   The temperatures at which the spline is defined. These
;                    are spaced at 5min intervals (note the ti are given in
;                    minutes).
;              .yi_vel  Not used, but tag kept for compatibility with earlier
;                       versions of routine.
;              .yi_cen  The spline points defined at the nodes ti. This is the
;                       spline that has been fit to the data.
;              .refwvl  The wavelength of the line used to derive the spline
;                       fit. If SPLINE is used as an input to eis_orbit_spline
;                       then refwvl will be the wavelength of the line
;                       originally used for the fit.
;              .cen_offset  If the /setzero keyword was set, then cen_offset
;                           is the offset applied to set the average centroid
;                           to zero. If /setzero not set, then cen_offset=0.
;              .vel_offset  If the /setzero keyword was set, then vel_offset
;                           is the offset applied to set the average velocity
;                           to zero. If /setzero not set, then vel_offset=0.
;
;             If you want to derive the centroid correction for arbitrary
;             times t, then you do the following in IDL:
;
;             IDL> ti=spline.ti
;             IDL> yi_vel=spline.yi_cen
;             IDL> y2=spl_init(ti,yi_cen)
;             IDL> y=spl_interp(ti,yi_cen,y2,t)
;
;             The Y are the centroids from the fit. To convert to velocity
;             do:
;
;             IDL> vel=lamb2v(y-refwvl,refwvl)
;
;             where refwvl is the reference wavelength of the line.
;
;    ERROR    The following values can be output:
;               0   everything's OK
;               1   the raster duration is less than 10 minutes
;               2   the end time is less than the start time (!)
;
; KEYWORDS
;
;    SETZERO  Normalize the centroid and velocity such that the average
;             centroid in the raster is set to zero km/s.
;
;    PLOTVEL  Does not do anything; kept for compatibility with older
;             versions. The routine will by default plot the orbit variation
;             - /quiet stops the plotting.
;
;    QUIET    If set, then the routine does not plot the orbital variation fit
;             or print information to the IDL window.
;
; OUTPUT
;
;    The output is a structure with the same tags as FITDATA, but the VEL
;    and CEN tags are corrected for the orbital variation.
;
; CALLS
;
;    EIS_CUBIC_SPLINE()
;
; HISTORY
;
;    Ver.1, 23-Feb-2007, Peter Young
;    Ver.2, 28-Feb-2007, Peter Young
;       corrected error when raster duration is less than 10 mins
;    Ver.3, 1-Mar-2007, Peter Young
;       corrected bug when there is no missing data
;    Ver.4, 11-Mar-2007, Peter Young
;       some data-sets have bad exposure start times. I now catch this problem.
;    Ver.5, 18-Jul-2007, Peter Young
;       added YRANGE= keyword
;    Ver.6, 8-Jan-2008, Peter Young
;       caught case where data-set only contains a single exposure
;    Ver.7, 5-Mar-2008, Peter Young
;       modified to deal correctly with missing data; also modified the
;       method for deriving the spline function - a call to eis_cubic_spline
;       is now used.
;    Ver.8, 11-Mar-2008, Peter Young
;       changed routine to fit centroid variation rather than velocity
;       variation; also fixed a bug introduced in ver. 7
;    Ver.9, 5-Nov-2008, Peter Young
;       the centroid offset (in angstroms) is added to the offset tag
;       in outdata; if the tag does not exist, then it is created.
;    Ver.10, 28-Oct-2009, Peter Young
;       corrected problem where very high time values (due to missing
;       data) messed up spline fit.
;-

IF n_params() LT 1 THEN BEGIN
  print,'Use:  IDL> fitdata2=eis_orbit_spline( fitdata [,/quiet, /setzero, '
  print,'                           yrange=, spline=, psym=, error= )'
  return,-1.
ENDIF


outdata=fitdata
missing=fitdata.missing

int=fitdata.int
vel=fitdata.vel
cen=fitdata.cen
tt=fitdata.exp_start_times

siz=size(vel)
nx=siz[1]
ny=siz[2]

IF nx EQ 1 THEN BEGIN
  print,'%EIS_ORBIT_SPLINE:  there is only one exposure in this data-set, so the orbit'
  print,'                    correction can not be performed. Returning...'
  return,-1
ENDIF


IF n_elements(yrange) EQ 0 THEN yrange=[0,ny-1]

t=tt/60.                        ; convert times to minutes
IF t[1] LT t[0] THEN BEGIN
  t_swtch=1
  t=reverse(t)
ENDIF ELSE BEGIN
  t_swtch=0
ENDELSE

error=0

;
; Do not proceed if the raster duration is less than 10 minutes
;
IF max(t)-min(t) LT 10 THEN BEGIN
  error=1
  return,outdata
ENDIF


;
; Collapse the 3D array along the slit (2nd dimension)
;
velx=average(vel[*,yrange[0]:yrange[1]],2,missing=missing)
cenx=average(cen[*,yrange[0]:yrange[1]],2,missing=missing)
;
; Check for T<1440 to filter out points with anomalously high time
; values.
;
i_notmiss=where(cenx NE missing AND t LE 1440.,n_notmiss)
IF n_notmiss NE 0 THEN BEGIN
  scenx=cenx[i_notmiss]
  tx=t[i_notmiss]
ENDIF ELSE BEGIN
  tx=t
ENDELSE


IF NOT keyword_set(quiet) THEN BEGIN
  plot,tx,scenx,/xsty,xra=[min(tx),max(tx)], /ynozero, $
       xtit='Time / minutes', $
       ytit='Centroid variation / '+string(197b), psym=psym, symsiz=1.5, $
       charsiz=1.4
ENDIF

IF n_tags(spline) EQ 0 THEN BEGIN
  y_spline=eis_cubic_spline(tx,scenx,nodes,yfit=yfit)
 ;
  IF NOT keyword_set(quiet) THEN BEGIN
    oplot,tx,yfit,th=2
  ENDIF
 ;
ENDIF ELSE BEGIN
   IF NOT keyword_set(quiet) THEN print,'% EIS_ORBIT_SPLINE: using SPLINE to correct for orbital variation'
  nodes=spline.ti
  y_spline=spline.yi_cen
 ;
  IF NOT keyword_set(quiet) THEN BEGIN
    y2=spl_init(nodes,y_spline)
    yfit=spl_interp(nodes,y_spline,y2,tx)
    oplot,tx,yfit-spline.refwvl+fitdata.refwvl,th=2
  ENDIF
ENDELSE



;
; Compute dcen and dv which are the orbital corrections to be applied to
; fitdata.cen and fitdata.vel, respectively.
; For dcen I need a reference wavelength to get the correction. This will be
; over-ridden later if /setzero is set.
; Note that I'm using T rather than TX now.
;
av_cenfit=average(yfit)
;
y2=spl_init(nodes,y_spline)
cenfit=spl_interp(nodes,y_spline,y2,t)
dcen=cenfit-av_cenfit
;
dv=lamb2v(dcen,fitdata.refwvl)


;
; make a note of missing data in fitdata 
;
i_miss=where(fitdata.cen EQ missing,nmiss)


;
; Correct the centroid and velocity for orbital variation
;
FOR i=0,nx-1 DO BEGIN
  IF t_swtch EQ 1 THEN j=nx-1-i ELSE j=i
  outdata.cen[j,*]=outdata.cen[j,*]-dcen[j]
  outdata.vel[j,*]=outdata.vel[j,*]-dv[j]
ENDFOR

IF NOT tag_exist(outdata,'offset') THEN outfit=add_tag(outdata,fltarr(nx,ny),'offset')
ident_y=fltarr(ny)+1.
outdata.offset=outdata.offset+dcen#ident_y

;
; Set any missing values back to missing
;
IF nmiss NE 0 THEN BEGIN
  outdata.cen[i_miss]=missing
  outdata.vel[i_miss]=missing
ENDIF


IF keyword_set(setzero) THEN BEGIN
  av_cen=average(outdata.cen,missing=missing)
  cen_offset=av_cen-fitdata.refwvl
  outdata.cen=outdata.cen-cen_offset
  outdata.offset=outdata.offset+cen_offset
 ;
  vel_offset=average(outdata.vel,missing=missing)
  outdata.vel=outdata.vel-vel_offset
 ;
  IF nmiss NE 0 THEN BEGIN
    outdata.cen[i_miss]=missing
    outdata.vel[i_miss]=missing
  ENDIF
ENDIF ELSE BEGIN
  cen_offset=0.
  vel_offset=0.
ENDELSE


IF n_tags(spline) EQ 0 THEN BEGIN
  ns=n_elements(nodes)
  spline={ ti: nodes, $
           yi_cen: y_spline, $
           yi_vel: fltarr(ns), $
           cen_offset: cen_offset, $
           vel_offset: vel_offset, $
           refwvl: fitdata.refwvl }
ENDIF


return,outdata

END