concavefalt5.m

平凹腔激光模式计算程序

clear  %方形平凹镜腔边界元迭代法高斯分布反射率双镜倾斜计算（平面镜光场分布、远场分布、本征值）
for t1=0.01:0.01:0.13
    for t2=0.01:0.01:0.13
L=340;  %腔长
a=1;  %镜宽
R1=500;    %凹镜曲率半径
W=10.6e-4;  %波长
N=a^2/(W*L)  %菲涅尔数
M=100;  %分元数
t=0.247   %倾斜角度
q1=t1*pi/180;  %倾斜角度的弧度数
q2=t2*pi/180;  %倾斜角度的弧度数
K=2*pi/W;   %K是波数。
G1X=zeros(M);  %X方向凹镜向平镜度越矩阵
G2X=zeros(M);    %X方向平镜向凹镜度越矩阵
G1Y=zeros(M);   %Y方向凹镜向平镜度越矩阵
G2Y=zeros(M);   %Y方向平镜向凹镜度越矩阵
for m=1:M
   for n=1:M
      G1X(n,m)=(2*a/M)*sqrt((i/(W*L))*exp(-i*K*L))*exp(-i*K*(((-a+2*a*n/M)-(-a+2*a*m/M))^2/(2*L)-(-a+2*a*m/M)^2/(2*R1)+(-a+2*a*n/M)*sin(q1)));
      G2X(m,n)=f1*exp(-(-a+2*a*n/M)^2/w1^2)*(2*a/M)*sqrt((i/(W*L))*exp(-i*K*L))*exp(-i*K*(((-a+2*a*m/M)-(-a+2*a*n/M))^2/(2*L)-(-a+2*a*m/M)^2/(2*R1)+(-a+2*a*n/M)*sin(q1)));
      G1Y(n,m)=(2*a/M)*sqrt((i/(W*L))*exp(-i*K*L))*exp(-i*K*(((-a+2*a*n/M)-(-a+2*a*m/M))^2/(2*L)-(-a+2*a*m/M)^2/(2*R1)+(-a+2*a*n/M)*sin(q2)));
      G2Y(m,n)=f1*exp(-(-a+2*a*n/M)^2/w1^2)*(2*a/M)*sqrt((i/(W*L))*exp(-i*K*L))*exp(-i*K*(((-a+2*a*m/M)-(-a+2*a*n/M))^2/(2*L)-(-a+2*a*m/M)^2/(2*R1)+(-a+2*a*n/M)*sin(q2)));
   end
end
GX=G1X*G2X;%算出为平镜光场
GY=G1Y*G2Y;
%save 0000 G1X G2X G1Y G2Y
clear G1X G1Y G2X G2Y
%以上为总矩阵。
%save 0001 GX GY
E1X=ones(M,1);  %平镜X方向光场分布
E1Y=ones(M,1);  %平镜Y方向光场分布
for d=1:5000
    E1X=GX*E1X;
    E1X=E1X/max(abs(E1X));
    E1Y=GY*E1Y;
    E1Y=E1Y/max(abs(E1Y));
end
MX=max(abs(GX*E1X));  %X向本征值
MY=max(abs(GY*E1Y));  %Y向本征值
M1=MX*MY;  %总本征值
E=E1X*E1Y.';   %平镜总光场分布
[maxcol,maxindexX]=max(E);  %极大值x坐标
[maxrow,maxindexY]=max(maxcol)  %极大值y坐标
for m2=1:M
indexX(m2)=-a+2*a*m2/M;
indexY(m2)=indexX(m2);
end
mesh(indexX,indexY,abs(E));  %光场幅度画图
pause   %程序暂停*******************************************
mesh(indexX,indexY,angle(E));  %画图
pause  %程序暂停*******************************************
for m2=1:M
    for n2=1:M
        QXG(m2,n2)=(2*a/M)*exp(-i*K*((-15*W/(2*a)+15*2*W*m2/(2*a*M))*(-a+2*a*n2/M)))*sqrt(exp(-i*K*D));
        QYG(m2,n2)=(2*a/M)*exp(-i*K*((-15*W/(2*a)+15*2*W*m2/(2*a*M))*(-a+2*a*n2/M)))*sqrt(exp(-i*K*D));
    end
end
QX=QXG*E1X;  %远场X向分布
QY=QYG*E1Y;  %远场Y向分布
Q=QX*QY.';   %远场总分布
indexx=zeros(1,M);   %远场分布坐标
 indexy=zeros(1,M);   %远场分布坐标
for m1=1:M
    indexx(m1)=-15+15*m2/M;
    indexy(m1)=indexx;
end
mesh(indexx,indexy,abs(Q));  %远场强度画图
end
end