function [z, zmax] = psf(pupil,w,zmax)


z = pupil.*exp(2*pi*j*w);
z = fftshift(z);
z = fft2(z);
z = fftshift(z);
z = z.*conj(z);
if (nargin<3)
    zmax = max(max(z));
end
z = z/zmax;
if (nargin>2)
    zmax = max(max(z));
end