reads,'psf1',a,nx,ny
tvwindow,0,xsize=10*nx,ysize=10*ny
a=a-median(a)
tvscl,congrid(a,10*nx,10*ny)
;
;
;
n=64
x=indgen(n)#replicate(1.,n)
y=replicate(1.,n)#indgen(n)
x=x-n/2
y=y-n/2
r=sqrt(x^2+y^2)*4./n
t=atan(x,y)
cx=r^3*cos(t)
cy=r^3*sin(t)
ax=r^2*cos(2*t)
ay=r^2*sin(2*t)
fo=r^2
pupil=(r lt 1.) and (r gt 0.4)
pix=24*0.55*3.58/2000.
mag=0.1/pix
tvwindow,1,xsize=10*nx,ysize=10*ny
input:
print,' input, cx,cy,ax,ay,fo'
read,cxx,cyy,axx,ayy,foo
i=complex(0.,1.)         
print,'computing psf'
psf=shift(abs(fft(pupil*exp(2*!pi*i*(foo*fo+cxx*cx+cyy*cy+axx*ax+ayy*ay))$
           ,1))^2,n/2,n/2)
print,'displaying psf'
psf=smooth(psf,mag)
tvscl,congrid(psf(n/2-mag*nx/2:n/2+mag*nx/2,n/2-mag*ny/2:n/2+mag*ny/2),10*nx,10*ny)
goto,input
stop
end