base = 'Filament*'
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 may be switched -> names change; sort on date/time
fc = fc(sort(strmid(fc, strlen(fc(0))-20)))

;;; 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)

 ;;; jumps?  only use part of the files
first = 0
last = 35

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

num = n_elements(fg)

 ;;; # wavelength points in Ha
n_ha = 5
 ;;; additional shift between RC and HA?
sh_ha = [0, -1]
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 = 'Filament'
dat = '20050919'
lam = ['-0700', '-0350', '-0000', '+0350', '+0700']
dt = '-dt35'

asc = 1  ;;; use assoced files?
pf = 0   ;;; polyfit subtract for shc?

par_bc = [-0.4269, 1.00211]
par_rc = [-0.9765, 1.00481]
par_ha = par_rc
par_ca = [-0.3487, 1.00232]  ;;; not measurable - use old value

;;; (for i in */*fits; do fits $i |grep NAXIS1|cut -c 27-30;done)|sort|tail -5
sx = 1238 & sy=991
IF keyword_set(asc) THEN BEGIN
    dummy = replicate(0, sx, sy)
    openw, 1, 'gb-cube.asc'
    gb = assoc(1, intarr(sx, sy))
;;    FOR i=0, num-1 DO gb(i) = dummy
    gb(num-1) = dummy
    openw, 2, 'bc-cube.asc'
    bc = assoc(2, intarr(sx, sy))
;;    FOR i=0, num-1 DO bc(i) = dummy
    bc(num-1) = dummy
    openw, 3, 'rc-cube.asc'
    rc = assoc(3, intarr(sx, sy))
;;    FOR i=0, num-1 DO rc(i) = dummy
    rc(num-1) = dummy
    openw, 4, 'ca-cube.asc'
    ca = assoc(4, intarr(sx, sy))
;;    FOR i=0, num-1 DO ca(i) = dummy
    ca(num-1) = dummy
    openw, 5, 'ha-cube.asc'
    ha = assoc(5, intarr(sx, sy, n_ha))
    dummy = replicate(0, sx, sy, n_ha)
;;    FOR i=0, num-1 DO ha(i) = dummy
    ha(num-1) = dummy
ENDIF ELSE BEGIN
    gb = intarr(sx, sy, num)
    bc = intarr(sx, sy, num)
    rc = intarr(sx, sy, num)
    ha = intarr(sx, sy, n_ha, num)
ENDELSE

sall = intarr(2, num)
nsum = 3
sh = intarr(2, num)
scl = 5000
scl_ha = 30
window, 1, xs=sx, ys=sy
tmp = fltarr(sx, sy)

FOR i=0, num-1 DO BEGIN
      ;;; read G band
    rdfits, p, fg(i)
    p1 = reform(p(*, *, 0))
    sgb = size(p1, /dim)
    IF keyword_set(asc) THEN BEGIN
        tmp(*) = 0
        tmp(0:sgb(0)-1, 0:sgb(1)-1) = fix(scl*p1)
        gb(i) = tmp
    ENDIF ELSE $
      gb(0:sgb(0)-1, 0:sgb(1)-1, i) = fix(scl*p1)
      ;;; read blue cont
    rdfits, p, fb(i)
    p1 = rot(reform(p(*, *, 0)), par_bc(0), par_bc(1), cub=-.5)
    sbc = size(p1, /dim)
    IF keyword_set(asc) THEN BEGIN
        tmp(*) = 0
        tmp(0:sbc(0)-1, 0:sbc(1)-1) = fix(scl*p1)
        bc(i) = tmp
    ENDIF ELSE $
      bc(0:sbc(0)-1, 0:sbc(1)-1, i) = fix(scl*p1)
      ;;; 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)
    IF keyword_set(asc) THEN BEGIN
        tmp(*) = 0
        tmp(0:src(0)-1, 0:src(1)-1) = fix(scl_ha*p1(xr0:xr0+src(0)-1, yr0:yr0+src(1)-1))
        rc(i) = tmp
    ENDIF ELSE $
      rc(0:src(0)-1, 0:src(1)-1, i) = fix(scl_ha*p1(xr0:xr0+src(0)-1, yr0:yr0+src(1)-1))
    IF keyword_set(asc) THEN BEGIN
        dummy(*) = 0
        FOR j=0, n_ha-1 DO BEGIN
            p1 = rot(reform(p(xh0:xh0+src(0)-1, yh0:yh0+src(1)-1, j+1)), $
                     par_rc(0), par_rc(1), cub=-.5)
            dummy(0:src(0)-1, 0:src(1)-1, j) = fix(p1*scl_ha)
        ENDFOR
        ha(i) = dummy 
    ENDIF ELSE BEGIN
        FOR j=0, n_ha-1 DO BEGIN
            p1 = rot(reform(p(xh0:xh0+src(0)-1, yh0:yh0+src(1)-1, j+1)), $
                     par_rc(0), par_rc(1), cub=-.5)
            ha(0:src(0)-1, 0:src(1)-1, j, i) = fix(p1*scl_ha)
        ENDFOR
    ENDELSE
    
      ;;; now shift all triples to match and cut the size to be the same
    
    ssh = intarr(2, 4)
    ssh(*, 1) = shc(gb(*, *, i), bc(*, *, i), /N2, /filt, POLY=pf)
    ssh(*, 2) = shc(gb(*, *, i), rc(*, *, i), /N2, /filt, POLY=pf)
    ssh(0, *) = ssh(0, *)-max(ssh(0, *))
    ssh(1, *) = ssh(1, *)-max(ssh(1, *))
    IF keyword_set(asc) THEN BEGIN
        gb(i) = shift(gb(*, *, i), ssh(0, 0), ssh(1, 0))
        bc(i) = shift(bc(*, *, i), ssh(0, 1), ssh(1, 1))
        rc(i) = shift(rc(*, *, i), ssh(0, 2), ssh(1, 2))
        dummy = ha(i)
        FOR j=0, n_ha-1 DO dummy(*, *, j) = shift(dummy(*, *, j), ssh(0, 2), ssh(1, 2))
        ha(i) = dummy
    ENDIF ELSE BEGIN
        gb(*, *, i) = shift(gb(*, *, i), ssh(0, 0), ssh(1, 0))
        bc(*, *, i) = shift(bc(*, *, i), ssh(0, 1), ssh(1, 1))
        rc(*, *, i) = shift(rc(*, *, i), ssh(0, 2), ssh(1, 2))
        FOR j=0, n_ha-1 DO ha(*, *, j, i) = shift(ha(*, *, j, i), ssh(0, 2), ssh(1, 2))
    ENDELSE
    
    sall(*, i) = (sgb+ssh(*, 0)) < (sbc+ssh(*, 1)) < (src+ssh(*, 2))

    tvscl, gb(*, *, i) & wait, .5
    tvscl, bc(*, *, i) & wait, .5
    tvscl, rc(*, *, i) & empty
    
    p1 = $
      gb(0:sall(0, i)-1, 0:sall(1, i)-1, i)/avg(gb(0:sall(0, i)-1, 0:sall(1, i)-1, i))+$
      bc(0:sall(0, i)-1, 0:sall(1, i)-1, i)/avg(bc(0:sall(0, i)-1, 0:sall(1, i)-1, i))+$
      rc(0:sall(0, i)-1, 0:sall(1, i)-1, i)/avg(rc(0:sall(0, i)-1, 0:sall(1, i)-1, i))
    
    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

IF keyword_set(asc) THEN BEGIN
    FOR i=0, num-1 DO BEGIN
        gb(i) = shift(gb(i), -sx(i), -sy(i))
        bc(i) = shift(bc(i), -sx(i), -sy(i))
        rc(i) = shift(rc(i), -sx(i), -sy(i))
        dummy = ha(i)
        FOR j=0, n_ha-1 DO dummy(*, *, j) = shift(dummy(*, *, j), -sx(i), -sy(i))
        ha(i) = dummy
    ENDFOR
ENDIF ELSE BEGIN
    FOR i=0, num-1 DO BEGIN
        gb(*, *, i) = shift(gb(*, *, i), -sx(i), -sy(i))
        bc(*, *, i) = shift(bc(*, *, i), -sx(i), -sy(i))
        rc(*, *, i) = shift(rc(*, *, i), -sx(i), -sy(i))
        FOR j=0, n_ha-1 DO ha(*, *, j, i) = shift(ha(*, *, j, i), -sx(i), -sy(i))
    ENDFOR
ENDELSE

;;; 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)
    ref = float(gb(0:resx-1, 0:resy-1, num/2))
    FOR i=num/2-5, num/2+5 DO ref = ref+gb(0:resx-1, 0:resy-1, i)
    tmp = total(ca(*, *, num/2-5:num/2+5), 3)
    print, (sh = shc(ref, tmp, /filt))
    sh = -manblinkshift(ref, tmp, ini=-sh)
    stop
    ca = ca(-sh(0):*, -sh(1):*, *)
ENDIF ELSE BEGIN
      ;;; Version 2
    ca = intarr(1280, 1024, 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*p1)
    ENDFOR
    sh = intarr(2, num)
    FOR i=1, num-1, 2 DO print, i, (sh(*, i)=shc(gb(*, *, i), ca(*, *, i), /filt, /n2, poly=pf))

    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, i), 0:tmp(1, i), 1), ca(*, *, 0), /filt, /n2, poly=1)
    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=1))

    stop

    FOR i=0, num-1, 2 DO ca(*, *, i) = shift(ca(*, *, i), sh(0, i), sh(1, i))

    newx = min(sca(0, *)+sh(0, *))
    newy = min(sca(1, *)+sh(1, *))
    ca = (temporary(ca))(0:resx-1, 0:resy-1, *)
ENDELSE

;;; same for all

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

get_lun, unit
spawn, 'mkdir ../Cubes'
openw, unit, '../Cubes/'+nam+'-gb-0000-'+dat+'-cube-'+resstring([resx, resy, num])+dt
FOR i=0, num-1 DO writeu, unit, gb(0:resx-1, 0:resy-1, i)
close, unit
openw, unit, '../Cubes/'+nam+'-bc-0000-'+dat+'-cube-'+resstring([resx, resy, num])+dt
FOR i=0, num-1 DO writeu, unit, bc(0:resx-1, 0:resy-1, i)
close, unit
openw, unit, '../Cubes/'+nam+'-rc-0000-'+dat+'-cube-'+resstring([resx, resy, num])+dt
FOR i=0, num-1 DO writeu, unit, rc(0:resx-1, 0:resy-1, i)
close, unit
openw, unit, '../Cubes/'+nam+'-ca-0000-'+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

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

free_lun, unit

END