题目：根据高斯光束数学模型，模拟仿真高斯光束在谐振腔中某一位置处的归一化强度分布并给出其二维、三维强度分布仿真图；用Matlab读取实际激光光斑照片中所记录的强度数据（读取照片中光斑的一个直径所记录的强度数据即可，Matlab读取照片数据命令为imread）,用该数据画出图片中激光光斑的强度二维分布图，与之前数学模型仿真图对比。（如同时考虑高斯光束光斑有效截面半径和等相位面特点，仿真高斯光束光强、光斑有效截面半径以及等相位面同时随传播距离z的变化并给出整体仿真图可酌情加分。）
原始光斑如图1所示，用imread命令读入matlab后直接用imshow命令读取即可，
CCD采集的高斯光束光强分布
图1 CCD采集的高斯光束强度分布
读入的数据是一个224 X 244的矩阵，矩阵中的数值代表光强分布。用读入的数据取中间一行（122行）画出强度分布如图2所示。
50
100150200
020406080100120140160
180实验测量高斯曲线
图2 实验测量高斯曲线
用理论上的高斯曲线公式画出理论高斯曲线如图3所示。
-40
-30-20-10010203040
00.2
0.4
0.6
0.8
1
理论高斯曲线
图3 理论高斯曲线
M文件如下：
A=imread('D:\documents\作业\激光原理与应用\高斯.bmp');
A1=A(:,122);
x1=1:1:224;
x2=-100:1:100;
a2=exp(-x2.^2/10);
figure
imshow(A);
axis off
title('\fontsize{12}CCD采集的高斯光束光强分布');
figure
plot(x2,a2,'linewidth',1,'color','b');
axis([-40 40 0 ])
title('\fontsize{12}实验测量高斯曲线')
figure
plot(x1,A1,'linewidth',1,'color','r')
title('\fontsize{12}理论高斯曲线')
axis([50 200 0 180])
画三维强度分布。取图片矩阵的中间层，用mesh命令画出三维图如图4所示。
图4 三维强度分布
由于读入的图片有一行白边，需要手动去除掉，否则三维图会有一边整体竖起来，影响观察。最终的M文件如下。
A=imread('D:\documents\作业\激光原理与应用\高斯.bmp');
[high, width, color] = size(A);
x=1:width;
y=1:high-1;
mesh(x', y', double(A(2:224,:,1)));
grid on
xlabel('x'),ylabel('y'),zlabel('z');
title('三维强度分布');
再用matlab仿真理论上传播过程中高斯光束的变化
这次先给出M文件：
%
%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^*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
程序主要根据高斯光束的传播规律计算传播过程中任意z处的高斯光强分布。运行结果：
Number of samples(enter from 100 to 500)=500
wavelength of light in [micrometers]=
Waist of Gaussian Beam in [mm]=1
ans =
Rayleigh range is [mm]
Propagation length (z) in [mm]100000
ans =
Initial peak amplitude is [mm]
ans =
Propagated peak amplitude is [mm]
ans =
Beamwidth(numerical) is [mm]
ans =
Beamwidth(theoretical) is [mm]
>>
束腰半径处的理想高斯光强分布
传播1m处的理想高斯光强分布
传播10m处的理想高斯光强分布
传播20m处的理想高斯光强分布
传播30m处的理想高斯光强分布
传播50m处的理想高斯光强分布
传播100m处的理想高斯光强分布
而用实验测得的光斑仿真的结果是：
原始光斑的光强分布
处
1m处
处
5m处
10m以后，已经基本是均匀强度的光斑