FUNCTION Readframe, unit, n
COMMON assocsize, sx0, sy0, imsiz
img = intarr(sx0, sy0)
point_lun, unit, 1l*imsiz*n
readu, unit, img
return, img
END

PRO Writeframe, img, unit, n
COMMON assocsize, sx0, sy0, imsiz
tmp = intarr(sx0, sy0)
s = size(img, /dim) < [sx0, sy0]
tmp(0:s(0)-1, 0:s(1)-1) = img(0:s(0)-1, 0:s(1)-1)
point_lun, unit, 1l*imsiz*n
writeu, unit, tmp
END


COMMON assocsize, sx0, sy0, imsiz

base = 'Plage-*'
fg = findfile('../FITS/GBand/'+base)
fb = findfile('../FITS/BlueCont/'+base)
fr = findfile('../FITS/Halpha/'+base)
ft = findfile('../FITS/RedCont/'+base)
fc = findfile('../FITS//CaIIH/'+base)
  ;;; ca and ha may be switched -> names change; sort on date/time
fc = fc(sort(strmid(fc, strlen(fc(0))-20)))
fr = fr(sort(strmid(fr, strlen(fr(0))-20)))

 ;;; jumps?  only use part of the files
first = 0
last = n_elements(fg)-1

fg = fg(first:last)
fb = fb(first:last)
fr = fr(first:last)
ft = ft(first:last)
fc = fc(first:last)

num = n_elements(fg)

;;; check timeline
s = strlen(fg(0))
tim = strmid(fg, s-11, 2)*3600l+strmid(fg, s-9, 2)*60l+strmid(fg, s-7, 2)
c = poly_fit(indgen(n_elements(tim)), tim, 1)
print, c(1)

 ;;; # wavelength points in Ha
n_ha = 5
 ;;; additional shift between RC and HA?
sh_ha = [0, 0]
xr0 = abs(sh_ha(0) > 0)
yr0 = abs(sh_ha(1) > 0)
xh0 = abs(sh_ha(0) < 0)
yh0 = abs(sh_ha(1) < 0)

 ;;; name settings for the cubes
nam = 'Plage-'
dat = '-20060424'
lam = ['gb-0000', 'bc-0000', 'rc-0000', 'ha-0700', 'ha-0350', 'ha-0000', $
       'ha+0350', 'ha+0700', 'ca-0000']
dt = '-dt30'

pf = 0   ;;; polyfit subtract for shc?

par_ca = [-0.343, 1.0040]
par_bc = [0.042, 1.0020]
par_rc = [-0.966, 1.0048]
par_ha = par_rc

;;; (for i in */*fits; do fits $i |grep NAXIS1|cut -c 27-30;done)|sort|tail -5

sx0 = 1236 & sy0=966
imsiz = 2l*sx0*sy0

FOR i=1, n_ha+3 DO openw, i, 'cube-'+strtrim(i, 2)+'.asc'

sall = intarr(2, num)
nsum = 3
sh = intarr(2, num)
scl = 5000
scl_ha = 25
scl_ca = 3000
imax = 32767

window, 1, xs=sx0, ys=sy0

FOR i=0, num-1 DO BEGIN
    imarr = intarr(sx0, sy0, n_ha+3)
      ;;; read G band
    rdfits, p, fg(i)
    p1 = reform(p(*, *, 0))
    sgb = size(p1, /dim)
    imarr(0, 0, 0) = fix(scl*p1 < imax)
      ;;; read blue cont
    rdfits, p, fb(i)
    p1 = rot(reform(p(*, *, 0)), par_bc(0), par_bc(1), cub=-.5)
    sbc = size(p1, /dim)
    imarr(0, 0, 1) = fix(scl*p1 < imax)
      ;;; read red cont/Halpha
    rdfits, p, fr(i)
    p1 = rot(reform(p(*, *, 0)), par_rc(0), par_rc(1), cub=-.5)
    src = size(p1, /dim)
      ;;; handle remaining ha-shift
    src = src-abs(sh_ha)
    p1 = fix(scl_ha*p1(xr0:xr0+src(0)-1, yr0:yr0+src(1)-1) < imax)
    imarr(0, 0, 2) = p1
    FOR j=1, n_ha DO BEGIN
        p1 = rot(reform(p(xh0:xh0+src(0)-1, yh0:yh0+src(1)-1, j)), $
                 par_rc(0), par_rc(1), cub=-.5)
        imarr(0, 0, 2+j) = fix(scl_ha*p1 < imax)
    ENDFOR
    
      ;;; now shift all triples to match and cut the size to be the same
    
    ssh = intarr(2, 3)
    FOR j=1, 2 DO ssh(*, j)=$
      shc(imarr(*, *, 0), imarr(*, *, j), /N2, /filt, POLY=pf)
    ssh(0, *) = ssh(0, *)-max(ssh(0, *))
    ssh(1, *) = ssh(1, *)-max(ssh(1, *))
    
    FOR j=0, 1 DO imarr(*, *, j) = shift(imarr(*, *, j), ssh(0, j), ssh(1, j))
    FOR j=2, n_ha+2 DO imarr(*, *, j) = shift(imarr(*, *, j), ssh(0, 2), ssh(1, 2))
    
    sall(*, i) = (sgb+ssh(*, 0)) < (sbc+ssh(*, 1)) < (src+ssh(*, 2))
    imarr = imarr(0:sall(0, i)-1, 0:sall(1, i)-1, *)
    
    tvscl, imarr(*, *, 0) & wait, .5
    tvscl, imarr(*, *, 1) & wait, .5
    tvscl, imarr(*, *, 2) & empty
    
    FOR j=0, n_ha+2 DO writeframe, imarr(*, *, j), j+1, i
    
    p1=$
      imarr(*, *, 0)/avg(imarr(*, *, 0)) + $
      imarr(*, *, 1)/avg(imarr(*, *, 1)) + $
      imarr(*, *, 2)/avg(imarr(*, *, 2))
    
    IF i EQ 0 THEN BEGIN
        pp = fltarr(nsum, sall(0, 0), sall(1, 0))
        pp(0, *, *) = p1
    ENDIF ELSE BEGIN
          ;;; find valid size
        spp = sall(*, i)
        FOR k=(i-num) > 0, i-1 DO spp = spp < sall(*, k)
        ref = total(pp(*, 0:spp(0)-1, 0:spp(1)-1), 1)
        sh(*, i) = shc(ref, p1, /filt, /N2, POLY=pf)
        sx = (sall(0, i) < sall(0, 0))
        sy = (sall(1, i) < sall(1, 0))
        p1 = p1(0:sx-1, 0:sy-1)
        pp(i MOD nsum, 0:sx-1, 0:sy-1) = shift(p1, sh(0, i), sh(1, i))
        print, i, sh(*, i)
        tvscl, ref
    ENDELSE

ENDFOR

sx = max(sh(0, *))-sh(0, *)
sy = max(sh(1, *))-sh(1, *)
resx = min(sall(0, *)-sx)
resy = min(sall(1, *)-sy)
  ;;; mave even dimensions
resx = resx/2*2
resy = resy/2*2

FOR j=1, n_ha+3 DO BEGIN
    FOR i=0, num-1 DO BEGIN
        writeframe, shift(readframe(j, i), -sx(i), -sy(i)), j, i
    ENDFOR
ENDFOR

;;; Now Ca
;;; either separate if only core (bad correlation w. G band)
;;; or align in two steps, first wing<>GB, then core<>wing
ca_version = 0

;stop

IF ca_version EQ 0 THEN BEGIN
      ;;; Version 1
    mkmov, ca, fc, ang=par_ca#replicate(1., num), scale=scl_ca, nsum=3
    ref = float(readframe(1, num/2-5))
    FOR i=num/2-4, num/2+5 DO ref = ref+readframe(1, i)
    ref = ref(0:resx-1, 0:resy-1)
    tmp = total(ca(*, *, num/2-5:num/2+5), 3)
    print, (sh = shc(ref, tmp, /filt))
    sh = -manblinkshift(ref, tmp, ini=-sh)
    stop
    newx0 = sh(0) > 0
    newy0 = sh(1) > 0
    ca = ca(-sh(0):*, -sh(1):*, *)
ENDIF ELSE BEGIN
      ;;; Version 2
    ca = intarr(sx0, sy0, num)
    sca = intarr(2, num)
    FOR i=0, num-1 DO BEGIN
        rdfits, p, fc(i)
        p1 = rot(reform(p(*, *, 0)), par_ca(0), par_ca(1), cub=-.5)
        sca(*, i) = size(p1, /dim)
        ca(0:sca(0, i)-1, 0:sca(1, i)-1, i) = fix(scl_ca*p1 < imax)
    ENDFOR
    sh = intarr(2, num)
    
    IF ca_version EQ 1 THEN BEGIN
        FOR i=1, num-1, 2 DO BEGIN
            gb = (readframe(1, i))(0:resx-1, 0:resy-1)
            print, i, (sh(*, i)=shc(gb, ca(*, *, i), /filt, /n2, poly=pf))
        ENDFOR
        stop
        FOR i=1, num-1, 2 DO $
          ca(*, *, i) = shift(ca(*, *, i), sh(0, i), sh(1, i))
        
        tmp = sca+sh
        sh(*, 0) = shc(ca(0:tmp(0, 1), 0:tmp(1, 1), 1), $
                       ca(*, *, 0), /filt, /n2, poly=pf)
        FOR i=2, num-1, 2 DO print, i, $
          (sh(*, i)=shc(total(ca(0:tmp(0, i), 0:tmp(1, i), [i-1, i+1]), 3), $
                        ca(*, *, i), /filt, /n2, poly=pf))
        stop
        
        FOR i=0, num-1, 2 DO $
          ca(*, *, i) = shift(ca(*, *, i), sh(0, i), sh(1, i))
    ENDIF ELSE BEGIN
        FOR i=0, num-1 DO BEGIN
            gb = (readframe(1, i))(0:resx-1, 0:resy-1)
            print, i, (sh(*, i)=shc(gb, ca(*, *, i), /filt, /n2, poly=pf))
        ENDFOR
        stop
        FOR i=0, num-1 DO ca(*, *, i) = shift(ca(*, *, i), sh(0, i), sh(1, i))
    ENDELSE
    
    newx0 = max(sh(0, *)) > 0
    newy0 = max(sh(1, *)) > 0
    
    newx = min(sca(0, *)+sh(0, *))
    newy = min(sca(1, *)+sh(1, *))
    ca = (temporary(ca))(newx0:newx-1, newy0:newy-1, *)
    
    IF newx0+newy0 GT 0 THEN BEGIN
        FOR j=1, n_ha+3 DO BEGIN
            FOR i=0, num-1 DO BEGIN
                writeframe, shift(readframe(j, i), -newx0, -newy0), j, i
            ENDFOR
        ENDFOR
    ENDIF
    
ENDELSE

;;; same for all

s = size(ca, /dim)/2*2
resx = ((resx-newx0) < s(0))/2*2
resy = ((resy-newy0) < s(1))/2*2

get_lun, unit
spawn, 'mkdir ../Cubes'

FOR j=0, n_ha+2 DO BEGIN
    openw, unit, '../Cubes/'+nam+lam(j)+dat+'-cube-'+resstring([resx, resy, num])+dt
    FOR i=0, num-1 DO writeu, unit, (readframe(j+1, i))(0:resx-1, 0:resy-1)
    close, unit
ENDFOR

openw, unit, '../Cubes/'+nam+lam(n_ha+3)+dat+'-cube-'+resstring([resx, resy, num])+dt
FOR i=0, num-1 DO writeu, unit, ca(0:resx-1, 0:resy-1, i)
close, unit

free_lun, unit
close, /all

END