; file: gb_ha_align.pro
; init: Oct 10 2011 
; note: not a .pro but the log of what I did 
;       needed consecutive alignments gb-rc, rc-ha5; then applied to hac 
; use xssw

; read gb image
rdfits,algb,'../fits/gb-ca2gb.fits' ; product from preceding alignment gb-ca

; read ha images
datapath='/home/rutten/data/DOT/2010-09-27/FITS/'
rdfits,ha,datapath+'E09N16mosS-ha-scan-20100927-135439.fits',header=haheader
rc=reform(ha[*,*,0])
hac=reform(ha[*,*,1])  ; size [850,850]
ha2=reform(ha[*,*,2])
ha3=reform(ha[*,*,3])
ha4=reform(ha[*,*,4])
ha5=reform(ha[*,*,5])

; coarse scaling  lowerleft [0,0] seems similar features
;    gb y=700 ~ ha y=500  ratio ha-px/gb-px = 0.71
; officially gb px 0.071 arcsec, ha px 0.109 arcsec, so OK, apply
oldn=850  ; for x and y; why not 2K/2?  strange 
newn=fix(oldn*0.109/0.071+0.5)
reshac=congrid(hac,newn,newn,cubic=-0.5)
resha4=congrid(ha4,newn,newn,cubic=-0.5)
resha5=congrid(ha5,newn,newn,cubic=-0.5)
resrc=congrid(rc,newn,newn,cubic=-0.5)

; check gb versus ha FOVs
window,0,xsi=1000,ysi=800
plot_image,gb           
window,1,xsi=800,ysi=1000
plot_image,ha4  

; first apply coarse manual shift alignment to get into ball park
window,0,xsi=1000,ysi=800
shiftrc=manblinkshift(algb,resrc)
  ; works very well! Recognize granules.  Result: shiftrc=[11,79]
shresrc=shift(resrc,-shiftrc)
flicker,algb,shresrc,1   ; indeed now mostly field rotation

; feature alignment, see Oct 10 2011 emails from Pit Suetterlin
  ref_pic=algb
  pic=shresrc
;  delvar, ref_points, f_points    ;; do this at restart 
;  .run find_points                ;; click manually on suited pairs
  fa={X:float(ref_points),Y:float(f_points)}
  p=[589., 453., 0., 1., 0., 0.]  ;; first pair = (rough) center as pivot
  pi = replicate({fixed:0b, limited: [1b, 1b], limits: [1.,1.]},6)
  pi(0:1).fixed=1
; set limits to search space: .18 = 10 deg angle; 5% in scale, 30 px
  pi.limits(0)=[589.,453.,-.18,0.95,-30,-30]  ; same pivot coordinates
  pi.limits(1)=[589.,453.,.18,1.05,30,30]
  p1=mpfit('rsf',p,functargs=fa,parinfo=pi,/quiet)
; check
  print,p1
; 589.000      453.000    0.0455231      1.01899      1.19042      4.74012
; p1=[589.000,453.000,0.0455231,1.01899,1.19042,4.74012]
    ; most luckily I copied the output here; always do, then no repeat
    ; selecion needed if you loose the resulting image (as I did)
  plot,f_points[0,*],f_points[1,*],psy=1
  plots,rsff(p1,ref_points),psy=4
; this is how to shift reference to pic, so reverse the result in applying
  pic1 = shift(rot(pic,-p1[2]*!radeg, 1./p1[3],p1[0],p1[1],/pivot),$
     round(-p1[4]),round(-p1[5]))
; check with flicker
  window,xsi=1000,ysi=800
  flicker,ref_pic,pic1,1  ; OK! minute shift, forget for now

; apply to ha4
shresha4=shift(resha4,-shiftrc)
alha4=shift(rot(shresha4,-p1[2]*!radeg, 1./p1[3],p1[0],p1[1],/pivot),$
     round(-p1[4]),round(-p1[5]))
flicker,algb,alha4,0.5  ; not so good

; apply to ha5 = shows granules better (red wing)
shresha5=shift(resha5,-shiftrc)
alha5=shift(rot(shresha5,-p1[2]*!radeg, 1./p1[3],p1[0],p1[1],/pivot),$
     round(-p1[4]),round(-p1[5]))
flicker,algb,alha5,0.5   ; small granule shift, can be done coorelation?

; find shift gb-ha5
shiftha5=shc(algb,alha5,/N2,/filt)   ; -7, 5
; trial=shift(alha5,-shiftha5)
; flicker,algb,trial,0.5  ; no good so filments disturb too much??
; manshiftha5=manblinkshift(algb,alha5)  ; 8,-5  ;so shc did well but sign
alha5n=shift(alha5,shiftha5)
flicker,algb,alha5n,0.5 ; excellent!

; apply all to hac
shreshac=shift(reshac,-shiftrc)
alhac=shift(rot(shreshac,-p1[2]*!radeg, 1./p1[3],p1[0],p1[1],/pivot),$
     round(-p1[4]),round(-p1[5]))
alha=shift(alhac,shiftha5)
flicker,algb,alha,1  ; no features whatsoever in common but must be it

; define edges to be removed (shift + rotation)
; seems to work but what if shifts are negative or rotation clockwise?
ref_pic=algb
pic1=alha
sizeref=size(ref_pic)
nx=sizeref[1]
ny=sizeref[2]

; plot_image,alha   ; big striped triangles but outside algb field of view
rotshift=fix(sin(p1[2])*nx/2+2)
rotshift=0 ; Halpha image much bigger so no problem tilted edges

xshift=round(p1[4])+rotshift
yshift=round(p1[5])+rotshift

; apply edge removal
if (xshift gt 0) then alref=ref_pic[0:nx-xshift-1,*] $
  else alref=ref_pic[-xshift:nx-1,*]
if (yshift gt 0) then alref=alref[*,0:ny-yshift-1,*] $
  else alref=alref[-yshift:ny-1,*]
if (xshift gt 0) then alpic=pic1[0:nx-xshift-1,*] $
  else alpic=pic1[-xshift:nx-1,*]
if (yshift gt 0) then alpic=alpic[*,0:ny-yshift-1,*] $ 
  else alpic=alpic[-yshift:ny-1,*]

; check edge removal
window,xsi=1000,ysi=800
flicker,alref,alpic,1  ; lower left
erase
flicker,alref[400:*,400:*],alpic[400:*,400:*],1  ; upper right
erase
flicker,alref[0:*,400:*],alpic[0:*,400:*],1      ; upper left
erase
flicker,alref[400:*,0:*],alpic[400:*,0:*],1      ; lower right
  
; write result
algb=alref
alha=alpic
writefits,'../fits/algb.fits',algb  ; was [1179, 907] , now [1178, 902]
writefits,'../fits/alha.fits',alha  ; was [1305, 1305], now [1178, 902]

; now cut ca image to the same new size
newsize=size(alref)
rdfits,alca,'../fits/ca-ca2gb.fits' ; product from preceding alignment gb-ca
alcanew=alca[0:newsize[1]-1,0:newsize[2]-1]   ; simplistic since upper edge
writefits,'../fits/alca.fits',alcanew

; check 
window,xsi=1000,ysi=800
flicker,algb,alcanew,0.3
flicker,algb,alha,1
flicker,alcanew,alha,1

; for fun
scatcont,alref*1000,alpic,/scatter,/dist,/xmom,/ymom,$
    xtitle='G band intensity',ytitle='Halpha intensity',xlabeloff=0.7,$
    /ps,plotfilename='../figs/gb_ha_scatter.ps' 
  ; bull's eye
  ; very slight bright-bright = network presence

;