
pro go_align, t0, t1, xrt=xrt, level0=level0, search_array=search_array, $
  cat_sot=cat_sot, files_sot=files_sot, refresh_cat=refresh_cat, $
  interactive=interactive, do_save=do_save, do_restore=do_restore, $
  every_nth=every_nth, do_ct_adjust=do_ct_adjust, do_shimizu=do_shimizu, $
  data_align=data_align, t_off_max=t_off_max, do_stop=do_stop, $
  _extra=_extra

;+
; name: go_align_corr
; call: go_align, t0, t1
; -

if not exist(level0) and not keyword_set(quicklook) then level0 = 1
if not exist(mode_select) then mode_select = 'fg'
if not exist(dir_top) then dir_top = '/net/castor/Users/slater/data/sot_aia_align'
if not exist(t_off_max) then t_off_max = 60

interactive = 1

align_str = build_str_arr()

if keyword_set(interactive) then begin

  if not exist(cat_sot) or keyword_set(refresh_cat) then $
    sot_cat, t0, t1, cat_sot, files_sot, xrt=keyword_set(xrt), search_array=search_array, $
             level0=level0, quicklook=quicklook, count=count, tcount=tcount

  if keyword_set(xrt) then $
    mode_tags = 'ec_fw1_, ec_fw2_ naxis1, naxis2' else $
    mode_tags = 'wave, waveid, obs_type, naxis1, naxis2'

  ans = ''
  while strupcase(strmid(ans,0,1)) ne 'N' do begin

if not keyword_set(do_restore) then begin

    ss_mode = sot_umodes(incat=cat_sot, /int)

    files_sot_select = files_sot[ss_mode] ; [last_nelem(ss_mode,200)]
    n_files = n_elements(files_sot_select)

    ; Optionally only process every 'nth' file:
    if keyword_set(every_nth) then begin
      ss_nth = indgen(n_files/every_nth) * every_nth
      files_sot_select = files_sot_select[ss_nth]
    endif

    read_sot, files_sot_select, index_sot, data_sot
    data_align = data_sot*0

;   ss_0 = where_arr(index_sot.date_obs, align_str.date_obs, n_match0)
;   if n_match0 gt 0 then begin
;     ss_01 = where_arr(index_sot.date_obs, align_str.date_obs, /map)
;     ss_1 = ss_01[where(ss_01 ne -1, n_match1)]

; Find closest matches in correlation alignment db:
;   ss_close = tim2dset(align_str.date_obs, index_sot.date_obs, offset=t_off_sec)
    ss_close = tim2dset(anytim(align_str.date_obs, /ints), $
                        anytim(index_sot.date_obs, /ints), offset=t_off_sec)
; Filter out poor matches:
    ss_close_good = where(abs(t_off_sec) le t_off_max, n_close_good)
    if n_close_good gt 0 then begin
      ss_0 = ss_close_good
      ss_1 = ss_close[ss_close_good]                                             

      if tag_exist(index_sot,'percentd') then begin
mom_peak = moment(align_str[ss_1].peak_str, sdev=sdev_peak)
sig_from_mean = abs((align_str[ss_1].peak_str - mom_peak[0])/sdev_peak)

        ss_qual = where( (index_sot[ss_0].percentd eq 100) and $
                         (sig_from_mean le 2), n_qual)

        if n_qual gt 1 then begin
          ss_good0 = ss_0[ss_qual]
          ss_good1 = ss_1[ss_qual]
        endif else begin
          ss_good0 = ss_0
          ss_good1 = ss_1
        endelse
      endif else begin
        ss_good0 = ss_0
        ss_good1 = ss_1
      endelse
    endif else begin
      stop, ' No matches to alignment db times.  Stopping.'
    endelse

    index_sot = index_sot[ss_good0]
    data_sot = data_sot[*,*,ss_good0]

    ndim_data = size(data_sot, /n_dim)
    if ndim_data eq 4 then $
      data_sot = reform(data_sot[*,*,1,*])
;     sot_iquv_3d, index_sot, data_sot, index_sot, data_sot, /i_only ;, /noscale
;   data_sot = data_sot[*,*,ss_good0]

    align_str_match = align_str[ss_good1]
    n_files = n_elements(index_sot)
;   prstr, index_sot.date_obs + '    ' + align_str_match.date_obs

; Read scaling information for both SOT and SDO images:
    wave_aia = aia_wave_match(index_sot, ss_axis3=ss_axis3, ss_axis4=ss_axis4, $
                              loval_sot=loval_sot, hival_sot=hival_sot, $
                              do_invert_sot=do_invert_sot, nsig_sot=nsig_sot_tab, $
                              loval_sdo=loval_sdo, hival_sdo=hival_sdo, $
                              do_invert_sdo=do_invert_sdo, nsig_sdo=nsig_sdo_tab, $
                              _extra=_extra)

; Determine clipping and scaling for SOT image:
    if not exist(nsig_sot) then nsig_sot = nsig_sot_tab

;    if ( (loval_sot eq 0) and (hival_sot eq 0) ) then begin
;      print, 'No limit values in database.'
;      data_sot = bytscl(data_sot)
;    endif else begin
;      if (do_invert_sot or keyword_set(force_invert_sot)) then $
;        data_sot = 255b-bytscl(loval_sot>data_sot<hival_sot) else $
;        data_sot =      bytscl(loval_sot>data_sot<hival_sot)
;    endelse

    data_align = data_sot*0

    x_shift = (align_str_match.xdisp_rot)/(align_str_match.cdelt1)
    y_shift = (align_str_match.ydisp_rot)/(align_str_match.cdelt2)
    x0 = (index_sot.naxis1 + 1)/2. - x_shift
    y0 = (index_sot.naxis2 + 1)/2. - y_shift

    if keyword_set(do_save) then $
      save, n_files, index_sot, data_sot, data_align, align_str_match, x0, y0, $
        file = './ct_adjust_test.save'

endif else begin

    restore, './ct_adjust_test.save'
    if keyword_set(verbose) then print, 'Restored previously saved data.'

endelse    

; Optionally replace outlier x0 and y0 values with interpolated values:
    if keyword_set(do_ct_adjust) then begin
      ct_adjust, align_str_match, x0, y0, nsig_bad=nsig_bad, do_plot=do_plot, $
        _extra=_extra
    endif
;STOP
    for i=0, n_files-1 do begin
; IDL Call:
;   Result = ROT( A, Angle, [Mag, X0, Y0] [, /INTERP] [, CUBIC=value{-1 to 0}] [, MISSING=value] [, /PIVOT] )
;     data_align[0,0,i] = shift(data_sot[*,*,i], x_shift[i], y_shift[i])
      data_align[0,0,i] =   rot(data_sot[*,*,i], 0, 1, x0[i], y0[i], missing=0)
    endfor

; Optionally make a shimizu-aligned cube as well:
    if keyword_set(do_shimizu) then begin
      data_shimizu = data_sot*0

      t_sot = index_sot.date_obs
      sot_offsets = get_shimizu(t_sot)
      crval1_offset = reform(sot_offsets[0,*])
      crval2_offset = reform(sot_offsets[1,*])
      x_shift_shimizu = (crval1_offset)/(align_str_match.cdelt1)
      y_shift_shimizu = (crval2_offset)/(align_str_match.cdelt2)
      x0_shimizu = (index_sot.naxis1 + 1)/2. - x_shift_shimizu
      y0_shimizu = (index_sot.naxis2 + 1)/2. - y_shift_shimizu

      for i=0, n_files-1 do begin
        data_shimizu[0,0,i] = shift(data_sot[*,*,i], x_shift_shimizu[i], y_shift_shimizu[i])
        data_shimizu[0,0,i] =   rot(data_sot[*,*,i], 0, 1, x0_shimizu[i], y0_shimizu[i], missing=0)
      endfor

    endif

    naxis1 = index_sot[0].naxis1
    naxis2 = index_sot[0].naxis2

    if keyword_set(do_shimizu) then begin
      if naxis2 gt 512 then $
        cube = [[rebin(data_align,   naxis1/2, naxis2/2, n_files, /samp)], $
                [rebin(data_shimizu, naxis1/2, naxis2/2, n_files, /samp)]] else $
        cube = [[data_align], [data_shimizu]]
    endif else begin
      if naxis2 gt 512 then $
        cube = [[rebin(data_align, naxis1/2, naxis2/2, n_files, /samp)], $
                [rebin(data_sot,   naxis1/2, naxis2/2, n_files, /samp)]] else $
        cube = [[data_align], [data_sot]]
    endelse
    xstepper, cube

    if keyword_set(do_stop) then stop, 'MEASSAGE: GO_ALIGN: Stopping in GO_ALIGN.PRO.'

    read, ' Select another (def is yes)? ', ans
    
  endwhile

endif else begin

  t_grid = anytim(timegrid(t0, t1 , hours=1), /ccsds)
  n_grid = n_elements(t_grid)-1

  for i=0, n_grid-2 do begin
    sot_cat, t_grid[i], t_grid[i+1], cat_sot, files_sot, xrt=keyword_set(xrt), $
             search_array=search_array, level0=level0, quicklook=quicklook, $
             count=count, tcount=tcount

    ; Filter to select desired mode:
    mode = strlowcase(cat_sot.wave + ' ' + cat_sot.obs_type)
    ss_match = where(strpos(mode, mode_select) ne -1, n_match)

    if n_match gt 0 then begin

      files_sot = files_sot[ss_match]
      cat_sot = cat_sot[ss_match]

      ; Filter out non-existent files:
      ss_exist = where(file_exist(files_sot) eq 1, n_exist)

      if n_exist gt 0 then begin

        files_sot = files_sot[ss_exist]
        cat_sot = cat_sot[ss_exist]
        n_files = n_exist

        ; Optionally only process every 'nth' file:
        if keyword_set(every_nth) then begin
          ss_nth = indgen(n_files/every_nth) * every_nth
          files_sot = files_sot[ss_nth]
          cat_sot = cat_sot[ss_nth]
          n_files = n_elements(ss_nth)
        endif

        for j=0, n_files-1 do begin
          file_sot = files_sot[j]
          read_sot, file_sot, index_sot, data_sot

          sdo_corr, index_sot, data_sot, indexp_sot, datap_sot, $
            file_sot=file_sot, file_sdo=file_sdo, $
            norollcorrect=norollcorrect, nocrop2disk=nocrop2disk, $
            peak_str=peak_str, align_str=align_str, $
            do_img_plot=do_img_plot, do_hist_plot=do_hist_plot, $
            do_xstep=do_xstep, do_blink=do_blink, verbose=verbose, $
            _extra=_extra

          if exist(align_str) then begin
            if not exist(align_str_arr) then align_str_arr = align_str else $
              align_str_arr = concat_struct(align_str_arr, align_str)
            delvarx, align_str
          endif
        endfor

;       write_genxcat, align_str_arr, topdir=dir_top, /weeksub, /hours_round, $
;                       /nelements, /geny
        write_genxcat, align_str_arr, topdir=dir_top
        delvarx, align_str_arr

      endif

    endif

    print, ' Finished for '+ anytim(t_grid[i],   /yoh) + $
                             anytim(t_grid[i+1], /yoh)

  endfor

endelse

end

