?? 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 + -