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 = 'QS*'
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 = 16
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 = 2
 ;;; 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 = 'QS-'
dat = '-20060325'
lam = ['gb-0000', 'bc-0000', 'rc-0000', 'ha-0500', 'ha+0500', 'ca-flip']
dt = '-dt00'

pf = 0   ;;; polyfit subtract for shc?

par_bc = [ 0.0432, 1.00182]
par_ca = [-0.3273, 1.00421]
par_rc = [-0.9462, 1.00036]
par_ha = par_rc

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

sx0 = 1237 & sy0=975
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 = 30
scl_ca = 3000

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

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

sh = [-56, -69]
ca = ca(-sh(0):*, -sh(1):*, *)
    
;;; same for all

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

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