pro psfiter, n=n

compile_opt idl2

!except = 0
!quiet = 1

if n_elements(n) eq 0 then n = 10

dw
device, retain=2

print, 'Reading files'
data = float(readfits('images/data.movie'))
psfn = float(readfits('images/data0002.reg.sum0500.fits'))

print, 'Making arrays'
psf = make_array(/float, 64, 64, n+1, value=0)
a = make_array(/float, 64, 64, 256, value=0)
p = make_array(/float, 64, 64, 256, value=0)

psf[*,*,0] = float(readfits('images/data.reg.sum0256.fits'))

window, 0, xs=512, ys=512, title='"Best" psf'
tvscl, rebin(psfn, 512, 512, /sam)
window, 1, xs=512, ys=512, title='Current psf'
tvscl, rebin(psf[*,*,0], 512, 512, /sam)

for i=0,n-1 do begin
	print, 'i = ', i
	
	print, '  Calculating a'
	for j=0,255 do begin
		a[*,*,j] = lucy(data[*,*,j], psf[*,*,i])
	endfor
	print, '    min(a) = ', min(a)

	print, '  Calculating p'
	for j=0,255 do begin
		p[*,*,j] = lucy(data[*,*,j], a[*,*,j])
	endfor
	
	print, '  Calculating psf'
	psf[*,*,i + 1] = total(p, 3)
	psf[*,*,i + 1] = psf[*,*,i + 1] / total(psf[*,*,i + 1])
	tvscl, rebin(psf[*,*,i + 1], 512, 512, /sam)
endfor

print, 'Computing autocorrelation ensemble'
f = make_array(/float, 64, 64, 256, value=0)
for i=0,255 do begin
	f[*,*,i]=(abs(fft(data[*,*,i])))^2
endfor
f2 = total(f, 3)
f2 = shift(f2, 31, 31)
f2 = f2 / total(f2)

window, 2, xs=512, ys=512, title='ac ensemble'
tvscl, rebin(f2, 512, 512, /sam)

stop, 'psf and ac displayed, .c to continue'

dw
for i=0,n do begin
	window, i, xs=320, ys=320, title='PSF ' + string(i)
	tvscl, rebin(psf[*,*,i], 320, 320, /sam)
endfor

stop, 'psf sequence displayed, .c to continue'

dw
window, 0, xs=512, ys=512, title='data[0]'
tvscl, rebin(data[*,*,0], 512, 512, /sam)
window, 1, xs=512, ys=512, title='a[0]'
tvscl, rebin(a[*,*,0], 512, 512, /sam)

stop, 'All done, .c to return'

!except = 1
!quiet = 0

end