MATLAB激光谐振腔仿真代码

引用激光光束传输与谐振腔仿真代码包ABCDRez的内容

文档阅读：

激光光束传输与谐振腔仿真ABCDRez_abcdrez csdn-CSDN博客

文件下载：

Laser Beam Propagation and Resonator Simulation ABCDRez - File Exchange - MATLAB Central

https://download.csdn.net/download/qq_42712244/89533054

3.2 驻波谐振腔

可使用矩阵光学方法对多元件的稳定腔进行分析，采用“G参数等价腔”来描述腔内束宽。参见吕百达教授著《激光光学》第十章第一节。

3.2.1 两个厚透镜晶体加曲面腔镜的谐振腔

使用“G参数等价腔”相关内容，编写具体函数Rez4mThick。

![](https://i-blog.csdnimg.cn/direct/398e0740c4254ae292df89e28d0d9355.png)

图3.2.1.1 谐振腔示意图

图3.2.1.2 腔内束宽示意图

[www,wthetaL0,FlagRez]=Rez4mThick(lambda,RezPara,str)

变量

意义

lambda

波长

RezPara

谐振腔参数

str

若str为’plot’则绘图；否则，不绘制

www

有效的束宽（这里用来描述模式体积）

wthetaL0

靠近输出镜的光束的束腰、发散角、光腰位置信息

FlagRez

谐振腔的标志信息

其中：有效束宽www为振荡光斑在晶体四个端面的平均值。

略去绘图部分的代码，重新改写函数为

[www,wthetaL0,FlagRez]=Rez4mThick00(lambda,RezPara)

变量

意义

lambda

波长

RezPara

谐振腔参数

www

有效的束宽（这里用来描述模式体积）

wthetaL0

靠近输出镜的光束的束腰、发散角、光腰位置信息

FlagRez

谐振腔的标志信息

这样，我们就可以快速计算大量腔型，并根据结果进行筛选！

如查看某一腔型的有效束宽随热焦距的变化情况

next
%% 腔参数
lambda=1064*10^-9;
%
rho1=-1500*10^-3;
d1=(10+0)*10^-3;
% F1=500*10^-3;
lenF1=30*10^-3;nF1=1.8;
d3=15*10^-3;
% F2=F1;
lenF2=lenF1;nF2=nF1;
d2=(40+0)*10^-3;
rho2=-2500*10^-3;
%
limR=5e-3;
LL=d1+lenF1+d3+lenF2+d2;
%% 暴力求解，耗时短
F1x=linspace(250,2500,250).*10^-3;
www=zeros(size(F1x));
for ii=1:length(F1x)
    F1=F1x(ii);
    F2=F1;
    RezPara=[rho1,d1,F1,lenF1,nF1,d3,F2,lenF2,nF2,d2,rho2];
    [www(ii),~,~]=Rez4mThick00(lambda,RezPara);
end
% 剔除超过边界的数据
www( (www>limR) )=0;
%% 绘图
plot(F1x,www.*10^3);
grid on;hold on;
xlabel('焦距/m');
ylabel('有效束宽/mm')
title('有效束宽随热焦距变化情况')
%% 版本信息
% 作者:                Quincy Howard
% 联系方式：           quincy.hd@qq.com
% 文件信息更新平台为   matlab云文件、GitHub、CSDN博客、知乎
% 若使用请注明来源
% 最后编辑于           2024 年 07 月 10 日

如查看某一腔型的有效束宽随两腔镜曲率的变化情况

如查看某一腔型的稳定性随两腔镜曲率的变化情况

next
%%
lambda=1064e-9;
%
% rho1=-1500*10^-3;
d1=(10+0)*10^-3;
F1=500*10^-3;lenF1=30*10^-3;nF1=1.8;
d3=15*10^-3;
F2=F1;lenF2=lenF1;nF2=nF1;
d2=(40+0)*10^-3;
% rho2=-2500*10^-3;
%
limR=5e-3;
LL=d1+lenF1+d3+lenF2+d2;
%% 暴力求解，耗时较短
tic
Rho=[5000,4500,4000,3500,3000,2500,2000,1500,1250,1000,800,600,500,450,400,350,300,250]*10^-3;
Rhoall=[inf,Rho,flip(-Rho),-inf];
www=zeros(length(Rhoall),length(Rhoall));
FlagRez=zeros(length(Rhoall),length(Rhoall));
for ii=1:length(Rhoall)
    rho1=Rhoall(ii);
    for jj=1:length(Rhoall)
        rho2=Rhoall(jj);
        RezPara=[rho1,d1,F1,lenF1,nF1,d3,F2,lenF2,nF2,d2,rho2];
        [www(ii,jj),~,FlagRez(ii,jj)]=Rez4mThick00(lambda,RezPara);
    end
end
toc
% 剔除超过边界的数据
www( (www>limR) )=0;
%% 绘图
[XX,YY]=meshgrid(Rhoall,Rhoall);
figure;
mesh(XX.*10^3,YY.*10^3,www.*10^3);
xlabel('rho1/mm');ylabel('rho2/mm');zlabel('www/mm');
title('有效束宽随两腔镜曲率的变化情况');
figure;
mesh(XX.*10^3,YY.*10^3,FlagRez);
xlabel('rho1/mm');ylabel('rho2/mm');zlabel('FlagRez');
title('稳定性随两腔镜曲率的变化情况');
%% 版本信息
% 作者:                Quincy Howard
% 联系方式：           quincy.hd@qq.com
% 文件信息更新平台为   matlab云文件、GitHub、CSDN博客、知乎
% 若使用请注明来源
% 最后编辑于           2024 年 07 月 10 日

有效束宽随两腔镜曲率的变化情况:

稳定性随两腔镜曲率的变化情况:

当然，也可调用ABCDRez中其他谐振腔模拟函数查看相关变化情况，这里就不在赘述。

其中引用的函数如下：


function [www,wthetaL0,FlagRez]=Rez4mThick00(lambda,RezPara)
% RezPara=[rho1,d1,F1,lenF1,nF1,d3,F2,lenF2,nF2,d2,rho2];
rho1=    RezPara(01);d1=    RezPara(02);F1=     RezPara(03);
lenF1=   RezPara(04);nF1=   RezPara(05);d3=     RezPara(06);
F2=      RezPara(07);lenF2= RezPara(08);nF2=    RezPara(09);
d2=      RezPara(10);rho2=  RezPara(11);
%% 不绘图，用于大量计算
syms z
%总腔长
LL=d1+lenF1+d3+lenF2+d2;
[~,MF1,~,~]=deF2Rho(F1,lenF1,nF1);
if (F1==F2 && lenF1==lenF2 && nF1==nF2)
    MF2=MF1;
else
    [~,MF2,~,~]=deF2Rho(F2,lenF2,nF2);
end
%
MM=[1,d2;0,1]*MF2*[1,d3;0,1]*MF1*[1,d1;0,1];
a=MM(1,1);b=MM(1,2);c=MM(2,1);d=MM(2,2);
%令
G1=a-b/rho1;
G2=d-b/rho2;
FlagRez(1)=G1*G2;
%% 是否满足0<G1*G2<1稳定条件
% figure(22);
if ( FlagRez(1)<0  ||  FlagRez(1)>1 )
    www=0;
    wthetaL0=zeros(1,3);
else
    %成立稳定性条件
    %0<G1*G2<1
    %% 计算
    %靠近两腔镜束腰大小
    w001=abs(sqrt(+(lambda*b/pi)*(sqrt(G1*G2*(1-G1*G2))/abs(G1+a*a*G2-2*a*G1*G2))));
    w002=abs(sqrt(+(lambda*b/pi)*(sqrt(G1*G2*(1-G1*G2))/abs(G2+d*d*G1-2*d*G1*G2))));
    %束腰位置
    L001=(b*G2*(a-G1))/(G1+a*a*G2-2*a*G1*G2);
    L002=LL-(  (b*G1*(d-G2))/(G2+d*d*G1-2*d*G1*G2)  );
    %基模发散角
    theta001=sqrt((lambda/(pi*w001))^2);
    theta002=sqrt((lambda/(pi*w002))^2);
    %% 腔内束宽
    %腔内的束宽变换
    %由镜1全反镜向输出镜变换
    [w003,theta003,L003]=fLRMm(w001,theta001,L001,lambda,MF1,lenF1,d1);
    %
    w001_z=fwz(w001,theta001,L001);
    w003_z=fwz(w003,theta003,L003);
    w002_z=fwz(w002,theta002,L002);
    www=1/4*double(subs(w001_z,z,d1)+subs(w003_z,z,d1+lenF1)+...
        subs(w003_z,z,d1+lenF1+d3)+subs(w002_z,z,d1+lenF1+d3+lenF2));
    wthetaL0=double([w002,theta002,L002]);
   
end
end
%% 版本信息
% 作者:                Quincy Howard
% 联系方式：           quincy.hd@qq.com
% 文件信息更新平台为   matlab云文件、GitHub、CSDN博客、知乎
% 若使用请注明来源
% 最后编辑于           2024 年 07 月 10 日


function [Rho,MF,MFront,MBack]=deF2Rho(F,lenF,nF)
    % F 厚透镜有效焦距
    % lenF 厚透镜长度
    % nF 厚透镜介质折射率
    n1=1;
    % 曲率表示
    syms phox
    MFront=[1,0;(nF-n1)/(nF*-phox),n1/nF];
    Mmedian=[1,lenF;0,1];
    MBack=[1,0;(n1-nF)/(n1*phox),nF/n1];
    Mf1=MBack*Mmedian*MFront;
    %求解代换
    phoxx=vpasolve(Mf1(2,1)==-1/F);
    [~,nn]=max(abs(phoxx));
    Rho=phoxx(nn);
    % % % % % % % %
    % 多个解？？！！
    % % % % % % % %
    
    MFront=double([1,0;(nF-n1)/(nF*-Rho),n1/nF]);
    Mmedian=double([1,lenF;0,1]);
    MBack=double([1,0;(n1-nF)/(n1*Rho),nF/n1]);
    MF=MBack*Mmedian*MFront;
end
%% 版本信息
% 作者:                Quincy Howard
% 联系方式：           quincy.hd@qq.com
% 文件信息更新平台为   matlab云文件、GitHub、CSDN博客、知乎
% 若使用请注明来源
% 最后编辑于           2024 年 07 月 10 日

%% 高阶高斯光束由左向右传输
%% [w02m,theta02m,L02m]=fLRMm(w01m,theta01m,L01m,lambda,cM,lencM,dcM)
function [w02m,theta02m,L02m]=fLRMm(w01m,theta01m,L01m,lambda,cM,lencM,dcM)
%         高阶转基模
        [w01,theta01,L01,M2]=fM2B(w01m,theta01m,L01m,lambda);
%         基模由左向右传输
        [w02,theta02,L02]=fLRM(w01,theta01,L01,lambda,cM,lencM,dcM);
%         基模转高阶
        [w02m,theta02m,L02m,~]=fB2M(w02,theta02,L02,M2);
end
%% 版本信息
% 作者:                Quincy Howard
% 联系方式：           quincy.hd@qq.com
% 文件信息更新平台为   matlab云文件、GitHub、CSDN博客、知乎
% 若使用请注明来源
% 最后编辑于           2024 年 07 月 10 日

%% 高阶模转基模
%% [w01,theta01,L01,M2]=fM2B(w01m,theta01m,L01m,lambda)
function [w01,theta01,L01,M2]=fM2B(w01m,theta01m,L01m,lambda)
    % w01m为高阶光束光腰 半径
    % theta01m为高阶光束  半角发散角
    % L01m为高阶光束光腰位置
    % M2为光束质量
    % w01为基模输出光腰 半径
    % theta01为基模光束  半角发散角
    % L01为基模光束光腰位置
    
    %默认使用1064nm
    % lambda=1064*10^-9;
    %默认在空气中折射率nn=1
    nn=1;
    
    M2=(w01m*theta01m*pi*nn)/(lambda);
    w01=(w01m)/sqrt(M2);
    theta01=(theta01m)/sqrt(M2);
    L01=L01m;
    
end
%% 版本信息
% 作者:                Quincy Howard
% 联系方式：           quincy.hd@qq.com
% 文件信息更新平台为   matlab云文件、GitHub、CSDN博客、知乎
% 若使用请注明来源
% 最后编辑于           2024 年 07 月 10 日


%% 基模由左向右传输
%% [w02,theta02,L02]=fLRM(w01,theta01,L01,lambda,cM,lencM,dcM)
function [w02,theta02,L02]=fLRM(w01,theta01,L01,lambda,cM,lencM,dcM)
% 默认使用基模计算,由左向右传输
% w01,theta01,L01 为入射光束参数
% lambda 为波长
% cM 为光学系统的传输矩阵
% lencM 为光学系统长度
% dcM 为 入射光束与光学系统的第一接触面 相对于 原点所处位置
% w02,theta02,L02 为出射光束参数
%% 物方光束
Z0=w01/theta01;
syms z;
w_1_z=w01*sqrt(1+((z-L01)/Z0)^2);
R_z1=Z0*((z-L01)/Z0+Z0/(z-L01));
R01=subs(R_z1,z,eps);
w_001=subs(w_1_z,z,0);
%% 传输矩阵相关定义及计算
syms X1 Y1 z
% syms R_1
% q参数定义
% syms s1 s2 a b c d
% 取n1=n2=1;
n1=1;n2=1;
%传输矩阵元素
s2=z;
a=cM(1,1);b=cM(1,2);c=cM(2,1);d=cM(2,2);
s1=dcM;
%其他特征参量
M=[1 s2;0,1]*[a b;c d]*[1,s1;0,1];
A=M(1,1);B=M(1,2);
% C=M(2,1);D=M(2,2);
Y2=(n1/n2)*Y1/...
    (A^2+2*X1*A*B+(X1^2+Y1^2)*B^2);
%% 像方光束
w_2=sqrt(lambda/(pi*Y2));
w_2=subs(w_2,[X1 Y1],[1/R01 lambda/(pi*w_001^2)]);
%光腰参考点平移
w_2=subs(w_2,(z),(z-s1));
%% 输出参数
% disp('w_2*w_2=');pretty(w_2*w_2);
T1=vpa((taylor((w_2*w_2),z,'Order',1)));
T2=vpa((taylor((w_2*w_2),z,'Order',2)));
T3=vpa((taylor((w_2*w_2),z,'Order',3)));
C=T1;
B=(T2-T1)/z;
A=(T3-T2)/(z*z);
w02=vpa(sqrt(C-B*B/(4*A)));
theta02=vpa(sqrt(A));
L02=vpa(-B/(2*A));
L02=L02+lencM;
end
%% 版本信息
% 作者:                Quincy Howard
% 联系方式：           quincy.hd@qq.com
% 文件信息更新平台为   matlab云文件、GitHub、CSDN博客、知乎
% 若使用请注明来源
% 最后编辑于           2024 年 07 月 10 日

%% 基模转高阶模
%% [w01m,theta01m,L01m,M2]=fB2M(w01,theta01,L01,M2)
function [w01m,theta01m,L01m,M2]=fB2M(w01,theta01,L01,M2)
    % w01m为高阶光束光腰 半径
    % theta01m为高阶光束  半角发散角
    % L01m为高阶光束光腰位置
    % M2为光束质量
    % w01为基模输出光腰 半径
    % theta01为基模光束  半角发散角
    % L01为基模光束光腰位置
    %     M2=(w01m*theta01m*pi*nn)/(lambda)
    % 默认在空气中折射率nn=1
    %     nn=1;
    %     lambda=(w01m*theta01m*pi*nn)/(M2);
    w01m=(w01)*sqrt(M2);
    theta01m=(theta01)*sqrt(M2);
    L01m=L01;
    
end
%% 版本信息
% 作者:                Quincy Howard
% 联系方式：           quincy.hd@qq.com
% 文件信息更新平台为   matlab云文件、GitHub、CSDN博客、知乎
% 若使用请注明来源
% 最后编辑于           2024 年 07 月 10 日

%% 光束随距离变化
%% w_z = fwz(w0,theta0,L0,z)
function w_z = fwz(w0,theta0,L0,z)
    if nargin==3
        syms z;
    end
    
    w_z=sqrt(w0^2+theta0^2*(z-L0)^2);
end
%% 版本信息
% 作者:                Quincy Howard
% 联系方式：           quincy.hd@qq.com
% 文件信息更新平台为   matlab云文件、GitHub、CSDN博客、知乎
% 若使用请注明来源
% 最后编辑于           2024 年 07 月 10 日

** next脚本:**

close all;
clear;
clc;
%% 版本信息
% 作者:                Quincy Howard
% 联系方式：           quincy.hd@qq.com
% 文件信息更新平台为   matlab云文件、GitHub、CSDN博客、知乎
% 若使用请注明来源
% 最后编辑于           2024 年 07 月 10 日

标签：算法 matlab github

本文转载自: https://blog.csdn.net/qq_42712244/article/details/140730397
版权归原作者 QuincyHoward 所有，如有侵权，请联系我们删除。

MATLAB激光谐振腔仿真代码

引用激光光束传输与谐振腔仿真代码包ABCDRez的内容

文档阅读：

文件下载：

3.2 驻波谐振腔

3.2.1 两个厚透镜晶体加曲面腔镜的谐振腔

如查看某一腔型的有效束宽随热焦距的变化情况

如查看某一腔型的有效束宽随两腔镜曲率的变化情况

如查看某一腔型的稳定性随两腔镜曲率的变化情况

发表评论

“MATLAB激光谐振腔仿真代码”的评论:

关于作者

overfit同步小助手

相关阅读

文章导航