; file: nisp_sdo_tilefigs.idl = test differentially for heightdiff
; init: Mar 30 2022  Rob Rutten  Deil from gong_solis_tilefigs.idl
; last: Apr 11 2022  Rob Rutten  Deil
; note: only the first gives good heightdiff curve

;; pro dummy

; select data    ## must have nisp_sdp.pro temp files AIA + GONG/SOLIS
cd,'/home/rutten/data/NISP/2012-08-15-quiet'
datetime='2012.08.15_23:00'

; select which pair
combo=1 ; SOLIS to GONG
;; combo=2 ; GONG to AIA 304
;; combo=3 ; SOLIS to AIA 304

; set box size in arcsec for trimboxim2 coalignment around disk center
;RR if not small then for SOLIS sumwing upset by heightdiff > 0
;RR so that non-small-box align gives wrong scale and tiling goes bad
if (combo eq 1) then begin
  boxsize=100   ; 100 OK for align quiet center
  pxtile=50     ; default (30 SDO arcsec, case here)
  tilecut=2     ; default
  minleft=0.2   ; default
  ; beautiful, result heightdiff 2700 km
endif
if (combo eq 2) then begin
  boxsize=300   ; 300 OK for align quiet center but not enough valid tiles
  pxtile=100
  tilecut=2
  minleft=0.1
  ; no good
endif
if (combo eq 3) then begin
  boxsize=200   ; 200 OK for align quiet center but not enough valid tiles
  pxtile=100    ; default 50
  tilecut=3     ; default 2
  ; weird large scale vortex
endif 

; choose heightdiff
heightdiff=0.     ;RR solis>gong wow!  Indeed needs small trimbox in coalign.
;; heightdiff=2700.

; output choices
blink=0
checkmetcalf=0
show=1
verbose=1


; ============ rest should be parameter free ===============================

; ----- get best images and parameters

; NB: in tmode=2 these best images are averspan averages with
; derotation for the trimboxcoal center (averspan 6 min = 1 arcsec) so
; they have corresponding near-limb tile shifts during the averspan
; that get time-averaged

; clean on-screen plots from previous runs
spawn,'killgv'

; get GONG "best" image
gongim=readfits('gong/temp/nispbest.fits')

; get GONG parameters
readcol,'gong/temp/nisp_coalign.dat',px2,angle,xpx2asym,ypx2asym,$
  xshiftfor,yshiftfor,xshiftrev,yshiftrev,nxfor,nyfor,nxrev,nyrev
px_gong=px2[0]
angle_gong=angle[0]

; get GONG-corresponding "best" AIA 304 image
gongaiaim=readfits('gong/temp/aiabest.fits')

; get SOLIS "best" image
solisim=readfits('solis/hacw/temp/nispbest.fits')

; get SOLIS parameters
readcol,'solis/hacw/temp/nisp_coalign.dat',px2,angle,xpx2asym,ypx2asym,$
  xshiftfor,yshiftfor,xshiftrev,yshiftrev,nxfor,nyfor,nxrev,nyrev
px_solis=px2[0]
angle_solis=angle[0]

; get SOLIS-corresponding "best" AIA 304 image
solisaiaim=readfits('solis/hacw//temp/aiabest.fits')

; use findalignimages for trimboxcenter (two images must get same px scale)
print,' -----'
print,' ----- start findalignimages trimboxim2center'
print,' -----'

if (combo eq 1) then begin
  wavpair=['solis','gong'] ; lower first
  im1=gongim   ; blink A = GONG
  im2=solisim
  px1=px_gong
  px2=px_solis ; the lower one
  angle=+(angle_solis-angle_gong)
endif
if (combo eq 2) then begin
  wavpair=['gong','aia304']
  im1=gongaiaim ; blink A = 304
  im2=gongim
  px1=0.6
  px2=px_gong
  angle=angle_gong
endif
if (combo eq 3) then begin
  wavpair=['solis','aia304']
  im1=solisaiaim   ; blink A = 304
  im2=solisim
  px1=0.6
  px2=px_solis
  angle=angle_solis
endif

; boxsize above > trimboxim2
boxpx=fix(boxsize/px2)
szim2=size(im2)
nx_im2=szim2[1]
ny_im2=szim2[2]
trimboxim2=[nx_im2/2-boxpx/2,ny_im2/2-boxpx/2,$
            nx_im2/2+boxpx/2,ny_im2/2+boxpx/2]

findalignimages,im1,im2,px1,px2,angle,px2asym,$
  shiftfor,shiftrev,nxfor,nyfor,nxrev,nyrev,$    
  nxmuckfor=nxmuckfor,nymuckfor=nymuckfor,$
  nxmuckrev=nxmuckrev,nymuckrev=nymuckrev,$
  skipmetcalf=skipmetcalf,maxmetcalf=maxmetcalf,minmetqual=minmetqual,$
  applypx2asym=applypx2asym,trimboxim2=trimboxim2,inshift12=inshift12,$
  smearim1=smearim1,smearim2=smearim2,$
  histopim1=histopim1,histopim2=histopim2,$
  muckdarkim1=muckdarkim1,muckdarkim2=muckdarkim2,$
  muckbrightim1=muckbrightim1,muckbrightim2=muckbrightim2,$
  flatten=flatten,flatim1=flatim1,flatim2=flatim2,$
  metqual=metqual,nmet=nmet,finalrinse=finalrinse,$
  checkmetcalf=checkmetcalf,blink=blink,show=0,verbose=verbose

;; print,' ----- px2/px_im2 ='+trimd(px2/px_im2,4)
;; print,' ----- angle ='+trimd(angle)+$
;;   ' angle_im2-angle_gong ='+trimd(angle_im2-angle_gong) 

; update scales to final values (@@@@@ demo = copy for pxscale12)
px2d=px2*px2asym          ; refresh 2D version to newest Metcalfed px2
pxratio=float(px1)/px2    ; refresh to newest Metcalfed px2, >1 for 1>2
pxrat2d=float(px1)/px2d 
pxscale1=[1.0,1.0]
pxscale2=[1.0,1.0]
if (pxratio lt 1.0) then begin ; px1 < px2 (non-standard CHROMIS > CRISP)
  if (applypx2asym eq 0) then pxscale2=1./[pxratio,pxratio] $
  else pxscale2=1./pxrat2d     ; pxscale1 remains [1,1]
endif
if (pxratio ge 1.0) then begin ; px1 > px2 (standard SDO > SST)
  pxscale1=[pxratio,pxratio]  
  if (applypx2asym eq 1) then pxscale2=px2asym
endif

; start showex blink nonmucked reverse (@@@@@ demo = copy content) 
;   reform im2 with latest values but no shift; instead shift im1
;   to center the rotated im2rev within a larger im1 field
reformimage,im1,im1rev,congridfactor=pxscale1,$
  shift=-shiftrev,$ 
  nxlarge=nxrev,nylarge=nyrev,missingvalue=-1,$ 
  cutcentralx=nxrev,cutcentraly=nyrev,splinip=1
reformimage,im2,im2rev,congridfactor=pxscale2,$
  rotate=-angle,$
  nxlarge=nxrev,nylarge=nyrev,missingvalue=-1,$
  cutcentralx=nxrev,cutcentraly=nyrev,splinip=1
if (show eq 1 or show eq 3) then showex,im1rev,im2rev,/blink

; make tilefigs 
shift=$
  findimshift_tiled(im2rev,im1rev,$
                    outdir='/tmp/',timediff=0.,heightdiff=heightdiff,$
                    arcsecpx=px2,pxtile=pxtile,tilecut=tilecut,$
                    precision=precision,$
                    writetileshifts=1,wavpair=wavpair,$
                    datetime=datetime,$
                    shiftchart=1,minleft=minleft,blink=0,verbose=1)

;RR STOP not needed since this is .idl

end

