📄 tmmcore.m
字号:
% 本程序用TMM求解Bragg光纤
clear all;
format long
tic;
% Constant Parameter
% 参数选择见文献“S. Guo, S. Albin, R. Rogowski,"Comparative analysis of Bragg fibers" Opt. Exp. 12 (2004) 198.”
L=15; % 包层周期数
N=2*L+1; % 边界数
epsilon0=8.854187818e-12; % 真空的介电常数
mu0=4*pi*10^(-7); % 真空的磁导率
Z=120*pi; % 自由空间波阻抗
period=2e-6; % 包层的周期
ncore=1; % 纤芯折射率
nhigh=2; % 包层中高折射率介质的折射率
nlow=1; % 包层中低折射率介质的折射率
nextern=1; % 包层外的折射率
rcore=5e-6; % 纤芯半径
rhigh=1e-6; % 包层中高折射率介质的厚度
rlow=1e-6; % 包层中低折射率介质的厚度
thick=(rhigh+rlow)*L+rcore; % 整个bagg光纤的厚度
flag=1; % flag=1时包层中第一层是高折射率层；flag=0时包层中第一层时低折射率层
m=0; % m=0为TE或者TM模，m>1为混合模
M=50; % 求解前M个模式的传播常数和场分布
interface=zeros(N,1);
interface(1)=rcore;
if(flag==1) % 第一层是高折射率层
for i=2:N
interface(i)=interface(1)+i/2*rhigh+(i-i/2-1)*rlow;
end
elseif(flag==0) % 第一层是低折射率层
for i=2:N
interface(i)=interface(1)+i/2*rlow+(i-i/2-1)*rhigh;
end
end
U=100;
lambda=linspace(1e-6,6e-6,U); % 工作波长
c=299792458; % 真空中的光速
omega=2.*pi.*c./lambda; % 角频率
k=2.*pi./lambda; % 空间波数
V=20;
neff=linspace(0+1e-4,1-1e-4,V); % 不能取1是为了避免以后除零
row1=[1;0;1;0];
% T(i)为第i层的传输矩阵
% TT为总传输矩阵
for u=1:U
jj=1;
kk=1;
beta=k(u).*neff;
kcore=k(u).*sqrt(ncore.^2-neff.^2); % 纤芯中的横向波数
khigh=k(u).*sqrt(nhigh.^2-neff.^2); % 高折射率层中的横向波数
klow=k(u).*sqrt(nlow.^2-neff.^2); % 低折射率层中的横向波数
kextern=k(u).*sqrt(nextern.^2-neff.^2); % 包层外的横向波数
% 纤芯中的贝塞尔函数和系数矩阵
rrcore=kcore.*rcore;
Jm_core=besselj(m,rrcore);
Ym_core=bessely(m,rrcore);
Jmdiff_core=-besselj(m+1,rrcore)+m./rrcore.*besselj(m,rrcore(1:V));
Ymdiff_core=-bessely(m+1,rrcore)+m./rrcore.*bessely(m,rrcore(1:V));
% 最外层的贝塞尔函数
rrextern=kextern.*interface(N);
Jm_extern=besselj(m,rrextern);
Ym_extern=bessely(m,rrextern);
Jmdiff_extern=-besselj(m+1,rrextern)+m./rrextern.*besselj(m,rrextern);
Ymdiff_extern=-bessely(m+1,rrextern)+m./rrextern.*bessely(m,rrextern);
for v=1:V
% 纤芯的系数矩阵
Mcore=[Jm_core(v),Ym_core(v),0,0;
i*omega(u)*epsilon0*ncore^2*Jmdiff_core(v)/kcore(v),i*omega(u)*epsilon0*ncore^2*Ymdiff_core(v)/kcore(v),-m*beta(v)*Jm_core(v)/(kcore(v)^2*rcore),-m*beta(v)*Ym_core(v)/(kcore(v)^2*rcore);
0,0,Jm_core(v),Ym_core(v);
-m*beta(v)*Jm_core(v)/(kcore(v)^2*rcore),-m*beta(v)*Ym_core(v)/(kcore(v)^2*rcore),i*omega(u)*mu0*Jmdiff_core(v)/kcore(v),i*omega(u)*mu0*Ymdiff_core(v)/kcore(v)];
% 最外层的系数矩阵
Mextern=[Jm_extern(v),Ym_extern(v),0,0;
i*omega(u)*epsilon0*nextern^2*Jmdiff_extern(v)/kextern(v),i*omega(u)*epsilon0*nextern^2*Ymdiff_extern(v)/kextern(v),-m*beta(v)*Jm_extern(v)/(kextern(v)^2*interface(N)),-m*beta(v)*Ym_extern(v)/(kextern(v)^2*interface(N));
0,0,Jm_extern(v),Ym_extern(v);
-m*beta(v)*Jm_extern(v)/(kextern(v)^2*interface(N)),-m*beta(v)*Ym_extern(v)/(kextern(v)^2*interface(N)),i*omega(u)*mu0*Jmdiff_extern(v)/kextern(v),i*omega(u)*mu0*Ymdiff_extern(v)/kextern(v)];
TT=Mcore;
for j=2:N
if(mod(j,2)==0) % 高折射率层中的贝塞尔函数
% 高折射率层左边的贝塞尔函数
rrhigh1=khigh.*interface(j-1);
Jm_high1=besselj(m,rrhigh1);
Ym_high1=bessely(m,rrhigh1);
Jmdiff_high1=-besselj(m+1,rrhigh1)+m./rrhigh1.*besselj(m,rrhigh1);
Ymdiff_high1=-bessely(m+1,rrhigh1)+m./rrhigh1.*bessely(m,rrhigh1);
% 高折射率层左边的系数矩阵
Mhigh1=[Jm_high1(v),Ym_high1(v),0,0;
i*omega(u)*epsilon0*nhigh^2*Jmdiff_high1(v)/khigh(v),i*omega(u)*epsilon0*nhigh^2*Ymdiff_high1(v)/khigh(v),-m*beta(v)*Jm_high1(v)/(khigh(v)^2*rhigh),-m*beta(v)*Ym_high1(v)/(khigh(v)^2*rhigh);
0,0,Jm_high1(v),Ym_high1(v);
-m*beta(v)*Jm_high1(v)/(khigh(v)^2*rhigh),-m*beta(v)*Ym_high1(v)/(khigh(v)^2*rhigh),i*omega(u)*mu0*Jmdiff_high1(v)/khigh(v),i*omega(u)*mu0*Ymdiff_high1(v)/khigh(v)];
% 高折射率层右边的贝塞尔函数
rrhigh2=khigh.*interface(j);
Jm_high2=besselj(m,rrhigh2);
Ym_high2=bessely(m,rrhigh2);
Jmdiff_high2=-besselj(m+1,rrhigh2)+m./rrhigh1.*besselj(m,rrhigh2);
Ymdiff_high2=-bessely(m+1,rrhigh2)+m./rrhigh1.*bessely(m,rrhigh2);
% 高折射率层右边的系数矩阵
Mhigh2=[Jm_high2(v),Ym_high2(v),0,0;
i*omega(u)*epsilon0*nhigh^2*Jmdiff_high2(v)/khigh(v),i*omega(u)*epsilon0*nhigh^2*Ymdiff_high2(v)/khigh(v),-m*beta(v)*Jm_high2(v)/(khigh(v)^2*rhigh),-m*beta(v)*Ym_high2(v)/(khigh(v)^2*rhigh);
0,0,Jm_high2(v),Ym_high2(v);
-m*beta(v)*Jm_high2(v)/(khigh(v)^2*rhigh),-m*beta(v)*Ym_high2(v)/(khigh(v)^2*rhigh),i*omega(u)*mu0*Jmdiff_high2(v)/khigh(v),i*omega(u)*mu0*Ymdiff_high2(v)/khigh(v)];
T=Mhigh1*inv(Mhigh2);
TT=TT*T;
else % 低折射率层中的贝塞尔函数
% 低折射率层左边的贝塞尔函数
rrlow1=klow.*interface(j-1);
Jm_low1=besselj(m,rrlow1);
Ym_low1=bessely(m,rrlow1);
Jmdiff_low1=-besselj(m+1,rrlow1)+m./rrlow1.*besselj(m,rrlow1);
Ymdiff_low1=-bessely(m+1,rrlow1)+m./rrlow1.*bessely(m,rrlow1);
% 低折射率层左边的系数矩阵
Mlow1=[Jm_low1(v),Ym_low1(v),0,0;
i*omega(u)*epsilon0*nlow^2*Jmdiff_low1(v)/klow(v),i*omega(u)*epsilon0*nlow^2*Ymdiff_low1(v)/klow(v),-m*beta(v)*Jm_low1(v)/(klow(v)^2*rlow),-m*beta(v)*Ym_low1(v)/(klow(v)^2*rlow);
0,0,Jm_low1(v),Ym_low1(v);
-m*beta(v)*Jm_low1(v)/(klow(v)^2*rlow),-m*beta(v)*Ym_low1(v)/(klow(v)^2*rlow),i*omega(u)*mu0*Jmdiff_low1(v)/klow(v),i*omega(u)*mu0*Ymdiff_low1(v)/klow(v)];
% 低折射率层右边的贝塞尔函数
rrlow2=klow.*interface(j);
Jm_low2=besselj(m,rrlow2);
Ym_low2=bessely(m,rrlow2);
Jmdiff_low2=-besselj(m+1,rrlow2)+m./rrlow2.*besselj(m,rrlow2);
Ymdiff_low2=-bessely(m+1,rrlow2)+m./rrlow2.*bessely(m,rrlow2);
% 低折射率层右边的系数矩阵
Mlow2=[Jm_low2(v),Ym_low2(v),0,0;
i*omega(u)*epsilon0*nlow^2*Jmdiff_low2(v)/klow(v),i*omega(u)*epsilon0*nlow^2*Ymdiff_low2(v)/klow(v),-m*beta(v)*Jm_low2(v)/(klow(v)^2*rlow),-m*beta(v)*Ym_low2(v)/(klow(v)^2*rlow);
0,0,Jm_low2(v),Ym_low2(v);
-m*beta(v)*Jm_low2(v)/(klow(v)^2*rlow),-m*beta(v)*Ym_low2(v)/(klow(v)^2*rlow),i*omega(u)*mu0*Jmdiff_low2(v)/klow(v),i*omega(u)*mu0*Ymdiff_low2(v)/klow(v)];
T=Mlow1*inv(Mlow2);
TT=TT*T;
end
end
TT=TT*Mextern;
rowN=TT*row1;
if(sqrt(rowN(1)^2+rowN(2)^2)>100) % TM模带隙
nTM(jj)=neff(v);
jj=jj+1;
elseif(sqrt(rowN(1)^2+rowN(2)^2)>10^5) % TE模带隙
nTE(kk)=neff(v);
kk=kk+1;
end
end
% plot(lambda(u),n,'r')
% hold on
% plot(lambda(u),nTE,'r')
% hold on
% plot(lambda(u),nTM,'b')
% hold on
end
toc⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -