基于matlab的激光谐振腔光场分布模拟和分析

一．课程设计的主要任务：
1.任务总述：用计算机模拟激光谐振腔的光场分布。

2.设计要求： 1）编程语言不限；
2）腔型包括：条形腔，矩形平平腔，圆形平平腔，矩形共焦腔，圆形共焦腔，倾斜腔等。

二．我个人完成的情况：
1.已经完成的：
1）用基本的循环迭代法：模拟了条形腔，矩形平平腔，圆形平平腔，矩形共焦腔，圆形共焦腔的光场的振幅和相位分布：
2）用传输矩阵结合分离变量的方法：模拟了条形腔，矩形平平腔，矩形共焦腔的光场的振幅和相位分布。

三，基本原理：
1.一般的迭代法的基本原理：
1)基于菲涅尔衍射积分的基本原理：
设左右镜面的任意两个点P 和P ’点，光场分别为),(y x u 和)','(y x u ，θ是PP ’连线和光轴的夹角，ρ为PP ’之间的距离，则：
⎰⎰+=
-S
ik dS e
y x u ik
y x u '
)cos 1()
','(4),(θρ
π
ρ
同理：
因此，左右通过上两式可以把激光谐振腔的左右有效地联系在一起，给出一个面的初始光场分布，经过往返迭代，可以得出如下的光场分布特性： j j y x u y x u ）','(1
),(1γ
=
+ 12,(1
)','(++=
j j y x u y x u ）γ
则说明激光谐振腔达到了自再现的条件，也是镜面上的场分布的稳定性条件。

2）网格化的思想：虽然实际的腔镜面上的光场分布是连续的，但考虑到用计算机计算的离散的特性，需要把腔镜分割成网格，以网格上离散的节点的光场值去拟合实际的镜面的光场。

根据镜面的几何结构的特点，分割方法不尽相同，具如下： A.条形腔：等间距取点，（示意图略）：
B,矩形镜面：如下图左所示的方法进行等间隔分割与取点； C,圆形镜面：如下图右所示的方法进行等间距等角度离散。

3）化积分的运算为求和的思想：结果加和存于一个二维数组中，通过循环，完成每一点的求和，具体的见代码（附有详细的注释）。

4
A,条形腔：
B ，平平矩形腔：原始的菲涅尔积分公式： 式中：ρ
θL
=)cos(，2^2)^21(2)^21L y y x x +-+-=（ρ 。

C,平平圆心腔：积分公式同平平矩形腔，但是不同的是在极坐标下的方程，且： 2^2))^sin(2)sin(1(2))^cos(2)cos(1L r r r r +-+-=φϕφϕρ（
⎰⎰+=
-S
ik dS
e y x u ik
y x u )cos 1()
,(4)','(θρ
πρ
()()()2
21exp 2a
ikL a x x i u x e ik u x dx L L λ--⎡⎤'-''⎢⎥
=-⎢⎥⎣⎦
⎰⎰+=
-S ik dS e y x u ik y x u ')cos 1()','(4),(θρ
πρ
ϕrdrd ds ='（左边）或 φrdrd ds ='（右边）。

D,共焦矩形腔：dxdy L
yy xx ik y x u ikL L i y x u a a b b )2'
'exp(),())exp(()','(+-=⎰⎰--λ
E,圆形共焦腔：
θ
θθθθθλθπrdrd L
rr ik r u ikL L i r u R )2))'sin()sin()'cos()(cos('exp(),())exp(()','(1020+-=⎰⎰
2，基于传输矩阵（结合分离变量）的基本原理：
1）考虑到matlab 有着强大的矩阵和数组的计算功能，也为了克服一般循环迭代的速度低的问题，我把光的传输过程一个传输矩阵，问题就归结于求矩阵的特征向量的问题，以x 轴方向为例：设)]'(),1()...2(),1([n Ua n Ua Ua Ua -是左边腔镜面离散点上光场的x 轴分布，而且：)]'(),1()...2(),1([n Ub n Ub Ub Ub -是右边镜面上的光场分布，
不妨以三阶矩阵为例说明， ⎥⎥⎥⎦
⎤
⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡3213332312322
21131211321'''U U U T T T T T T T T T U U U ，且[]))2/()2)^(exp(()exp(L x x ik ikL L i T j i ij --⎪⎭
⎫ ⎝⎛-=λ
且
[]
[]
[]
32
13212322
21
23
22212
3
222
123
23
23222222
21212
13
3223
32
2-2)(y y y x x x y y y y y y y y y x x x x x x x x x y x y x y x j
i j i j
i
⎥⎥⎥⎦⎤⎢
⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡+⎥⎥⎥⎦⎤⎢⎢⎢⎣
⎡=-+=-⨯⨯
而上述的矩阵可以直接有x,y 列向量直接生成。

因此直接代入上式，得出传输矩阵。

再由初始额x 轴方向的光场分布往返做矩阵的乘法，完成迭代，而且此种方法还可以直接求出特征值和特征向量。

四．仿真结果：
1.一般的迭代法的仿真结果： 1）条形腔：
运行： A=rand(1,201);
for i=1:200 bar_type end 200次后得带稳定的幅度和相位分布如下：
且菲涅尔数是10，迭代次数为200次，基本上达到了稳定状态
2）矩形平平腔：
菲涅尔数目任然是10，以上由于计算速度比较慢，因此只迭代了5次，没有达到稳定状态。

3）圆形平平腔：菲涅尔数为10，迭代次数为3次。

4）矩形共焦腔：菲涅尔数为：0.5，迭代次数为6次。

5）球面共焦腔（圆心镜共焦腔）：菲涅尔数为4，迭代4次。

2，基于传数矩阵（结合分离变量法）的仿真方法：
1）条形腔传输矩阵法：
菲涅尔数为10，迭代了300次，达到了稳定自再现状态：
2）矩形平平腔的传输矩阵法：
菲涅尔数为1，迭代次数为200次，达到了稳定的状态。

3）矩形共焦腔的传输矩阵法：
菲涅尔数为0.2，迭代次数为10次，达到了自在线的状态。

五，结果的分析：
1.Fox_li的迭代法和特征向量矩阵方法的比较如下：
从我前面的仿真对比来看，两者各有优势和不足之处：
（1）Fox_Li的衍射积分法，：
1）优点：特点是思路简单，符合腔的模式形成的逻辑顺序，而且方法具有普适性，原则上对于任何几何形状的开腔都可以计算，而且还可以计算不对等，倾斜等的腔型：
2）缺点是运算量比较大，因此迭代的速度比较慢。

（2）特征向量法：
1）优点：由于整个运算的过程是基于向量和数组的，只需要得出单程传输的传输矩阵，就可以很快得出矩阵的特征值和特征向量：
2）缺点：由于要确定单程传输的矩阵，必须确定矩阵元，但是即使是规则额腔型，传输矩阵元也不易得出，因此该方法的关键就是确定每个矩阵元。

另外，还有很多的方法可以进行激光谐振腔的光场模拟，还有FFT算法（变换到频域上），有限元光束传播法，有限差分法等等。

2.不同的初始分布如：均匀分布，随机分布，三角分布等对最终稳定场分布的影响：
我每次进行模拟时，选择不用的初值分布，得出：这对于最终的光场的分布没有影响，不懂的初始状态，改变的只是从初始值到稳定值的迭代的次数，不会影响最终的结果。

3.具体的编程细节对结果的影响：
1）所取的点的多少的影响：点数取得越多，得到的曲线曲面越平滑，越接近真实的光场分布，但是点取得过多，会直接做成计算量的急剧上升，造成不能耗时过长或基本算不出来的情。

我的经验是：对于基于传输矩阵的算法（速度快很多）可以比Fox_li迭代法取点对多些，一维的条形腔可以比二维的矩阵或圆形腔的一个方向上的点数可以适当多些。

2）菲涅数的选取对于结果的影响：
原因分析：由于共焦腔左右都是凹面腔镜，对光束有一定的汇聚作用，特别是对于近轴光线，几何偏折损耗几乎为0，只有衍射损耗：而对于平平腔，只要是与光轴有夹角的光线，经过有限次数的传输后会偏折出腔外，因此，相同的菲涅尔数的情况下，平平腔比共焦腔的损耗要大些。

B.对于同样是平平腔或共焦腔，圆形比方形的损耗略大些，但是从理论上分析，方形共焦和圆形共焦腔的基的振幅分布，光斑尺寸，等相位面的曲率半径，光束发散角等都相同，这一点在仿真时也得到了体现。

3）其他参数，如波长，腔的尺寸等对场分布的影响，都可以归结到对菲涅尔数的影响，从而影响光场的分布。

4.分离变量法和直接考虑二维的情形对编程的影响：
1）直接考虑二维的情形：假设每个面上取得点数均为N*N个（矩形镜为例），则传输举证的阶数就非常高，高达N*N阶，因此相应的矩阵元数目为N^4个，数目非常大，不仅计算量非常大，而且矩阵元的对应也非常复杂：但是其优点是精度更高，更加接近真实场分布。

2）考虑分离变量法：由于对于直角坐标系，分离变量非常容易，且同样假设左右各N*N 个点，则考虑分离变量后，只需要求出x,和方向的N 阶矩阵，且x 和y 方向的矩阵式相同的，只有N^2个矩阵元，计算量比二维的小很多。

缺点是：要能分离变量，对原始的菲涅尔积分公式（二维）必须进行简化和近似，因此会引入一定的误差。

六．参考文献
[1] 徐银新.激光谐振腔模式研究.硕士学位论文，TN248,10701
[2] 程愿应，江超，王又青，胡进，李家熔.光腔模式及传输的特征向量算法.2005,5 [3] 周炳琨，高以智，陈峥嵘.激光原理.北京：国防出版社.2002.32~37.
七．设计体会
经过了一个星期左右的鏖战，我终于完成了这次激光课程设计，从开始学习matlab 中，到开始独立地编写程序，再到程序的调试，到结果的分析，到文档的撰写，虽然过程比较艰辛，但是觉得从中学到了很多东西。

我虽然没有和别人合作，但是在编写和调试的过程中，特别是对于一些matlab 指令的调用，我得到了女朋友（自动化专业，精通matlab)的帮助，在此特别表示感谢。

我最开始用最基本的循环迭代的方法，得到的结果有一个计算速度慢的特点，于是我决定寻求其他的解决方案，结合分离变量法的传输矩阵法就是我寻求到的运算速度大幅度提高的运算方法。

我通过看matlab 的指导书，对于matlab 的核心就是基于矩阵和数组的计算方法。

其基本的运算单元就是数组和矩阵，而非其他语言中的单一的变量。

我从matlab 的基本的指令开始学起，从符号运算到数值计算，从数值阵列到向量化运算，再到函数和数据的可视化（绘图），再到M 文件和函数句柄，觉得自己通过这次课程算是真正的对于matlab 入了门。

我虽然没有来得及编写GUI 图像界面，但是我追求的是用更多更优化的方法去仿真更多的腔型的场分布。

总结起来，我用两种方法编写了条形腔，矩形平面腔，圆形平面腔，矩形共焦腔和圆形共焦腔，总的m 文件时8个，总的代码行数接近一千，虽然遗憾没有时间完成GUI 图形界面，但还是觉得有很多的收获。

在编程的过程中，我还遇到了一些问题，虽然上网查了很多资料，但是目前还没有解决，就是在用传输矩阵法进行模拟时，如果迭代的次数到达一定的值，会出现NaN 的情况，在可视化的图形上表现为图形无法显示（坐标系上是空的），我查看变量是出现了inf 等情况，我上网查了一下，这是在IEEE 规定中，0/0,0-0,
∞∞∞
∞
-，等都会返回NaN(Not a Number),可能是我计算过程中出现上述的情形导致计算中断，但是更加奇怪的是：我可以求出传输矩阵的特征值和特征向量，却不能一直通过迭代得到特征向量，这是我不明白的地方，目前正在思考这个问题。

我追求的是更好的方法和更多的腔型，我目前正在努力的方向是：（1）用FFT 把空间域变换到频域进行模拟；（2）用时域有限元的方法，目前正在学习该方法；（3）能否对于一般的自由曲面对于任意的初始分布得到最终的光场分布。

这些事情可能要等到报告上交之后才有时间去完成，报告中无法体现。

八．附录：所有的程序代码及其相应的输入指令：
注意：由于我没有图形化界面，所以老师要看我的仿真结果，要按照我所给的输入指令行输入，主要是初始参数的形式（向量还是二维数组）和函数名。

1.一般的循环迭代方法：
1）条形腔：bar_type.m
源代码如下：
%对称条形谐振腔的自再现模
M=201; %条形腔左边的总的分割数目
N=201; %条形腔右边的总的分割数目
lamda=6.328e-7;
L=1e5*lamda;
a=1e3*lamda;
b=1e3*lamda;
N1=a^2/(L*lamda);
N2=b^2/(L*lamda);
k=2*pi/lamda;
co=sqrt(1i*(exp(-1i*k*L))/(lamda*L));
B=zeros(1,201);
for n=-(N-1)/2:1:(N-1)/2
x2=(2*b)*n/N;
for m=-(M-1)/2:1:(M-1)/2
x1=(2*a)*m/M;
B(n+(N+1)/2)=B(n+(N+1)/2)+exp(-1i*k*(x1-x2)^2/(2*L))*A(m+(M+1)/2) *(2*a/M)*co;
end
end
B=B/abs((max(B)));
A=zeros(1,201);
for m=-(M-1)/2:1:(M-1)/2
x1=(2*a)*m/M;
for n=-(N-1)/2:1:(N-1)/2
x2=(2*b)*n/N;
A(m+(M+1)/2)=A(m+(M+1)/2)+exp(-1i*k*(x1-x2)^2/(2*L))*B(n+(N+1)/2) *(2*b/N)*co;
end
end
A=A/abs(max(A));
x1=(-(M-1)/2:1:(M-1)/2)*(2*a/M); x2=(-(N-1)/2:1:(N-1)/2)*(2*b/N); subplot(2,2,1);
plot(x1,abs(A));
subplot(2,2,3);
plot(x1,angle(A));
subplot(2,2,2);
plot(x2,abs(B));
subplot(2,2,4);
plot(x2,angle(B));
输入的指令和得到的结果如下：
clear all;
A=rand(1,201);
tic;
for i=1:200
bar_type;
end
toc;
N1,N2
Elapsed time is 155.402364 seconds.
N1 =
9.999999999999998
N2 =
9.999999999999998
2）矩形平平腔：pin_pin_rectangular.m
源代码如下：
%平平矩形腔的自再现模
M = 31; %第一面镜的横向分割数N = 31; %第一面镜的纵向分割数E = 31; %第二面镜的横向分割数
F = 31; %第二面镜的纵向分割数
lamda = 6.328e-7;
L = 1e5*lamda;
a = 1e3*lamda;
b = 1e3*lamda;
c = 1e3*lamda;
d = 1e3*lamda;
k = 2*pi/lamda;
N1=(a*b)/(L*lamda);
N2=(c*d)/(L*lamda);
B=zeros(E,F);
for e=-(E-1)/2:1:(E-1)/2
for f=-(F-1)/2:1:(F-1)/2
x2=2*c*e/E;
y2=2*d*f/F;
for m=-(M-1)/2:1:(M-1)/2
for n=-(N-1)/2:1:(N-1)/2
x1=2*a*m/M;
y1=2*b*n/N;
[r,costheta]=zihanshu(x1,y1,x2,y2,L);
co=exp(-i*k*r)*(1+costheta)/r;
B(e+(E+1)/2,f+(F+1)/2)=B(e+(E+1)/2,f+(F+1)/2)+i*k/(4*pi)*(4*a*b/( M*N))*A(m+(M+1)/2,n+(N+1)/2)*co;
end
end%solve B(e+(E+1)/2,f+(F+1)/2)
end
end
B=B/abs(max(max(B)));
A=zeros(M,N);
for m=-(M-1)/2:1:(M-1)/2
for n=-(N-1)/2:1:(N-1)/2
x1=2*a*m/M;
y1=2*b*n/N;
for e=-(E-1)/2:1:(E-1)/2
for f=-(F-1)/2:1:(F-1)/2
x2=2*c*e/E;
y2=2*d*f/F;
[r,costheta]=zihanshu(x1,y1,x2,y2,L);
co=exp(-i*k*r)*(1+costheta)/r;
A(m+(M+1)/2,n+(N+1)/2)=A(m+(M+1)/2,n+(N+1)/2)+i*k/(4*pi)*(4*c*d/(
E*F))*B(e+(E+1)/2,f+(F+1)/2)*co;
end
end%solve B(e+(E+1)/2,f+(F+1)/2) end
end
A=A/abs(max(max(A)));
x1 = (-(M-1)/2:1:(M-1)/2)'*a/M;
y1 = (-(N-1)/2:1:(N-1)/2)'*b/N;
x2 = (-(E-1)/2:1:(E-1)/2)'*c/E;
y2 = (-(F-1)/2:1:(F-1)/2)'*d/F;
subplot(1,2,1);
mesh(x1,y1,abs(A));
subplot(1,2,2);
mesh(x2,y2,abs(B));
子函数的源代码如下：
function [r,costheta] = zihanshu( x1,y1,x2,y2,L ) %UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
r=L+((x1-x2)^2+(y1-y2)^2)/(2*L);
costheta=L/r;
end
函数的输入命令和得到的如下：
clear all;
A=rand(31);
pin_pin_rectangular
pin_pin_rectangular
pin_pin_rectangular
pin_pin_rectangular
pin_pin_rectangular
N1,N2
N1 =
9.999999999999998
N2 =
9.999999999999998
3）圆形平平腔；circle_type.m
源代码如下：
%圆形平平腔的自再现模
M = 48; %第一面镜的径向分割数
N = 48; %第一面镜的角向分割数
E = 48; %第二面镜的径向分割数
F = 48; %第二面镜的角向分割数
lamda = 6.328e-7;
L =1e5*lamda;
R1=1e3*lamda;
R2=1e3*lamda;
k=2*pi/lamda;
co1=1i*k/(4*pi);
N1=(R1^2)/(L*lamda);
N2=(R2^2)/(L*lamda);
B=zeros(E,F);
for e=1:E
for f=1:F
r2(e,f)=(e*R2)/E;
theta2(e,f)=(2*f*pi)/F;
x2(e,f)=r2(e,f)*cos(theta2(e,f));
y2(e,f)=r2(e,f)*sin(theta2(e,f));
for m=1:M
for n=1:N
r1(m,n)=(m*R1)/M;
theta1(m,n)=(2*n*pi)/N;
x1(m,n)=r1(m,n)*cos(theta1(m,n));
y1(m,n)=r1(m,n)*sin(theta1(m,n));
rou(m,n)=sqrt((x1(m,n)-x2(e,f))^2+(y1(m,n)-y2(e,f))^2+L^2);
costheta(m,n)=L/rou(m,n);
co(m,n)=exp(-1i*k*rou(m,n))*(1+costheta(m,n))/rou(m,n);
B(e,f)=B(e,f)+A(m,n)*(pi*R1^2/(M*N))*co(m,n)*co1;
end
end%solve B(e+(E+1)/2,f+(F+1)/2)
end
end
B=B/abs(max(max(B)));
A=zeros(M,N);
for n=1:N
r1(m,n)=(m*R1)/M;
theta1(m,n)=(2*n*pi)/N;
x1(m,n)=r1(m,n)*cos(theta1(m,n));
y1(m,n)=r1(m,n)*sin(theta1(m,n));
for e=1:E
for f=1:F
r2(e,f)=(e*R2)/E;
theta2(e,f)=(2*f*pi)/F;
x2(e,f)=r2(e,f)*cos(theta2(e,f));
y2(e,f)=r2(e,f)*sin(theta2(e,f));
rou(e,f)=sqrt((x1(m,n)-x2(e,f))^2+(y1(m,n)-y2(e,f))^2+L^2);
costheta(e,f)=L/rou(e,f);
co(e,f)=exp(-1i*k*rou(e,f))*(1+costheta(e,f))/rou(e,f);
A(m,n)=A(m,n)+B(e,f)*(pi*R2^2/(E*F))*co(e,f)*co1;
end
end%solve B(e+(E+1)/2,f+(F+1)/2)
end
end
A=A/abs(max(max(A)));
r1=(1:1:M)*(R1/M);
theta1=2*pi/N:2*pi/N:2*pi;
r2=(1:1:E)*(R2/E);
theta2=2*pi/F:2*pi/F:2*pi;
[Rr1,THETA1]=meshgrid(r1,theta1);
[Rr2,THETA2]=meshgrid(r2,theta2);
Z1=abs(A)';
Z2=abs(B)';
[x1,y1,z1]=pol2cart(THETA1,Rr1,Z1);
[x2,y2,z2]=pol2cart(THETA2,Rr2,Z2);
subplot(1,2,1);
mesh(x1,y1,z1);
subplot(1,2,2);
mesh(x2,y2,z2);
输入的指令和得到的结果如下：
clear all;
A=rand(48);
circle_type;
circle_type
N1,N2
N1 =
9.999999999999998
N2 =
9.999999999999998
4）矩形共焦腔：rec_confocs.m
源代码如下：
%方形镜共焦腔的自再现模
M = 21; %第一面镜的横向分割数N = 21; %第一面镜的纵向分割数
E = 21; %第二面镜的横向分割数
F = 21; %第二面镜的纵向分割数
lamda = 6.327e-7;
L = 2e6*lamda;
R1 = L;
R2 = L;
a = 1e3*lamda;
b = 1e3*lamda;
c = 1e3*lamda;
d = 1e3*lamda;
k = 2*pi/lamda;
cons=1i*exp(-1i*k*L)/(L*lamda);
N1=(a*b)/(L*lamda);
N2=(c*d)/(L*lamda);
B=zeros(E,F);
for e=-(E-1)/2:1:(E-1)/2
for f=-(F-1)/2:1:(F-1)/2
x2=(2*c)*(e/E);
y2=(2*d)*(f/F);
for m=-(M-1)/2:1:(M-1)/2
for n=-(N-1)/2:1:(N-1)/2
x1=(2*a)*(m/M);
y1=(2*b)*(n/N);
co2=exp(1i*k*(x1*x2+y1*y2)/L);
B(e+(E+1)/2,f+(F+1)/2)=B(e+(E+1)/2,f+(F+1)/2)+(4*a*b/(M*N))*A(m+( M+1)/2,n+(N+1)/2)*co2*cons;
end
end%solve B(e+(E+1)/2,f+(F+1)/2)
end
end
B=B/abs(max(max(B)));
A=zeros(M,N);
for m=-(M-1)/2:1:(M-1)/2
for n=-(N-1)/2:1:(N-1)/2
x1=(2*a)*(m/M);
y1=(2*b)*(n/N);
for e=-(E-1)/2:1:(E-1)/2
for f=-(F-1)/2:1:(F-1)/2
x2=(2*c)*(e/E);
y2=(2*d)*(f/F);
co1=exp(1i*k*(x1*x2+y1*y2)/L);
A(m+(M+1)/2,n+(N+1)/2)=A(m+(M+1)/2,n+(N+1)/2)+(4*e*f/(E*F))*B(e+( E+1)/2,f+(F+1)/2)*co1*cons;
end
end%solve B(e+(E+1)/2,f+(F+1)/2)
end
end
A=A/abs(max(max(A)));
x1 = (-(M-1)/2:1:(M-1)/2)'*2*a/M;
y1 = (-(N-1)/2:1:(N-1)/2)'*2*b/N;
x2 = (-(E-1)/2:1:(E-1)/2)'*2*c/E;
y2 = (-(F-1)/2:1:(F-1)/2)'*2*d/F;
subplot(1,2,1);
mesh(x1,y1,abs(A));
subplot(1,2,2);
mesh(x2,y2,abs(B));
输入的指令及其所得到的的结果如下：
clear all;
A=rand(31);
rec_confocs
rec_confocs
rec_confocs
rec_confocs
rec_confocs
rec_confocs
N1,N2
N1 =
0.500000000000000
N2 =
0.500000000000000
5）圆形共焦腔：
源代码如下: sphe_type.m
%共焦球面腔自再现模
M = 48; %第一面镜的径向分割数
N = 48; %第一面镜的角向分割数
E = 48; %第二面镜的径向分割数
F = 48; %第二面镜的角向分割数
lamda = 6.328e-7; %激光器的谐振波长
L =1e6*lamda; %两球面镜的中心的距离a1=2e3*lamda; %左边腔镜的横截面尺寸a2=2e3*lamda; %右边腔镜的横截面尺寸
k=2*pi/lamda; %波矢大小
co1=1i*exp(-1i*k*L)/(L*lamda); %中间系数
N1=a1^2/(L*lamda); %菲涅尔数
N2=a2^2/(L*lamda);
B=zeros(E,F); %右边的初始光场分布
for e=1:E
for f=1:F
r2(e,f)=(e*a2)/E;
theta2(e,f)=(2*f*pi)/F;
x2(e,f)=r2(e,f)*cos(theta2(e,f));
y2(e,f)=r2(e,f)*sin(theta2(e,f));
for m=1:M
for n=1:N
r1(m,n)=(m*a1)/M;
theta1(m,n)=(2*n*pi)/N;
x1(m,n)=r1(m,n)*cos(theta1(m,n));
y1(m,n)=r1(m,n)*sin(theta1(m,n));
co(m,n)=exp(1i*k*(x2(e,f)*x1(m,n)+y2(e,f)*y1(m,n))); B(e,f)=B(e,f)+A(m,n)*(pi*a1^2/(M*N))*co(m,n)*co1;
end
end%solve B(e+(E+1)/2,f+(F+1)/2)
end
end
B=B/abs(max(max(B)));
A=zeros(M,N);
for m=1:M
for n=1:N
r1(m,n)=(m*a1)/M;
theta1(m,n)=(2*n*pi)/N;
x1(m,n)=r1(m,n)*cos(theta1(m,n));
y1(m,n)=r1(m,n)*sin(theta1(m,n));
for e=1:E
for f=1:F
r2(e,f)=(e*a2)/E;
theta2(e,f)=(2*f*pi)/F;
x2(e,f)=r2(e,f)*cos(theta2(e,f));
y2(e,f)=r2(e,f)*sin(theta2(e,f));
co(e,f)=exp(1i*k*(x1(m,n)*x2(e,f)+y1(m,n)*y2(e,f))); A(m,n)=A(m,n)+B(e,f)*(pi*a2^2/(E*F))*co(e,f)*co1;
end
end%solve B(e+(E+1)/2,f+(F+1)/2)
end
end
A=A/abs(max(max(A)));
r1=(1:1:M)*(a1/M);
theta1=2*pi/N:2*pi/N:2*pi;
r2=(1:1:E)*(a2/E);
theta2=2*pi/F:2*pi/F:2*pi;
[Rr1,THETA1]=meshgrid(r1,theta1);
[Rr2,THETA2]=meshgrid(r2,theta2);
Z1=abs(A)';
Z2=abs(B)';
[x1,y1,z1]=pol2cart(THETA1,Rr1,Z1);
[x2,y2,z2]=pol2cart(THETA2,Rr2,Z2); subplot(1,2,1);
mesh(x1,y1,z1);
subplot(1,2,2);
mesh(x2,y2,z2);
输入的指令和得到的结果如下：
clear all;
A=rand(48,48);
sphe_type
sphe_type
sphe_type
sphe_type
N1,N2
N1 =
4
N2 =
4
2.基于传输矩阵的仿真方法：bar_typem.m
1）条形腔的传输矩阵方法：
%对称条形谐振腔的自再现模
M=701; %条形腔左边的总的分割数目N=701; %条形腔右边的总的分割数目
lamda=6.328e-7;
L=1e5*lamda;
a=1e3*lamda;
b=1e3*lamda;
k=2*pi/lamda;
co=sqrt(1i*(exp(-1i*k*L))/(lamda*L)); N1=a^2/(L*lamda);
N2=b^2/(L*lamda);
x1=(-(M-1)/2:1:(M-1)/2)'*(2*a/(M-1)); x2=(-(N-1)/2:1:(N-1)/2)'*(2*b/(N-1));
x1_f=x1.^2;
x2_f=x2.^2;
x1_fm=repmat(x1_f,1,M);
x2_fm=repmat(x2_f,1,N);
det_xf=x1_fm+x2_fm'-2.*x1*x2';
T=co.*exp(-1i*k*det_xf/(2*L));
B=T*A;
B=B/abs(max(B));
A=T*B;
A=A/abs(max(A));
x1=(-(M-1)/2:1:(M-1)/2)*(2*a/(M-1)); x2=(-(N-1)/2:1:(N-1)/2)*(2*b/(N-1)); subplot(2,2,1);
plot(x1,abs(A));
subplot(2,2,3);
plot(x1,angle(A));
subplot(2,2,2);
plot(x2,abs(B));
subplot(2,2,4);
plot(x2,angle(B));
输入的指令和得到的结果：
clear all;
A=rand(701,1);
clear all;
A=rand(701,1);
tic;
for i=1:300
bar_typem;
end
toc
Elapsed time is 27.890094 seconds.
>> N1,N2
N1 =
9.999999999999998
N2 =
9.999999999999998
2）矩形平平腔的传输矩阵法仿真：pin_recm.m
源代码如下：
%传输矩阵法模拟对称矩形平平腔的自再现模
M=201; %左边镜面的横向离散点数
N=201; %左边镜面的纵向离散点数
E=201; %右边镜面的横向离散点数
F=201; %右边镜面的纵向离散点数
lamda=6.328e-7; %中心波长
L=1e6*lamda; %腔长
a=1e3*lamda; %左边镜面的横向尺寸
b=1e3*lamda; %左边镜面的纵向尺寸
c=1e3*lamda; %右边镜面的横向尺寸
d=1e3*lamda; %右边镜面的纵向尺寸
N1=(a*b)/(L*lamda); %左边镜面的菲涅尔数
N2=(c*d)/(L*lamda); %右边镜面的菲涅尔数
k=2*pi/lamda; %波数大小
co=sqrt(1i*(exp(-1i*k*L))/(lamda*L)); %中间系数
x1=(-(M-1)/2:1:(M-1)/2)'*(2*a/(M-1)); %左边镜面的离散点的横坐标
x2=(-(E-1)/2:1:(E-1)/2)'*(2*c/(E-1)); %右边镜面的离散点的横坐标
y1=(-(N-1)/2:1:(N-1)/2)'*(2*b/(N-1)); %左边镜面的离散点的纵坐标
y2=(-(F-1)/2:1:(F-1)/2)'*(2*d/(F-1)); %右边镜面的离散点的纵坐标
x1_f=x1.^2; %每个点坐标平方
x2_f=x2.^2;
y1_f=y1.^2;
y2_f=y2.^2;
x1_fm=repmat(x1_f,1,M); %左边镜面所有点的横坐标平方矩阵的转置x2_fm=repmat(x2_f,1,E); %右边镜面所有点的横坐标平方矩阵
y1_fm=repmat(y1_f,1,N); %左边镜面所有点的纵坐标平方矩阵的转置y2_fm=repmat(y2_f,1,F); %右边镜面所有点的横坐标平方矩阵det_xf=x1_fm+x2_fm'-2.*x1*x2'; %左右镜面对应点横坐标的差的平方det_yf=y1_fm+y2_fm'-2.*y1*y2'; %左右镜面对应点横坐标的差的平方TX=co.*exp(-1i*k*det_xf/(2*L)); %左右镜面横坐标传输矩阵
TY=co.*exp(-1i*k*det_yf/(2*L)); %左右镜面纵坐标传输矩阵
Bx=TX*Ax; %右边光场的横向分布
By=TY*Ay; %右边光场的纵向分布
Ax=TX*Bx; %左边光场的横向分布
Ay=TY*By; %右边光场的横向分布
A=Ay*Ax'; %左边光场的二维分布
A=A/abs(max(max(A))); %归一化
B=By*Bx'; %右边光场的二维分布
B=B/abs(max(max(B))); %归一化
[X1,Y1]=meshgrid(x1',y1'); %左边镜面二维坐标网格[X2,Y2]=meshgrid(x2',y2'); %右边镜面二维坐标网格subplot(2,2,1);
mesh(X1,Y1,abs(A)); %左边镜面光场的幅度分布subplot(2,2,3);
mesh(X1,Y1,angle(A)); %左边镜面光场的相位分布subplot(2,2,2);
mesh(X2,Y2,abs(B)); %右边镜面光场的幅度分布subplot(2,2,4);
mesh(X2,Y2,angle(B)); %右边镜面光场的相位分布
输入的指令和得到的结果如下：
clear all;
Ax=rand(201,1);
Ay=rand(201,1);
tic;
for i=1:10
pin_recm;
end
>> N1,N2
N1 =
1
N2 =
1
>> toc
Elapsed time is 21.508314 seconds.
3）矩形共焦腔的传输矩阵仿真：recfoc_matric.m
源代码如下：
%传输矩阵法模拟对称矩形平平腔的自再现模
N=201; %左边镜面的纵向离散点数
E=201; %右边镜面的横向离散点数
F=201; %右边镜面的纵向离散点数
lamda=6.328e-7; %中心波长
L=5e6*lamda; %腔长
a=1e3*lamda; %左边镜面的横向尺寸
b=1e3*lamda; %左边镜面的纵向尺寸
c=1e3*lamda; %右边镜面的横向尺寸
d=1e3*lamda; %右边镜面的纵向尺寸
N1=(a*b)/(L*lamda); %左边镜面的菲涅尔数
N2=(c*d)/(L*lamda); %右边镜面的菲涅尔数
k=2*pi/lamda; %波数大小
co=sqrt(1i*(exp(-1i*k*L))/(lamda*L)); %中间系数
x1=(-(M-1)/2:1:(M-1)/2)'*(2*a/(M-1)); %左边镜面的离散点的横坐标x2=(-(E-1)/2:1:(E-1)/2)'*(2*c/(E-1)); %右边镜面的离散点的横坐标y1=(-(N-1)/2:1:(N-1)/2)'*(2*b/(N-1)); %左边镜面的离散点的纵坐标y2=(-(F-1)/2:1:(F-1)/2)'*(2*d/(F-1)); %右边镜面的离散点的纵坐标TX=co.*exp(1i*k*x1*x2'/L); %左右镜面横坐标传输矩阵TY=co.*exp(1i*k*y1*y2'/L); %左右镜面纵坐标传输矩阵
Bx=TX*Ax; %右边光场的横向分布
By=TY*Ay; %右边光场的纵向分布
Ax=TX*Bx; %左边光场的横向分布
Ay=TY*By; %右边光场的横向分布
A=Ay*Ax'; %左边光场的二维分布
A=A/abs(max(max(A))); %归一化
B=By*Bx'; %右边光场的二维分布
B=B/abs(max(max(B))); %归一化
[X1,Y1]=meshgrid(x1',y1'); %左边镜面二维坐标网格
[X2,Y2]=meshgrid(x2',y2'); %右边镜面二维坐标网格subplot(2,2,1);
mesh(X1,Y1,abs(A)); %左边镜面光场的幅度分布subplot(2,2,3);
mesh(X1,Y1,angle(A)); %左边镜面光场的相位分布subplot(2,2,2);
mesh(X2,Y2,abs(B)); %右边镜面光场的幅度分布subplot(2,2,4);
mesh(X2,Y2,angle(B)); %右边镜面光场的相位分布
输入的指令和得到的结果如下：clear all;
Ax=ones(201,1);
Ay=ones(201,1);
for i=1:10
recfoc_matric
end
>> N1,N2
N1 =
0.2000
N2 =
0.2000。