pro sdo_corr, index_sot, data_sot, indexp_sot, datap_sot, file_sot=file_sot, $
  index_sdo, data_sdo, indexp_sdo, datap_sdo, file_sdo=file_sdo, $
  ss_axis3=ss_axis3, ss_axis4=ss_axis4, $
  do_restore_test=do_restore_test, $ 
  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_plot_map=do_plot_map, do_xstep=do_xstep, do_blink=do_blink, $
  verbose=verbose, _extra=_extra

;+
; NAME:
;	sdo_corr
; PURPOSE:
;	Cross correlate input image with fov-and-time-matched SDO cutout image
;       of appropriate wavelength, to determine translational offsets.
; CATEGORY:
;	Image registration
; CALL:
;	sdo_corr, index_sot, data_sot
; OPTIONAL KEYWORDS INPUTS:
;	verbose
;	_extra
; OUTPUTS:
;	None
; OPTIONAL KEYWORDS OUTPUTS:
;	sot_align_struct - structure containing various parameters related to the
;                          IRIS and AIA images and the cross correlation results:
;	  date_obs_sot
;	  date_obs_sdo
;	  file_sdo
;	  x_off
;	  y_off
;	  wave_sot
;	  wave_sdo
;	  raster_xsize
;	  raster_ysize
;	  peak_str
; HISTORY:
;       2013-11-21 - JPW : - added logic to correct for AIA CROTA2
;                          - removed loop around interpolate
;                          - changed logic to rotate xdisp with PC matrix
;                            and to properly add it to the old CRVALi
;       2014-01-13 - JPW : - now uses an on-disk subfield for images at the limb
;                            turn off with /nocrop2disk
;                          - also fixed xcen,ycen and removed scale stuff and
;                            some other unused code
; TODO:
;       - Change instrument-specific parameter names into generic ones:
;       - Handle vectors of images
;       - Caching of reusables (fit coefficients, etc)
;       - Auto-scale SDO images using same method as for SOT
;-

if keyword_set(do_restore_test) then restore, './sdo_corr_20100628.sav'

if not exist(max_t_diff_sec) then max_t_diff_sec = 60
if not exist(n_smooth) then n_smooth = 5
if not exist(peak_str_min) then peak_str_min = 3
if not exist(wait_glitch) then wait_glitch = 0.4
if not exist(wait_sec) then wait_sec = 0 ; 0.4
if not exist(n_blink) then n_blink = 10
if not exist(wait_blink) then wait_blink = 0 ; 0.1

t_sot = index_sot.date_obs
sat_rot = index_sot.sat_rot
sunpb0r = pb0r(index_sot.date_obs)
rsun = sunpb0r[2]*60.

; Determine which image to use in s 3+ dimensional image cube:
; (Use temporary kluge table to determine when to use first image in cube)
; TODO: Make this logic less klugey.) 

; Totally Ad Hoc Kluge 1:
mode_sot = index_sot.wave + ' ' + index_sot.obs_type
if mode_sot eq 'TF Na I 5896 FG shutterless I and V with 0.2s intervals 1' then ss_axis3 = 0

if tag_exist(index_sot, 'naxis3') then begin
  if not exist(ss_axis3) then ss_axis3 = index_sot.naxis3-1
endif
if tag_exist(index_sot, 'naxis4') then begin
  if not exist(ss_axis4) then ss_axis4 = index_sot.naxis4-1
endif

; Reduce the cube to a single 2D image:
if ( tag_exist(index_sot, 'naxis3') and (not tag_exist(index_sot, 'naxis4')) ) then $
  data_sot = reform(data_sot[*,*,ss_axis3])
if tag_exist(index_sot, 'naxis4') then $
  data_sot = reform(data_sot[*,*,ss_axis3,ss_axis4])

if mode_sot eq 'TF Na I 5896 FG shutterless I and V with 0.2s intervals 1' then begin
  ss_neg = where(data_sot lt 0, n_neg)
  data_sot[ss_neg] = -data_sot[ss_neg]
endif

; Optionally prep SOT image:
;if keyword_set(do_prep_sot) then begin
;  fg_prep, index_sot, data_sot, indexp_sot, datap_sot
;endif else begin
  indexp_sot = index_sot
  datap_sot = data_sot
;  ; In case of no fg_prep, then update pointing keywords.
;  ; For now, use Shimizu table values, eventually to be replaced by
;  ;   continuously and iteratively improved table from the correlation results:
;  sot_offsets = get_shimizu(t_sot)
;PRINT, SOT_OFFSETS
;  crval1_offset = sot_offsets[0]
;  crval2_offset = sot_offsets[1]
;  indexp_sot.crval1 = indexp_sot.crval1 + crval1_offset
;  indexp_sot.crval2 = indexp_sot.crval2 + crval2_offset
;  indexp_sot.xcen = indexp_sot.crval1
;  indexp_sot.ycen = indexp_sot.crval2
;endelse

; Add CROTA2 tag to IRIS header if not present:
;if tag_exist(indexp_sot, 'crota2') then indexp_sot.crota2 = sat_rot 

; Update rotation angle from fit to measured data:
if not keyword_set(norollcorect) then begin
  restore, '/net/topaz/Users/slater/soft/idl/idl_startup/sot/coeff_fit_sot_delta_rotang.sav'
; indexp_sot.crota1 = indexp_sot.crota1 + (10*sine_func(anytim(t_sot), coeff_fit_da))
; indexp_sot.crota1 = indexp_sot.crota1 + (05*sine_func(anytim(t_sot), coeff_fit_da))
  indexp_sot.crota1 = indexp_sot.crota1 - 0.8
  indexp_sot.crota2 = indexp_sot.crota1
; indexp_sot.crota1 = indexp_sot.crota1 + 5
; indexp_sot.crota2 = indexp_sot.crota2 + 5
endif

; Update platescale from fit to measured data:
restore, '/net/topaz/Users/slater/soft/idl/idl_startup/sot/coeff_fit_sot_delta_platescale.sav'
;indexp_sot.cdelt1 = indexp_sot.cdelt1 / (1.005*poly_func(anytim(t_sot), reform(coeff_fit_ds)))
;indexp_sot.cdelt2 = indexp_sot.cdelt1
;indexp_sot.cdelt1 = indexp_sot.cdelt1 / 1.02
 indexp_sot.cdelt1 = indexp_sot.cdelt1 / 1.015
 indexp_sot.cdelt2 = indexp_sot.cdelt1

; Generate regulation WCS structure from SSW index structure for SOT image:
wcsp_sot = fitshead2wcs(indexp_sot)

; Define size and FOV of IRIS image (CCD readout region):
tsc1 = 1
tec1 = indexp_sot.naxis1
tsr1 = 1
ter1 = indexp_sot.naxis2
iNx = tec1-tsc1+1
iNy = ter1-tsr1+1

; Define arrays to contain x and y Helioprojective (arcsec) coordinates of all image pixels:
solx = fltarr(iNx,iNy)
soly = fltarr(iNx,iNy)

for i=0,iNx-1 do begin
	for j=0,iNy-1 do begin
		solx[i,j] = (tsc1 + i - indexp_sot.crpix1)*wcsp_sot.pc[0,0] + $
					(tsr1 + j - indexp_sot.crpix2)*wcsp_sot.pc[0,1]
		soly[i,j] = (tsc1 + i - indexp_sot.crpix1)*wcsp_sot.pc[1,0] + $
					(tsr1 + j - indexp_sot.crpix2)*wcsp_sot.pc[1,1]
	endfor
endfor

solx *= indexp_sot.cdelt1
soly *= indexp_sot.cdelt2
solx += indexp_sot.crval1
soly += indexp_sot.crval2

; Optionally check if limb crosses IRIS cutout field.  If so, crop to on-disk portion:
crop_flag = 0
if ~keyword_set(nocrop2disk) then begin
  fovmar = 5.0
  fovmin = 20.0
  ; 5 arcsec margin:
  mx_solr = rsun - fovmar
  solr = sqrt(solx*solx+soly*soly)
  ; corners of IRIS image :
  solr_corn = solr[[[0,iNx-1],[0,iNx-1]],[[0,0],[iNy-1,iNy-1]]]
  if min(solr_corn) lt mx_solr and max(solr_corn) gt mx_solr then begin
    ; corner closest to suncenter "cctsc" :
    wcorn = where(solr_corn eq min(solr_corn))
    cornx = wcorn[0] mod 2
    corny = wcorn[0] / 2
    ; find points where solar limb (minus margin) crosses edges of image
    ; edges adjacent to cctsc (only need nxe1,nye1 = # of edge pixels inside limb)
    dmy = where(solr[*,corny*(iNy-1)] lt mx_solr,nxe1)
    dmy = where(solr[cornx*(iNx-1),*] lt mx_solr,nye1)
    ; edges not adjacent to cctsc
    dmy = where(solr[*,(1-corny)*(iNy-1)] lt mx_solr,nxe2)
    dmy = where(solr[(1-cornx)*(iNx-1),*] lt mx_solr,nye2)
    ; Use limb crossing pts to determine optimum subfield size.
    ; Assuming straight limb for this calculation, which provides
    ; a somewhat conservative result (draw 4 cases to derive eq's)
    case 1 of                                  ; 4 cases
      (nxe1 eq inx) and (nye1 eq iny): begin   ; crossing only non-adjacent edges
        inx1 = nxe2 > fix((float(inx)+float(nye2)*float(inx-nxe2)/float(iny-nye2))/2.0) < inx
        iny1 = nye2 > fix((float(iny)+float(nxe2)*float(iny-nye2)/float(inx-nxe2))/2.0) < iny
        end
      (nxe1 eq inx) and (nye1 lt iny): begin   ; crossing both y edges
        inx1 = fix(float(nye1)*float(inx)/float(2*(nye1-nye2) > nye1))
        iny1 = nye1/2 > nye2
        end
      (nxe1 lt inx) and (nye1 eq iny): begin   ; crossing both x edges
        inx1 = nxe1/2 > nxe2
        iny1 = fix(float(nxe1)*float(iny)/float(2*(nxe1-nxe2) > nxe1))
        end
      else: begin                              ; crossing only adjacent edges
        inx1 = nxe1/2
        iny1 = nye1/2
        end
    endcase
    ; update IRIS image cutout if new subimage large enough (i.e., at least 20 " on the side)
    iFOVx1 = inx1*indexp_sot.cdelt1
    iFOVy1 = iny1*indexp_sot.cdelt2
    if iFOVx1 ge fovmin and iFOVy1 ge fovmin then begin
      otsc1 = cornx * (iNx-inx1)
      otsr1 = corny * (iNy-iny1)
      tsc1 = tsc1+otsc1
      tec1 = tsc1+inx1-1
      tsr1 = tsr1+otsr1
      ter1 = tsr1+iny1-1
      iNx = inx1
      iNy = iny1
      solx = solx[otsc1:otsc1+inx1-1,otsr1:otsr1+iny1-1]
      soly = soly[otsc1:otsc1+inx1-1,otsr1:otsr1+iny1-1]
      crop_flag = 1
    endif
  endif
endif

; Define SOT cutout:
;cutoutp_sot = datap_sot[tsc1-1:tec1-1,tsr1-1:ter1-1] / indexp_sot.exptime
 cutoutp_sot = datap_sot[tsc1-1:tec1-1,tsr1-1:ter1-1]

; Handle NaN's:
ss_fin = where(finite(cutoutp_sot) eq 1, n_fin, comp=ss_inf, ncomp=n_inf)
if n_inf gt 0 then cutoutp_sot[ss_inf] = min(cutout_sotp[ss_fin])

; Find SDO wave to match, as well as scaling information for both SOT and SDO images:
wave_aia = aia_wave_match(indexp_sot, mode=mode, 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)

; Return if no match found:
if wave_aia eq '' then return

; 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.'
  cutoutp_sot_scaled = bytscl(cutoutp_sot)
endif else begin
  if (do_invert_sot or keyword_set(force_invert_sot)) then $
    cutoutp_sot_scaled = 255b-bytscl(loval_sot>cutoutp_sot<hival_sot) else $
    cutoutp_sot_scaled =      bytscl(loval_sot>cutoutp_sot<hival_sot)
endelse

if keyword_set(do_img_plot) then begin
  wdef,0,indexp_sot.naxis1,indexp_sot.naxis2
  tvscl, cutoutp_sot_scaled
endif

hist = histogram(cutoutp_sot, loc=xloc)
ss_nz = where(hist ne 0, n_nz)

if n_nz gt 0 then begin
  hist_nz = hist[ss_nz]
  xloc_nz = xloc[ss_nz]
; For now, at least, set hist, xloc to nz versions:
  hist = hist_nz
  xloc = xloc_nz
  if n_elements(hist) gt n_smooth then $
    hist_integral_smooth = smooth(total(hist, /cum), n_smooth) / total(hist) else $
    return
;   hist_integral_smooth = total(hist, /cum) / total(hist)
; deriv_hist_integral_smooth =  hist_integral_smooth[1:*] - hist_integral_smooth
  deriv_hist_integral_smooth = [hist_integral_smooth[1:*] - hist_integral_smooth,0]
; Try to fit the dominant guassian to automatically determine lo and hi cutoff values:
  n_pts = n_elements(hist)
  n0_crop = fix(0.1*n_pts)
  x_crop = xloc[n0_crop:*]
; hist_crop = hist[n0_crop:*]
  hist_crop = deriv_hist_integral_smooth[n0_crop:*]
  y_crop = hist_crop
  yfit_crop = gaussfit(x_crop, y_crop, afit, chisq=chisqr, nterms=3, sigma=sigma, yerror=yerror)
  loval_sot_auto = afit[1] - nsig_sot*afit[2]
  hival_sot_auto = afit[1] + nsig_sot*afit[2]

  if (do_invert_sot or keyword_set(force_invert_sot)) then $
    cutoutp_sot_scaled_auto = 255b-bytscl(loval_sot_auto>cutoutp_sot<hival_sot_auto) else $
    cutoutp_sot_scaled_auto =      bytscl(loval_sot_auto>cutoutp_sot<hival_sot_auto)

  if keyword_set(do_hist_plot) then begin
    wdef,1,1024,768
    !p.multi = [0,1,3]
    plot, xloc, hist, charsize=2
    plot, xloc, hist_integral_smooth, charsize=2
    plot, xloc, deriv_hist_integral_smooth, charsize=2

    oplot, x_crop, yfit_crop, psym=2

    oplot, [loval_sot_auto, loval_sot_auto], [0, max(yfit_crop)], linestyle=4, thick=2
    oplot, [hival_sot_auto, hival_sot_auto], [0, max(yfit_crop)], linestyle=4, thick=2
  endif

endif else begin
  stop, ' No non-zero values of histogram. Stopping.'
endelse

if keyword_set(do_img_plot) then begin
  wset,0
  tvscl, cutoutp_sot_scaled_auto
endif

; Construct call to ssw_jsoc_time2data to get SDO image:
if not exist(index_sdo) then $
  ssw_jsoc_time2data, anytim(t_sot, /ccsds), anytim(anytim(t_sot)+100, /ccsds), $
    wave=wave_aia, index_sdo, file_sdo, /files_only
; ssw_jsoc_time2data, anytim(t_sot, /ccsds), anytim(anytim(t_sot)+300, /ccsds), $
;   wave=wave_aia, index_sdo, file_sdo, /files_only
; ssw_jsoc_time2data, t_sot, t_sot, wave=wave_aia, index_sdo, file_sdo, /files_only

  ; If no file found, try one more time, in case there was a system glitch:
if ( (file_sdo[0] eq '') or (file_sdo[0] eq 'Remote listing failed.') ) then begin
  wait, wait_glitch
  ssw_jsoc_time2data, anytim(t_sot, /ccsds), anytim(anytim(t_sot)+100, /ccsds), $
    wave=wave_aia, index_sdo, file_sdo, /files_only
; ssw_jsoc_time2data, anytim(t_sot, /ccsds), anytim(anytim(t_sot)+300, /ccsds), $
;   wave=wave_aia, index_sdo, file_sdo, /files_only
; ssw_jsoc_time2data, t_sot, t_sot, wave=wave_aia, index_sdo, file_sdo, /files_only
endif

if all_vals(file_exist(file_sdo)) eq [-1] then return

read_sdo, file_sdo, index_sdo, data_sdo, /use_index, /uncomp_delete

; If read_sdo failed, try once more:
if not exist(index_sdo) then begin
  print, 'read_sdo failed.  Waiting 100 seconds and trying again.'
  wait, wait_glitch
  read_sdo, file_sdo, index_sdo, data_sdo, /use_index, /uncomp_delete
endif

if ( (not exist(index_sdo)) or (not exist(data_sdo)) ) then return

; ssw_jsoc_time2data, t_sot, t_sot, wave=wave_aia, index_sdo, data_sdo
; ssw_jsoc_time2data, anytim(t_sot, /ccsds), anytim(anytim(t_sot)+300, /ccsds), $
;   wave=wave_aia, index_sdo, data_sdo

if not tag_exist(index_sdo, 'crota2') then return

file_sdo = file_sdo[0]
index_sdo = index_sdo[0]
data_sdo  = reform(data_sdo[*,*,0])

; Optionally prep SDO image:
if keyword_set(do_prep_sdo) then begin
  aia_prep, index_sdo, data_sdo, indexp_sdo, datap_sdo
endif else begin
  indexp_sdo = index_sdo
  datap_sdo = data_sdo
endelse

; Rotate SOT cutout coordinates (solx/soly) to the SDO roll:
solx_sdo =  solx*cos(indexp_sdo.crota2/!radeg) + soly*sin(indexp_sdo.crota2/!radeg)
soly_sdo = -solx*sin(indexp_sdo.crota2/!radeg) + soly*cos(indexp_sdo.crota2/!radeg)

; Translate SOT cutout coordinates to SDO pixel coords
solx_sdo = solx_sdo/indexp_sdo.cdelt1 + indexp_sdo.crpix1 - 1 
soly_sdo = soly_sdo/indexp_sdo.cdelt2 + indexp_sdo.crpix2 - 1

; Define matching SDO cutout using cubic interpolation on full SDO image:
cutoutp_sdo = interpolate(datap_sdo, solx_sdo, soly_sdo, cubic=-0.5)

; Determine clipping and scaling for SDO image:
if not exist(nsig_sdo) then nsig_sdo = nsig_sdo_tab

if ( (loval_sdo eq 0) and (hival_sdo eq 0) ) then begin
  print, 'No limit values for SDO image in database.'
  cutoutp_sdo_scaled = bytscl(cutoutp_sdo)
endif else begin
  if (do_invert_sdo or keyword_set(force_invert_sdo)) then $
    cutoutp_sdo_scaled = 255b-bytscl(loval_sdo>cutoutp_sdo<hival_sdo) else $
    cutoutp_sdo_scaled =      bytscl(loval_sdo>cutoutp_sdo<hival_sdo)
endelse

if keyword_set(do_img_plot) then begin
  wset,0
  tvscl, cutoutp_sdo_scaled
endif

hist = histogram(cutoutp_sdo, loc=xloc)
ss_nz = where(hist ne 0, n_nz)

if n_nz gt 0 then begin
  hist_nz = hist[ss_nz]
  xloc_nz = xloc[ss_nz]
; For now, at least, set hist, xloc to nz versions:
  hist = hist_nz
  xloc = xloc_nz
  if n_elements(hist) gt n_smooth then $
    hist_integral_smooth = smooth(total(hist, /cum), n_smooth) / total(hist) else $
    return
;   hist_integral_smooth = total(hist, /cum) / total(hist)
; deriv_hist_integral_smooth =  hist_integral_smooth[1:*] - hist_integral_smooth
  deriv_hist_integral_smooth = [hist_integral_smooth[1:*] - hist_integral_smooth,0]

; Try to fit the dominant guassian to automatically determine lo and hi cutoff values:
  n_pts = n_elements(hist)
  n0_crop = fix(0.00*n_pts)
  x_crop = xloc[n0_crop:*]
; hist_crop = hist[n0_crop:*]
  hist_crop = deriv_hist_integral_smooth[n0_crop:*]
  y_crop = hist_crop
  yfit_crop = gaussfit(x_crop, y_crop, afit, chisq=chisqr, nterms=3, sigma=sigma, yerror=yerror)

  loval_sdo_auto = afit[1] - nsig_sdo*afit[2]
  hival_sdo_auto = afit[1] + nsig_sdo*afit[2]

  if (do_invert_sdo or keyword_set(force_invert_sdo)) then $
    cutoutp_sdo_scaled_auto = 255b-bytscl(loval_sdo_auto>cutoutp_sdo<hival_sdo_auto) else $
    cutoutp_sdo_scaled_auto =      bytscl(loval_sdo_auto>cutoutp_sdo<hival_sdo_auto)

  if keyword_set(do_hist_plot) then begin
    wdef,1,1024,768
    !p.multi = [0,1,3]
    plot, xloc, hist, charsize=2
    plot, xloc, hist_integral_smooth, charsize=2
    plot, xloc, deriv_hist_integral_smooth, charsize=2

    oplot, x_crop, yfit_crop, psym=2

    oplot, [loval_sdo_auto, loval_sdo_auto], [0, max(yfit_crop)], linestyle=4, thick=2
    oplot, [hival_sdo_auto, hival_sdo_auto], [0, max(yfit_crop)], linestyle=4, thick=2
  endif

endif else begin
  stop, ' No non-zero values of histogram. Stopping.'
endelse

if keyword_set(do_img_plot) then begin
  wset,0
  tvscl, cutoutp_sdo_scaled_auto

  ; Optionally display cutouts prior to correlation:
  if crop_flag eq 1 then begin
    if iNx gt 1024 then begin
      wdef,2,2*(iNx/2),iNy/2
      tvscl, congrid(cutoutp_sot_scaled_auto,iNx/2,iNy/2)
      tvscl, congrid(cutoutp_sdo_scaled_auto,iNx/2,iNy/2), iNx/2, 0
    endif else begin
      wdef,2,2*iNx,iNy
      tvscl, cutoutp_sot_scaled_auto
      wait, wait_sec
      tvscl, cutoutp_sdo_scaled_auto, iNx, 0
    endelse
  endif
endif

; Create SDO/IRIS cutout cube and cross correlate to get offsets:
;cutout_cube = [[[cutoutp_sdo_scaled_auto]],[[cutoutp_sot_scaled_auto]]]
 cutout_cube = [[[cutoutp_sdo_scaled]],[[cutoutp_sot_scaled]]]

disp = tr_get_disp(cutout_cube, /shift, do_surf=do_surf, peak_str=peak_str, _extra=_extra)

xdisp = disp[0,1]*indexp_sot.cdelt1
ydisp = disp[1,1]*indexp_sot.cdelt2

; Rotate the X and Y offsets from the correlation to the roll angle of the SOT image:
xdisp_rot = xdisp*wcsp_sot.pc[0,0] + ydisp*wcsp_sot.pc[0,1]
ydisp_rot = xdisp*wcsp_sot.pc[1,0] + ydisp*wcsp_sot.pc[1,1]

xcen_sot_cor = indexp_sot.xcen + xdisp_rot
ycen_sot_cor = indexp_sot.ycen + ydisp_rot

crval1_sot_cor = $
  comp_fits_crval(xcen_sot_cor, indexp_sot.cdelt1, indexp_sot.naxis1, indexp_sot.crpix1)
crval2_sot_cor = $
  comp_fits_crval(ycen_sot_cor, indexp_sot.cdelt2, indexp_sot.naxis2, indexp_sot.crpix2)

crpix1_sot_cor = $
  comp_fits_crpix(xcen_sot_cor, indexp_sot.cdelt1, indexp_sot.naxis1, crval1_sot_cor)
crpix2_sot_cor = $
  comp_fits_crpix(ycen_sot_cor, indexp_sot.cdelt2, indexp_sot.naxis2, crval2_sot_cor)

; Update SOT X,YCEN, CRVAL1,2, CRPIX1,2 tags (add X,YCEN tags if they don't exist):
if not tag_exist(indexp_sot, 'xcen') then $
  indexp_sot = add_tag(indexp_sot, xcen_sot_cor, 'xcen') else $
  indexp_sot.xcen = xcen_sot_cor
if not tag_exist(indexp_sot, 'ycen') then $
  indexp_sot = add_tag(indexp_sot, ycen_sot_cor, 'ycen') else $
  indexp_sot.ycen = ycen_sot_cor

indexp_sot.crval1 = crval1_sot_cor
indexp_sot.crval2 = crval2_sot_cor
indexp_sot.crpix1 = crpix1_sot_cor
indexp_sot.crpix2 = crpix2_sot_cor

; Update SOT PC matrix tags:
; ( CODE HERE )

; Create minimmal sot structure and add correlation offsets and peak strength tags:
tag_list = 'naxis1, naxis2, date_obs, ctime, data_lev, orig_rf0, ver_rf0, ' + $
           'saa, hlz, waveid, obs_type, xscale, yscale, fgxoff, fgyoff, ' + $
           'fgccdix0, fgccdix1, fgccdiy0, fgccdiy1, crpix1, crpix2, ' + $
           'sc_attx, sc_atty, crval1, crval2, cdelt1, cdelt2, sat_rot, ' + $
           'inst_rot, crota1, crota2, xcen, ycen, fovx, fovy, fgbinx, fgbiny, ' + $
           'exptime, wave, darkflag, fgmode, fgnint, ctserv, ctmex, ctmey, ' + $
           'ctmode, t_spccd, t_fgccd, t_ctccd, focus, percentd'

align_str = str_subset(indexp_sot, tag_list)
align_str = add_tag(align_str, xdisp_rot, 'xdisp_rot')
align_str = add_tag(align_str, ydisp_rot, 'ydisp_rot')
align_str = add_tag(align_str, peak_str, 'peak_str')
align_str = add_tag(align_str, file_sot, 'file_sot')
align_str = add_tag(align_str, file_sdo, 'file_sdo')
align_str = add_tag(align_str, indexp_sot.date_obs, 't_sot')
align_str = add_tag(align_str, indexp_sdo.date_obs, 't_sdo')
dt_sdo = anytim(indexp_sdo.date_obs) - anytim(indexp_sot.date_obs)
align_str = add_tag(align_str, dt_sdo, 'dt_sdo')

if keyword_set(do_plot_map) then begin
  index2map, indexp_sot, datap_sot, mapp_sot
  plot_map, mapp_sot, center=[indexp_sot.crval1, indexp_sot.crval2], $
    fov=[02,02], grid=5, /limb, charsize=2
  img_sot = tvrd()
  index2map, indexp_sdo, datap_sdo, mapp_sdo
  plot_map, mapp_sdo, center=[indexp_sot.crval1, indexp_sot.crval2], $
    fov=[02,02], grid=5, /limb, charsize=2
  img_sdo = tvrd()
  img_cube = [[[img_sot]],[[img_sdo]]]
  xstepper, img_cube
endif

if ( keyword_set(do_xstep) and keyword_set(verbose) ) then xstepper, cutout_cube
if ( keyword_set(do_blink) and keyword_set(verbose) ) then begin
  siz_cutout = size(cutout_cube)
  wdef, 2, siz_cutout[1], siz_cutout[2]
  for i=0,n_blink-1 do begin
    tvscl, cutout_cube[*,*,0]
    wait, wait_blink
    tvscl, cutout_cube[*,*,1]
    wait, wait_blink
  endfor
endif

end