clear
N=input（‘Number?of?samples（enter?from?100?to?500）=‘）;
L=10*10^-3;
Ld=input（‘wavelength?of?light?in?[micrometers]=‘）;
Ld=Ld*10^-6;
ko=（2*pi）/Ld;
wo=input（‘Waist?of?Gaussian?Beam?in?[mm]=‘）;
wo=wo*10^-3;
z_ray=（ko*wo^2）/2*10^3;
sprintf（‘Rayleigh?range?is?%f?[mm]‘z_ray）
z_ray=z_ray*10^-3;
z=input（‘Propagation?length?（z）?in?[mm]=‘）;
z=z*10^-3;

dx=L/N;

for?n=1:N+1
????for?m=1:N+1
????????x（m）=（m-1）*dx-L/2;
????????y（n）=（n-1）*dx-L/2;
????????Gau（nm）=exp（-（x（m）^2+y（n）^2）/（wo^2））;
????????Kx（m）=（2*pi*（m-1））/（N*dx）-（（2*pi*（N））/（N*dx））/2;
????????Ky（m）=（2*pi*（n-1））/（N*dx）-（（2*pi*（N））/（N*dx））/2;
????????H（nm）=exp（j/（2*ko）*z*（Kx（m）^2+Ky（n）^2））;
????end
end

FGau=fft2（Gau）;
FGau=fftshift（FGau）;

FGau_pro=FGau.*H;

Peak_ini=max（max（abs（Gau）））;
sprintf（‘Initial?pea