
pro mk_corr_cat, t0, t1, days_chunk=days_chunk, files_corr_db=files_corr_db, $
  file_driven=file_driven, time_driven=time_driven, $
  dir_top=dir_top, search_array=search_array, refresh=refresh, $
  cat_sot=cat_sot, cat_corr=cat_corr, use_ctme=use_ctme, $
  p_degree=p_degree, spline_sigma=spline_sigma, $
  med_wid=med_wid, $
  diff_jump_corr=diff_jump_corr, diff_jump_ct=diff_jump_ct, $
  do_plot=do_plot, do_all_seg_plot=do_all_seg_plot, do_save=do_save, $
  charsize=charsize, $
  loud=loud, quiet_where=quiet_where, qstop=qstop, do_test=do_test, _extra=_extra

;+
; PHILOSOPHY:
;   possible criteria for adding a 'boundary':
;     - jump in unsmoothed, un-median-filtered correlation offset time series in either
;       x or y axis exceeds some value. That value can be:
;         - an absolute difference between temporaily adjacent offsets (e.g., 1 arcsec).
;         - some number of sigma of the smoothed and/or median-filtered correlation
;           offset time series.
;     - time difference between adjacent correlation offsets exceeds a maximum time difference
;     - a change in the obs_type of the image (presumably (possibly) associated with a
;       change in an observing mode).
;-

if not exist(dir_top) then begin
  if keyword_set(do_test) then $
    dir_top = '$SSWDB/hinode/gen/hinode_sdo_alignment/sot_aia_corr_db_test' else $
    dir_top = '/sanhome/slater/public_html/outgoing/slater/sot_corr_testing'
;   dir_top = '$SSWDB/hinode/gen/hinode_sdo_alignment/sot_aia_corr_db'
endif
if not exist(dir_corr_segments) then begin
  if keyword_set(do_test) then $
  dir_corr_segments = '$SSWDB/hinode/gen/hinode_sdo_alignment/sot_corr_segment_db_test' else $
  dir_corr_segments = '/sanhome/slater/public_html/outgoing/slater/sot_corr_testing'
; dir_corr_segments = '$SSWDB/hinode/gen/hinode_sdo_alignment/sot_corr_segment_db'
endif

if not exist(level0) then level0 = 1

if not exist(search_array) then $
  search_array = ['naxis1>511 && naxis2>511','wave=Ca*II*H*line']

if not exist(days_chunk) then days_chunk = 1
if not exist(delta_buff_minutes) then delta_buff_minutes = 1

; TODO: Define diff_jump as some factor * local (running) standard deviation:
if not exist(xdisp_diff_max) then xdisp_diff_max = 0.5 ; 1.0
if not exist(ydisp_diff_max) then ydisp_diff_max = 0.5 ; 1.0

if not exist(ctmex_diff_max) then ctmex_diff_max = 800.0 ; 1000.0
if not exist(ctmey_diff_max) then ctmey_diff_max = 800.0 ; 1000.0

if not exist(t_corr_diff_max) then t_corr_diff_max = 300 ; seconds

; Set default plotting parameters:
if not exist(yrange) then yrange = [-10,10]
if not exist(yrange2) then yrange2 = [-20000,20000]
if not exist(charsize) then charsize = 3

if not exist(loud) then loud = 1
if not exist(quiet_where) then quiet_where = 1

if keyword_set(time_driven) then begin
  if ( exist(t0) and exist(t1) ) then begin
    t_grid = anytim(timegrid(t0, t1, days=days_chunk), /ccsds)
    n_grid = n_elements(t_grid)-1
  endif else begin
    print, "One or both of 't0', 't1', not passed.  returning."
    return ; TODO - Set error here
  endelse
endif else begin
  n_grid = n_elements(files_corr_db)
  t_grid = strarr(n_grid+1)
endelse

for i=0, n_grid-1 do begin

  delvarx, t_seg_arr

  if keyword_set(time_driven) then begin

; Read all sdo correlation records in buffered grid interval, filter and skip if
; less then 2 records in the interval:
    t0_buff = anytim(anytim(t_grid[i])-delta_buff_minutes*60l, /ccsds)
    t1_buff = anytim(anytim(t_grid[i+1])+delta_buff_minutes*60l, /ccsds)
    read_genxcat, t0_buff, t1_buff, cat_corr, topdir=dir_top, count=count_corr

    if count_corr gt 1 then begin
      t_corr = cat_corr.date_obs

      ss_good_time = where( (anytim(t_corr) ge anytim(t_grid[i])) and $
                            (anytim(t_corr) le anytim(t_grid[i+1])), n_good_time)
      if n_good_time gt 1 then begin
        cat_corr = cat_corr[ss_good_time]

        ss_good = struct_where(cat_corr, search_array=search_array, count_good, $
                               quiet=quiet_where)
        if count_good gt 1 then begin
          cat_corr = cat_corr[ss_good]
          count_corr = count_good
        endif else begin
;         print, ' Interval: ' + t_grid[i] + ' to ' + t_grid[i+1] + ':'
          print, ' Less than 2 matches found for specified SEARCH_ARRAY. Skipping.'
          goto, SKIP_INTERVAL
        endelse
  
      endif else begin
;       print, ' Interval: ' + t_grid[i] + ' to ' + t_grid[i+1] + ':'
        print, ' Less than 2 matches found for this time interval. Skipping.'
        goto, SKIP_INTERVAL
      endelse

    endif else begin
;     print, ' Interval: ' + t_grid[i] + ' to ' + t_grid[i+1] + ':'
      print, ' Less than 2 matches found for this time interval. Skipping.'
      goto, SKIP_INTERVAL
    endelse

  endif else begin

    rd_genx, files_corr_db[i], buff
    cat_corr = buff.savegen0
    count_corr = n_elements(cat_corr)
    t_grid[i]   = anytim(cat_corr[0].date_obs, /ccsds)
    t_grid[i+1] = anytim(cat_corr[count_corr-1].date_obs, /ccsds)
    t0 = t_grid[i]
    t1 = t_grid[i+1]

  endelse

  t_fg_start = get_fg_start(t0, t1)
  if t_fg_start[0] ne 'NO_MATCHING_TIMES' then $
    n_fg_start = n_elements(t_fg_start) else n_fg_start = 0
STOP
; Sort and select uniq records:
  t_corr = cat_corr.date_obs
  tsec_corr = anytim(t_corr)
  ss_sort_corr = sort(tsec_corr)
  cat_corr = cat_corr[ss_sort_corr]

  t_corr = cat_corr.date_obs
  tsec_corr = anytim(t_corr)
  ss_uniq_corr = uniq(tsec_corr)
  cat_corr = cat_corr[ss_uniq_corr]

  t_corr = cat_corr.date_obs
  tsec_corr = anytim(t_corr)
  n_vals_corr = n_elements(cat_corr)
  if not exist(med_wid) then med_wid = n_vals_corr<9 ; <10

; Calculate jumps based upon absolute value of differences in adjacent correlation
; offsets in either x or y axes:
  xdisp = cat_corr.xdisp_rot
  xdisp_med = median(xdisp,med_wid)
  xdisp_diff = xdisp[1:*] - xdisp
  ss_jump_corr_x = where(abs(xdisp_diff) ge xdisp_diff_max, n_jump_corr_x)

  ydisp = cat_corr.ydisp_rot
  ydisp_med = median(ydisp,med_wid)
  ydisp_diff = ydisp[1:*] - ydisp
  ss_jump_corr_y = where(abs(ydisp_diff) ge ydisp_diff_max, n_jump_corr_y)

; Calculate jumps based upon absolute value of differences in adjacent correlation
; offsets in either x or y axes:
  ctmex = cat_corr.ctmex
  ctmex_med = median(ctmex,med_wid)
  ctmex_diff = ctmex[1:*] - ctmex
  ss_jump_ctmex = where(abs(ctmex_diff) ge ctmex_diff_max, n_jump_ctmex)

  ctmey = cat_corr.ctmey
  ctmey_med = median(ctmey,med_wid)
  ctmey_diff = ctmey[1:*] - ctmey
  ss_jump_ctmey = where(abs(ctmey_diff) ge ctmey_diff_max, n_jump_ctmey)

; Calculate jumps based upon time differences in adjacent correlation offsets:
  t_corr_diff_sec = tsec_corr[1:*] - tsec_corr
  ss_jump_t_corr_diff = where(abs(t_corr_diff_sec) ge t_corr_diff_max, n_jump_t_corr_diff)

; Calculate final jumps based up all correlation offset jumps, CT jumps, and
; time jumps:

  if keyword_set(use_corr) then begin
    if n_jump_corr_x gt 0 then ss_jump = ss_jump_corr_x
    if n_jump_corr_y gt 0 then begin
      if not exist(ss_jump) then ss_jump = ss_jump_corr_y else $
                                 ss_jump = [ss_jump, ss_jump_corr_y]
    endif
  endif

  if keyword_set(use_ct) then begin
    if n_jump_ctmex gt 0 then begin
      if not exist(ss_jump) then ss_jump = ss_jump_ctmex else $
                                 ss_jump = [ss_jump, ss_jump_ctmex]
    endif
    if n_jump_ctmey gt 0 then begin
      if not exist(ss_jump) then ss_jump = ss_jump_ctmey else $
                                 ss_jump = [ss_jump, ss_jump_ctmey]
    endif
  endif

  if keyword_set(use_time) then begin
    if n_jump_t_corr_diff gt 0 then begin
      if not exist(ss_jump) then ss_jump = ss_jump_t_corr_diff else $
                                 ss_jump = [ss_jump, ss_jump_t_corr_diff]
    endif
  endif

  ss_jump = UNIQ(ss_jump, SORT(ss_jump))
  n_jump = n_elements(ss_jump)
  
  if n_jump gt 0 then begin
;   ss_seg_st = [0,[ss_jump+1]]
    ss_seg_st = [0,[ss_jump  ]]
;   ss_seg_en = [[ss_jump  ],n_vals_corr-1]
    ss_seg_en = [[ss_jump-1],n_vals_corr-1]
    t0_seg = [t0,t_corr[ss_jump]]
    t1_seg = [t_corr[ss_jump],t1]

    for j=0, n_jump do begin

      cat_corr_seg = cat_corr[ss_seg_st[j]:ss_seg_en[j]]
      t_seg = t_corr[ss_seg_st[j]:ss_seg_en[j]]
      t_seg_sec = anytim(t_seg)

      n_pts_seg = n_elements(cat_corr_seg)
      med_wid_seg  = n_pts_seg<5 ; <3
      smoo_wid_seg = (n_pts_seg-1)<5 ; n_pts_seg<5

if n_pts_seg gt 1 then begin

      xdisp_seg = xdisp[ss_seg_st[j]:ss_seg_en[j]]
;     xdisp_med_seg = xdisp_med[ss_seg_st[j]:ss_seg_en[j]]
      xdisp_med_seg = median(xdisp_seg, med_wid_seg)
      xdisp_seg_smooth = smooth(xdisp_med_seg, smoo_wid_seg)

      ydisp_seg = ydisp[ss_seg_st[j]:ss_seg_en[j]]
;     ydisp_med_seg = ydisp_med[ss_seg_st[j]:ss_seg_en[j]]
      ydisp_med_seg = median(ydisp_seg, med_wid_seg)
      ydisp_seg_smooth = smooth(ydisp_med_seg, smoo_wid_seg)
      
      if keyword_set(do_plot) then begin
        wset, 0
        cleanplot, /silent
        !p.multi = [0,1,2]
        utplot,  t_seg, xdisp_seg, timerange=[t_corr[ss_seg_st[j]],t_corr[ss_seg_en[j]]], $
                 /xstyle, /ynoz, psym= 4, charsize=charsize, symsize=symsize
        outplot, t_seg, xdisp_seg_smooth, psym=-1, symsize=1, linestyle=2, thick=2
        outplot, t_seg, xdisp_med_seg, psym=-2, symsize=1, linestyle=0, thick=2

        utplot,  t_seg, ydisp_seg, timerange=[t_corr[ss_seg_st[j]],t_corr[ss_seg_en[j]]], $
                 /xstyle, /ynoz, psym= 4, charsize=charsize, symsize=symsize
        outplot, t_seg, ydisp_seg_smooth, psym=-1, symsize=symsize, linestyle=2, thick=2
        outplot, t_seg, ydisp_med_seg, psym=-2, symsize=symsize, linestyle=0, thick=2

        stop
      endif

      if not exist(t_seg_arr) then begin
        sttim_seg_arr = t0_seg[j]
        entim_seg_arr = t1_seg[j]
        t_seg_arr = t_seg
        xdisp_seg_arr = xdisp_seg
        xdisp_med_seg_arr = xdisp_med_seg
        xdisp_seg_smooth_arr = xdisp_seg_smooth
        ydisp_seg_arr = ydisp_seg
        ydisp_med_seg_arr = ydisp_med_seg
        ydisp_seg_smooth_arr = ydisp_seg_smooth
        cat_corr_seg_arr = cat_corr_seg
      endif else begin
        sttim_seg_arr = [sttim_seg_arr, t0_seg[j]]
        entim_seg_arr = [entim_seg_arr, t1_seg[j]]
        t_seg_arr = [t_seg_arr, t_seg]
        xdisp_seg_arr = [xdisp_seg_arr, xdisp_seg]
        xdisp_med_seg_arr = [xdisp_med_seg_arr, xdisp_med_seg]
        xdisp_seg_smooth_arr = [xdisp_seg_smooth_arr, xdisp_seg_smooth]
        ydisp_seg_arr = [ydisp_seg_arr, ydisp_seg]
        ydisp_med_seg_arr = [ydisp_med_seg_arr, ydisp_med_seg]
        ydisp_seg_smooth_arr = [ydisp_seg_smooth_arr, ydisp_seg_smooth]
        cat_corr_seg_arr = concat_struct(cat_corr_seg_arr, cat_corr_seg)
      endelse

endif else begin
      print, 'Only a single point in this segment. Skipping.'
      print, strtrim(t_seg[0],2) + ' to ' + strtrim(t_seg[n_pts_seg-1],2) + '.'
endelse

    endfor

    if keyword_set(do_all_seg_plot) then begin

      wset, 1
      cleanplot, /silent
      !p.multi = [0,1,2]
      trange = [t_grid[i], t_grid[i+1]]
      utplot, t_corr, xdisp, timerange=trange, yrange=yrange, $
              /xstyle, /ystyle, psym=1
      outplot, t_corr, xdisp_med

; Draw boundary lines for each segment:
      if n_jump gt 0 then begin
        for j=0, n_jump do begin
          outplot, [t_corr[ss_seg_en[j]], t_corr[ss_seg_en[j]]], yrange, psym=0, thick=2
        endfor
      endif

; Plot times of FG_STARTs;
if n_fg_start gt 0 then $
  for j=0, n_fg_start-1 do outplot, [t_fg_start[j], t_fg_start[j]], yrange, psym=0, thick=2

      utplot, t_corr, ydisp, timerange=trange, yrange=yrange, $
              /xstyle, /ystyle, psym=1
      outplot, t_corr, ydisp_med

; Draw boundary lines for each segment:
      if n_jump gt 0 then begin
        for j=0, n_jump do begin
          outplot, [t_corr[ss_seg_en[j]], t_corr[ss_seg_en[j]]], yrange, psym=0, thick=2
        endfor
      endif

; Plot times of FG_STARTs;
if n_fg_start gt 0 then $
  for j=0, n_fg_start-1 do outplot, [t_fg_start[j], t_fg_start[j]], yrange, psym=0, thick=2

    endif

    if keyword_set(do_all_seg_plot) then begin
 
      wdef,2,1024+512,512+512+256
      cleanplot, /silent
      !p.multi = [0,1,5]

xdisp_med_mom = moment(xdisp_med, mean=xdisp_med_mean, sdev=xdisp_med_sdev)
yrange_xdisp = 4*[-xdisp_med_sdev, xdisp_med_sdev] + xdisp_med_mean

      trange = [t_grid[i], t_grid[i+1]]
      utplot, t_corr, xdisp, timerange=trange, yrange=yrange_xdisp, $
              /xstyle, /ystyle, charsize=charsize, ytit='x_corr (arcsec)', psym=1
      outplot, t_corr, xdisp_med
      if n_jump gt 0 then begin
        for j=0, n_jump do begin
          outplot, [t_corr[ss_seg_en[j]], t_corr[ss_seg_en[j]]], $
                   yrange, psym=0, thick=1, linestyle=2
        endfor
      endif

      utplot, t_corr, ctmex, timerange=trange, $ ;yrange=yrange2, $
              /xstyle, /ystyle, charsize=charsize, ytit='ctmex', psym=1
      outplot, t_corr, ctmex_med
      if n_jump gt 0 then begin
        for j=0, n_jump do begin
          outplot, [t_corr[ss_seg_en[j]], t_corr[ss_seg_en[j]]], $
                   yrange2, psym=0, thick=1, linestyle=2
        endfor
      endif

; Plot times of FG_STARTs;
if n_fg_start gt 0 then $
  for j=0, n_fg_start-1 do outplot, [t_fg_start[j], t_fg_start[j]], yrange, psym=0, thick=2

ydisp_med_mom = moment(ydisp_med, mean=ydisp_med_mean, sdev=ydisp_med_sdev)
yrange_ydisp = 4*[-ydisp_med_sdev, ydisp_med_sdev] + ydisp_med_mean

      utplot, t_corr, ydisp, timerange=trange, yrange=yrange_ydisp, $
              /xstyle, /ystyle, charsize=charsize, ytit='y_corr (arcsec)', psym=1
      outplot, t_corr, ydisp_med
      if n_jump gt 0 then begin
        for j=0, n_jump do begin
          outplot, [t_corr[ss_seg_en[j]], t_corr[ss_seg_en[j]]], $
                   yrange, psym=0, thick=1, linestyle=2
        endfor
      endif

; Plot times of FG_STARTs;
if n_fg_start gt 0 then $
  for j=0, n_fg_start-1 do outplot, [t_fg_start[j], t_fg_start[j]], yrange, psym=0, thick=2

      utplot, t_corr, ctmey, timerange=trange, $ ;yrange=yrange2, $
              /xstyle, /ystyle, charsize=charsize, ytit='ctmey', psym=1
      outplot, t_corr, ctmey_med
      if n_jump gt 0 then begin
        for j=0, n_jump do begin
          outplot, [t_corr[ss_seg_en[j]], t_corr[ss_seg_en[j]]], $
                   yrange2, psym=0, thick=1, linestyle=2
        endfor
      endif

; Plot times of FG_STARTs;
if n_fg_start gt 0 then $
  for j=0, n_fg_start-1 do outplot, [t_fg_start[j], t_fg_start[j]], yrange, psym=0, thick=2

      ct_xy_rms = sqrt(ctmex*ctmex + ctmey*ctmey)
      ct_xy_rms_med = median(ct_xy_rms, med_wid)
      utplot, t_corr, ct_xy_rms, timerange=trange, $ ;yrange=yrange2, $
              /xstyle, /ystyle, charsize=charsize, ytit='ct_xy_rms', psym=1
      outplot, t_corr, ct_xy_rms_med
      if n_jump gt 0 then begin
        for j=0, n_jump do begin
          outplot, [t_corr[ss_seg_en[j]], t_corr[ss_seg_en[j]]], $
                   yrange2, psym=0, thick=1, linestyle=2
        endfor
      endif

; Plot times of FG_STARTs;
for j=0, n_fg_start-1 do outplot, [t_fg_start[j], t_fg_start[j]], yrange, psym=0, thick=2

    endif

; Optionally write out save file for fitted segment data:
    if keyword_set(do_save) then begin
;     t0_filnam = time2file(t_corr[ss_seg_st[0]])
;     t1_filnam = time2file(t_corr[ss_seg_en[n_jump]])
;     filnam_sav = 'sot_corr_segments_' + t0_filnam + '_' + t1_filnam
      t0_filnam = time2file(t_grid[i], /date)
      filnam_sav = 'sot_corr_segments_' + t0_filnam
      filnam_sav = concat_dir(dir_corr_segments, filnam_sav)
      save, sttim_seg_arr, entim_seg_arr, t_seg_arr, xdisp_seg_arr, xdisp_med_seg_arr, $
            xdisp_seg_smooth_arr, ydisp_seg_arr, ydisp_med_seg_arr, ydisp_seg_smooth_arr, $
            cat_corr_seg_arr, filename=filnam_sav
    endif

  endif else begin
;   print, ' Interval: ' + t_grid[i] + ' to ' + t_grid[i+1] + ':'
    print, ' No jumps found in this time interval. Skipping.'
  endelse

  SKIP_INTERVAL:

  if keyword_set(loud) then $
    print, ' Finished for ' + anytim(t_grid[i], /yoh) + ' to ' + anytim(t_grid[i+1], /yoh) 

  if keyword_set(qstop) then stop

endfor

if keyword_set(qstop) then stop

end
