filreg = '../FITS/Prom_09jul00_08[345]* ../FITS/Prom_09jul00_09*'
f = findfile(filreg)
tim = 3600l*strmid(f, 21, 2)+60l*strmid(f, 23, 2)+strmid(f, 25, 2)
tim0 = tim(0)
tim = tim-tim0

  ;;; align the frames, equalize intensity

img = 0

n = n_elements(f)
sh = intarr(2, n)
ss = intarr(2, n)

rdfits, p, f(0), /swap
p = reverse(reform(p(*, *, img)))
ss(*, 0) = (size(p))(1:2)

FOR i=1, n-1 DO BEGIN
    rdfits, p1, f(i), /swap
    p1 = reverse(reform(p1(*, *, 2)))
    ss(*, i) = (size(p1))(1:2)
    sh(*, i) = shc(p(0:255, 0:127), p1(0:255, 0:127), /filt)
;    sh(*, i) = sh(*, i)+sh(*, i-1)
    tvscl, p1
    print, i, sh(*, i)
;    p = p1
ENDFOR

sx = max(sh(0, *))-sh(0, *)
sy = max(sh(1, *))-sh(1, *)

mov = intarr(256, 128, n)
FOR i=0, n-1 DO BEGIN
    rdfits, p, f(i), /swap
    p = reverse(reform(p(*, *, 0)))
    tmp = p(sx(i):sx(i)+255, sy(i):sy(i)+127)
    mov(*, *, i) = fix(1000.*tmp/avg(tmp))
ENDFOR

ref = total(mov, 3)/n

FOR i=0, n-1 DO BEGIN
    sh(*, i) = sh(*, i)+shc(ref, mov(*, *, i), /filt)
ENDFOR

sx = max(sh(0, *))-sh(0, *)
sy = max(sh(1, *))-sh(1, *)

resx = min(ss(0, *)-sx)
resy = min(ss(1, *)-sy)

IF resx/2 NE resx/2. THEN resx = resx-1
IF resy/2 NE resy/2. THEN resy = resy-1

erase

mov = intarr(resx, resy, n)

FOR i=0, n-1 DO BEGIN
    rdfits, p, f(i), /swap

    p = reverse(reform(p(*, *, 0)))

    tmp = p(sx(i):sx(i)+resx-1, sy(i):sy(i)+resy-1)
    mov(*, *, i) = fix(1000.*tmp/avg(tmp))
    tvscl, mov(*, *, i)
ENDFOR

END















mkmov, mov, filreg

  ;;; equalize power spectra

enhance, mov

  ;;; remove distortions

rubb, mov, av=5, cell=24, step=16, range=3, /plot, /same

  ;;; replace one bad frame by mean of suroounding ones

mov(*, *, 28) = (mov(*, *, 29)+mov(*, *, 27))/2

openw, 1, 'mov_668x472x37.raw'
writeu, 1, mov
close, 1

s = size(mov)
sx = s(1)
sy = s(2)
num = s(3)

delvar, mov

conefilt, 'mov_668x472x37.raw', sx, sy, num, 420

mov = intarr(sx, sy, num)
openr, 1, 'sum.dat'
readu, 1, mov
close, 1

ssx = 16*(sx/16)
ssy = 16*(sy/16)

mm1 = scale_array(mov(0:ssx-1, 0:ssy-1, *))

tt = tim+tim0

hr = tt/3600
mn = (tt MOD 3600)/60
sc = tt MOD 60
tt = nnumber(hr, 2)+":"+nnumber(mn, 2)+":"+nnumber(sc, 2)

show, mm1(*, *, 0), /cent

FOR i=0, num-1 DO BEGIN
    tv, mm1(*, *, i)
    clock, tt(i), /dev, siz=100, pos=[10, 10]
    mm1(*, *, i) = tvrd()
ENDFOR

write_mpeg, '../MPEG/AR9169_20sep00_sp.mpg', mm1

END