I wrote a short matlab script file that is suppose to run Fresnel's propagation (diffraction), such that given a certain input field U0, it will tell you how the field looks after distance z0. I compared the result to textbook results, and it seems like my program works fine. The problem is if I try to take two propagation steps instead of just one. i.e, instead of taking a single iteration of the program to propagate a distance z0, I take two iterations of the program to propagate a distance z0/2 each. Then I get complete nonsense, and I can't figure out what's the problem. Any advice will be accepted with great gratitude. Here is the code:
function U = fresnel_advance (U0, dx, dy, z, lambda)
% The function receives a field U0 at wavelength lambda
% and returns the field U after distance z, using the Fresnel
% approximation. dx, dy, are spatial resolution.
k=2*pi/lambda;
[ny, nx] = size(U0);
Lx = dx * nx;
Ly = dy * ny;
dfx = 1./Lx;
dfy = 1./Ly;
u = ones(nx,1)*((1:nx)-nx/2)*dfx;
v = ((1:ny)-ny/2)'*ones(1,ny)*dfy;
O = fftshift(fft2(U0));
H = exp(1i*k*z).*exp(-1i*pi*lambda*z*(u.^2+v.^2));
U = ifft2(O.*H);
After calling fft2, you call also fftshift to have the DC frequency at the middle.
But when you call ifft2, the function assumes that you have still the DC frequency at (1,1). So you have to come back to this format before doing the inverse FFT.
So changing the final line to U = ifft2(fftshift(O.*H)) may solve the problem.
EDIT
I've just seen that Matlab advises to use ifftshift after fftshift insteaf of twice ifftshift (can't find the version in which it is introduced). According to the documentation, the sequences of calls ifftshift(fftshift(X)) and ifftshift(fftshift(X)) are not equivalent in case of odd sizes.
So I think it is better to do: U = ifft2(ifftshift(O.*H)) in the last step of your code.