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 = 'Spicules-*'
fg = findfile('../FITS/GBand/'+base)
fb = findfile('../FITS/BlueCont/'+base)
fr = findfile('../FITS/RedCont/'+base)
fh = findfile('../FITS/Halpha/'+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)
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 = 0
 ;;; 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 = 'Spicules-'
dat = '-20060424'
lam = ['gb-0000', 'bc-0000', 'rc-0000', 'ha-0300', 'ca-0000']
dt = '-dt12'

pf = 0   ;;; polyfit subtract for shc?

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

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

sx0 = 1240 & sy0=978
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

stop

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)
ca = ca(-sh(0):*, -sh(1):*, *)

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

mkmov, ha, fh, ang=par_ha#replicate(1., num), scale=scl, nsum=3
tmp1 = total(ha(*, *, num/2-5:num/2+5), 3)
print, (sh=shc(tmp, tmp1, /filt))
sh = -manblinkshift(tmp, tmp1, ini=-sh)
ha = ha((-sh(0) > 0):*, (-sh(1) > 0):*, *)

newx0 = sh(0) > 0
newy0 = sh(1) > 0

ca = ca(newx0:*, newy0:*, *)
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
s = size(ha, /dim)/2*2
resx = ((resx-newx0) < s(0))/2*2
resy = ((resy-newy0) < s(1))/2*2

get_lun, unit
cubedir='../Cubes/'
spawn, 'mkdir '+cubedir

FOR j=0, n_ha+2 DO BEGIN
    openw, unit, cubedir+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, cubedir+nam+lam(n_ha+3)+dat+'-cube-'+resstring([resx, resy, num])+dt
FOR i=0, num-1 DO writeu, unit, ha(0:resx-1, 0:resy-1, i)
close, unit

openw, unit, cubedir+nam+lam(n_ha+4)+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