
function sot_dc_align, index_in, files_sot=files_sot, use_ts_data=use_ts_data, use_gls_data=use_gls_data, $
  off_sot2sc_ew=off_sot2sc_ew, t_sot2sc_ew=t_sot2sc_ew, $
  off_sot2sc_ns=off_sot2sc_ns, t_sot2sc_ns=t_sot2sc_ns, $
  offsets_xcen=offsets_xcen, offsets_ycen=offsets_ycen, $
  use_spline=use_spline, use_linear=use_linear, quiet=quiet, qstop=qstop

;+
;   Name: SOT_DC_ALIGN
;
;   Purpose: Correct XCEN and YCEN tags of SOT index records using Shimizu-generated Hinode alignment results.
;            Return either in full SOT index record array (if input index array passed in), or in reduced tag
;            index array (DATE_OBS, XCEN, YCEN) if SOT file list passed in.
;
;   Input Parameters:
;      INDEX_IN (Positional) - If passed replicate the index array and update XCEN and YCEN tags in the new array.
;                              If not passed, then it is assumed that keyword parameter FILES_SOT is passed.
;                              In this case, use MREADFITS_HEADER to read a reduced set of tags from all headers
;                              (DATE_OBS, XCEN, YCEN) and update the XCEN and YCEN tags in this index array.
;   Output Paramters:
;      Index array (either full or reduced) with XCEN and YCEN tags corrected,
;
;   Keyword Parameters:
;      FILES_SOT - See above.
;
;   Calling Sequence:
;      IDL> index_out = sot_dc_align(index_in [,files_sot=files_sot])
;
;   History:
;      15-mar-2009 - G.L.Slater
;
;-

; Default is to use Shimizu offset data:
if not keyword_set(use_gls_data) then use_ts_data = 1 else use_ts_data = 0
; Default is spline fit:
use_linear = 1

; If SOT FITS file list passed instead of index array, then read headers from files.
; We use mreadfits_header.pro because it's faster.  However, for as yet undetermined reasons, its speed
; drops off dramaticially when n_files_sot gets large.  Therefore, we arbitrarily 'chunk' the headers in
; 100 file reads:

if not exist(index_in) then begin
  if exist(files_sot) then begin
    n_files_sot = n_elements(files_sot)
    n_reads = long(ceil(float(n_files_sot)/100.))
    tag_list = ['date_obs,xcen,ycen']

    for i=0l, n_reads-1 do begin
      if i eq (n_reads-1) then $
        mreadfits_header, files_sot[(100l*i):*], index0, only_tags=tag_list else $
        mreadfits_header, files_sot[(100l*i):(100l*(i+1l)-1)], index0, only_tags=tag_list
      if i eq 0 then index_in = index0 else index_in = concat_struct(index_in, index0)
    endfor
  endif
endif

index_out = index_in
n_index = n_elements(index_out)

t_sec_sot = anytim(index_out.date_obs)
t_int_sot = anytim(index_out.date_obs, /int)

; Read in either Shimizu offset data or Slater offset data if not passed in:

if ( (not exist(off_sot2sc_ew)) or (not exist(t_sot2sc_ew)) or $
     (not exist(off_sot2sc_ns)) or (not exist(t_sot2sc_ns)) ) then begin

  print, ' Offset data not passed.  Reading it in."

  if not keyword_set(use_gls_data) then begin
; Read in Shimizu offset data:
    rd_align_data, sot_fovcen2scpnt_ew=off_sot2sc_ew, t_sot_fovcen2scpnt_ew=t_sot2sc_ew, $
                   sot_fovcen2scpnt_ns=off_sot2sc_ns, t_sot_fovcen2scpnt_ns=t_sot2sc_ns
  endif else begin
; Read in Slater offset data:
    dum = tim2sc_offsets( sc2sot_ew_off=off_sot2sc_ew, t_ew_sot=t_sot2sc_ew, $
                          sc2sot_ns_off=off_sot2sc_ns, t_ns_sot=t_sot2sc_ns)
  endelse

endif

; Define interpolation regions for either Shimizu or slater offset data:
if not keyword_set(use_gls_data) then $
  sttim_region = ['01-oct-2006','23-apr-2007','12-aug-2007','06-may-2008','06-aug-2008','01-jan-2020'] else $
  sttim_region = ['01-oct-2006','23-apr-2007','12-aug-2007','06-may-2008','06-aug-2008','01-jan-2020']
n_region = n_elements(sttim_region)

sttim_reg_sec = anytim(sttim_region)
; TODO - Remove this line:
; ss_sttim_reg = tim2dset(anytim(sttim_region,/int))
; Convert offset table times to seconds:
t_sec_sot2sc_ew = anytim(t_sot2sc_ew)
t_sec_sot2sc_ns = anytim(t_sot2sc_ns)

; Do piecewise interpolation against the Shimizu offset data:

; Define arrays to contain offset values:
offsets_xcen = fltarr(n_index)
offsets_ycen = fltarr(n_index)

for i=0, n_region-2 do begin
  ss_sot_reg0 = where( ((t_sec_sot ge sttim_reg_sec[i]) and (t_sec_sot lt sttim_reg_sec[i+1])), n_match_sot_reg0 )
  if (n_match_sot_reg0 gt 0) then begin
    ss_off_ew_reg0 = where( ((t_sec_sot2sc_ew ge sttim_reg_sec[i]) and (t_sec_sot2sc_ew lt sttim_reg_sec[i+1])), n_match_off_ew_reg0 )
    ss_off_ns_reg0 = where( ((t_sec_sot2sc_ns ge sttim_reg_sec[i]) and (t_sec_sot2sc_ns lt sttim_reg_sec[i+1])), n_match_off_ew_reg0 )

    t_sec_sot_reg0 = t_sec_sot[ss_sot_reg0]
    t_sec_sot2sc_ew_reg0 = t_sec_sot2sc_ew[ss_off_ew_reg0]
    t_sec_sot2sc_ns_reg0 = t_sec_sot2sc_ns[ss_off_ns_reg0]
    off_sot2sc_ew_reg0 = off_sot2sc_ew[ss_off_ew_reg0]
    off_sot2sc_ns_reg0 = off_sot2sc_ns[ss_off_ns_reg0]

    if keyword_set(use_spline) then begin
      off_sot_ew = spline(t_sec_sot2sc_ew_reg0, off_sot2sc_ew_reg0, t_sec_sot_reg0)
      off_sot_ns = spline(t_sec_sot2sc_ns_reg0, off_sot2sc_ns_reg0, t_sec_sot_reg0)
    endif

    if keyword_set(use_linear) then begin
      off_sot_ew = interpol(off_sot2sc_ew_reg0, t_sec_sot2sc_ew_reg0, t_sec_sot_reg0)
      off_sot_ns = interpol(off_sot2sc_ns_reg0, t_sec_sot2sc_ns_reg0, t_sec_sot_reg0)
    endif

    offsets_xcen[ss_sot_reg0] = off_sot_ew
    offsets_ycen[ss_sot_reg0] = off_sot_ns
  endif

endfor

; Make adjustsments to SOT xcen, ycen:

index_out.xcen = index_out.xcen + offsets_xcen
index_out.ycen = index_out.ycen + offsets_ycen

if keyword_set(qstop) then stop, ' Stopping on request.'

return, index_out

end
