function [Uxy,x,y] = FieldPropagation2023(Uin,DeltaZeta,z,wavelength)
%FieldPropagation simulates the propagation of an initial Field Distribution
%Uin in space. It computes the electrical field distribution Uxy in a
%plane with coordinates x and y. This plane is placed at distance behind an
%aperture which is illuminated with light of wavelenth wavelength.

%input arguments:
%   Uin(complex double): initial field distribution:  NxN (N= 2^p)
%   DeltaZeta (double): distance between pixels in Uin in meters
%   z (double): propagation distance in meters
%   wavelength (double): in meters

%output arguments:
%   Uxy(NxN complex double): electrical filed distribution in x-y-plane after propagation
%   x(NxN double): x coordinate of x-y-plane
%   y(NxN double): y coordinate of x-y-plane

% wavenumber
k = 2*pi/wavelength;
N = length(Uin);

%Scale the aperture plane with coordinates Zeta and Eta such that the 1 mm slit occupy the central 50 points
[Zeta,Eta] = meshgrid(-(N-1)/2*DeltaZeta:DeltaZeta:(N-1)/2*DeltaZeta);

%scale the coordinates in the frequency domain
[fZeta,fEta] = meshgrid(linspace(-1/(2*DeltaZeta),1/(2*DeltaZeta),N));

%define the critical distance to choose between convolution and fourier procedure
zc = N*DeltaZeta^2/wavelength;

if z <= zc
    %convolution procedure for distances smaller than the critical distance (lecture notes, slide 171)
   
    %take fourier transform of initial field distribution and shift the zero frequency to the center
    Aperture = fftshift(fft2(Uin));
    %multiply with the phase function in frequency domain
    F = Aperture.*exp(-1i*2*pi^2/k*z*(fZeta.^2 + fEta.^2));
    %take the inverse fourier transform to obtain field distribution
    Uxy = ifft2(F);
    %spatial coordinates remain unchanged
    x = Zeta;
    y = Eta;
else
    %fourier procedure for distances larger than the critical distance (lecture notes, slide 172)
    
    %multiply initial distribution with phase function in space
    f = Uin.*exp(1i*k/(2*z)*(Zeta.^2 + Eta.^2));
    %take fourier transform and shift the zero frequency to the center
    F = fftshift(fft2(f));
    %scale the spatial coordinates (x=y=fZeta*wavelength*z)
    x = fZeta*wavelength*z;
    y = fEta*wavelength*z;
    %multiply with phase function to obtain the field distribution
    Uxy = 1/(1i*wavelength*z)*F.*exp(1i*k/(2*z)*(x.^2 + y.^2));
end
end