changing circular shift in hologram reconstruction
The problem is this line:

recon_image = ifftshift(ifft2(psi));


It should be

recon_image = ifft2(ifftshift(psi));


You need to shift the 0 frequency back to the corners before doing ifft2

I find fftshift / ifftshift and creating the frequency coordinates causes so many programming errors. This is my code for creating the coordinates (with 0 frequency in the corner):

function [ux,uy,r,th] = fftMesh(sizeV)
%FFTMESH - Create DFT cartesian and polar coordinates
%
% Inputs:
% 
% sizeV     [width,height] of spectrum (same as image size)
%
% Outputs:
%
% ux        x in cartesian coords
% uy        y in cartesian coords
% r         radius in polar coords  
% th        angle (rad CCW) in polar coords

rows = sizeV(1);
cols = sizeV(2);

[ux, uy] = meshgrid(...
    ([1:cols]-(fix(cols/2)+1))/(cols-mod(cols,2)), ...
    ([1:rows]-(fix(rows/2)+1))/(rows-mod(rows,2))...
    );
ux = ifftshift(ux);   % Quadrant shift to put 0 frequency at the corners
uy = ifftshift(uy);

% Convert to polar coordinates
th = atan2(uy,ux);
r = sqrt(ux.^2 + uy.^2);

end


If you use this, you dont have to worry about using fftshift / ifftshift. Alternatively, just remove the ifftshift lines in the function to keep the 0 freq in the middle.

Updated code:

clear;

lambda = 632.8e-9/1.333; %wavelength
pix_size = 7.4e-6; %camera's pixel size
dx = pix_size;
dy = pix_size;

hologram = double(rgb2gray(imread('test.jpg')));
N = size(hologram);

tukey_win = tukeywin(N(1),0.075); %tukey window
% multiply image (single precision) with the tukey window to darken edges
trimmed_holo = single(hologram) .* (tukey_win * tukey_win'); 

prop_distance = 0.01:0.005:0.2; % propagation distances vector
fft_hologram = fft2(trimmed_holo); %Fourier transform of hologram

for k = 1:length(prop_distance) % do multiple propagation distances
% obtain transfer function 
G = trans_func(lambda, N(1), dx, dy, prop_distance(k));

%multiply FT of impulse response times FT of hologram 
psi = G .* fft_hologram; 

% get slice of reconstructed optical field through inverse FT
recon_image = ifft2(psi); 

%write image to file
imagesc(real(recon_image)); pause;
imagesc(imag(recon_image)); pause;
imagesc(abs(recon_image))

end

function G = trans_func(lambda, N, Delta_xi, Delta_eta, z)

M = N;
[m, n] = meshgrid(...
    ([1:M]-(fix(M/2)+1))/(M-mod(M,2)), ...
    ([1:N]-(fix(N/2)+1))/(N-mod(N,2))...
    );
m = ifftshift(m);   % Quadrant shift to put 0 frequency at the corners
n = ifftshift(n);

first_term  = (lambda^2 * (n + (N^2 * Delta_xi^2)  ./ (2 * z * lambda)).^2) ...
    ./ (N^2 * Delta_xi^2);
second_term = (lambda^2 * (m + (M^2 * Delta_eta^2) ./ (2 * z * lambda)).^2) ...
    ./ (M^2 * Delta_eta^2);
arg_root = 1 - first_term - second_term;
arg_root(arg_root&lt;0) = 0;