高斯光束传播 matlab

%Simulation of diffraction of Gaussian Beam

clear;

%Gaussian Beam

%N:sampling number

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:step size

dx=L/N;

for n=1:N+1

    for m=1:N+1

        %Space axis

        x(m)=(m-1)*dx-L/2;

                    y(n)=(n-1)*dx-L/2;

        %Gaussian Beam in space domain

                    Gau(n,m)=exp(-(x(m)^2+y(n)^2)/(wo^2));

        %Frequency axis

        Kx(m)=(2*pi*(m-1))/(N*dx)-((2*pi*(N))/(N*dx))/2;

        Ky(n)=(2*pi*(n-1))/(N*dx)-((2*pi*(N))/(N*dx))/2;

        %Free space transfer function

        H(n,m)=exp(j/(2*ko)*z*(Kx(m)^2+Ky(n)^2));

    end

end

%Gaussian Beam in Frequency domain

FGau=fft2(Gau);

FGau=fftshift(FGau);

%Propagated Gaussian beam in Frequency domain

FGau_pro=FGau.*H;

%Peak amplitude of the initial Gaussian beam

Peak_ini=max(max(abs(Gau)));

sprintf('Initial peak amplitude is %f [mm]',Peak_ini)

%Propagated Gaussian beam in space domain

Gau_pro=ifft2(FGau_pro);

Gau_pro=Gau_pro;

%Peak amplitude of the propagated Gaussian beam

Peak_pro=max(max(abs(Gau_pro)));

sprintf('Propagated peak amplitude is %f [mm]',Peak_pro)

%Calculated Beam Width

[N M]=min(abs(x));

Gau_pro1=Gau_pro(:,M);

[N1 M1]=min(abs(abs(Gau_pro1)-abs(exp(-1)*Peak_pro)));

Bw=dx*abs(M1-M)*10^3;

sprintf('Beam width(numerical) is %f[mm]',Bw)

%Theoretical Beam Width

W=(2*z_ray)/ko*(1+(z/z_ray)^2);

W=(W^0.5)*10^3;

sprintf('Beam width(theoretical) is %f[mm]',W)

%axis in mm scale

x=x*10^3;

y=y*10^3;

figure(1);

mesh(x,y,abs(Gau))

title('Initial Gaussian Beam')

xlabel('x [mm]')

ylabel('y [mm]')

axis([min(x) max(x) min(y) max(y) 0 1])

axis square

figure(2);

mesh(x,y,abs(Gau_pro))

title('propagated Gaussian Beam')

xlabel('x [mm]')

ylabel('y [mm]')

axis([min(x) max(x) min(y) max(y) 0 1])

axis square