fiber_pml.m

fiber trasmission for PML boundary

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                       1D标量4th有限差分算法求解光纤                                            %
% 
% 创建日期：    2006年8月21日                                                                                    %
% 最近修改时间：2006年7月28日                                                                                    %
% 创建人：      laoer                                                                                            %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
clc;
format long
% Constant Parameter
% 参数选择见文献“Numerical analysis of bragg fiber using a compact 1D finite-difference frequency-domain menthod”
tic;
Z=120*pi;                            %  自由空间波阻抗
lambda=1.5e-6;                       %  工作波长
k0=2*pi/lambda;                      %  空间波数     
ncore=1.45;                          %  纤芯折射率
nclad=1;                             %  包层的折射率
rcore=3e-6;                          %  纤芯半径
rclad=3e-6;                          %  包层的厚度
thick=rcore+rclad;                   %  整个光纤的厚度
m=1;                                 %  m=0为TE或者TM模，m>0为混合模
t=100;                               %  纤芯内的离散点数为t+1
h=rcore/t;                           %  步长,这样选择步长使得纤芯和包层交界面上有离散点
T=t+round(rclad/h);                  %  半径中沿径向的总离散点数为T+1                          
M=1;                                 %  求解前M个模式的传播常数和场分布
L_PML=10;                            %  PML中的点数   
N=T+L_PML;                           %  计算区域的总点数
h_PML=h;                             %  PML中的步长
thick_PML=L_PML*h_PML;               %  PML的厚度
Thick=thick+thick_PML;               %  计算区域的总厚度
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 每点的相对介电常数和坐标
[epsilon,r1,r2]=e_PML(ncore,nclad,rcore,rclad,thick,T,h,h_PML,N,t)
% 将相对介电常数向量变为对角阵
epsilon=spdiags(epsilon,0,N,N);
% r1是每个离散点的坐标，r2是每个离散点前面半个步长的坐标
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% PML边界条件,这里采用复坐标伸缩PML
% PML参数
n=3;
alfi=200;        % 衰减系数
clear i;
for j=T+1:N
    r1(j)=r1(j)+i*alfi*(r1(j)-thick)^(n+1)/((n+1)*thick_PML^n);
    r2(j)=r2(j)+i*alfi*(r2(j)-thick)^(n+1)/((n+1)*thick_PML^n);
end
h1=zeros(N,1);   % PML中T+0.5...N-0.5 的步长r1(i)-r1(i-1)
h2=zeros(N,1);   % PML中T+1...N 的步长r2(i+1)-r2(i)
h1(1:T)=h;
h2(1:T-1)=h;
for i=T+1:N
    h1(i)=r1(i)-r1(i-1);
end
for i=T:N-1
    h2(i)=r2(i+1)-r2(i);
end 
clear i;  % 避免下面的i用到上面的值
% h2(N)要用到PML外面半个步长的点
h2(N)=Thick+h_PML/2+i*alfi*(Thick+h_PML/2-thick)^(n+1)/((n+1)*thick_PML^n)-r2(N);
%                               |
%                               v
%                         相当于r2(N+1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 矩阵及本征值
A=matrix_PML(rcore,ncore,nclad,T,h,N,h1,h2,r1,r2,t,k0,epsilon)
guess=(k0*ncore)^2;
options.tol = 1e-8;
options.disp = 0;
[Vector,D]=eigs(A,M,guess,options);  % D1为特征值，V的列向量为对应特征向量。
j=1;
for i=1:M
    if(real(sqrt(D(i,i))/k0)<ncore&real(sqrt(D(i,i))/k0)>nclad)
        neff(j)=sqrt(D(i,i))/k0;
        vector(:,j)=Vector(:,i);    % 取出对应与某个本征值的特征向量
        j=j+1;
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 场分布
S=5;    % 每隔S个点画一次
n=size(1:S:T,2);
r3=zeros(T,1);
for i=1:T
    r3(i)=r1(i,i); 
end
rr(1:n)=r3(1:S:T);  % 横坐标
r4=spdiags(r1);     % 取出对角元素
r5=zeros(T,1);
r5(1:T)=r4(1:T);


% 电场分布
col1=size(vector1,2);  % 得到列数
% Er(T,N,S,n,rr,r5,Z,vector1,col1)
Efi(T,N,S,n,rr,r5,Z,vector1,col1)
% Ez(k0,T,N,S,n,rr,Z,m,epsilon_z,Vfi,Uz,R2,vector1,neff1,col1)
% 
% % 磁场分布
% col2=size(vector2,2);  % 得到列数
% Hr(T,N,S,n,rr,r5,vector2,col2)
% Hfi(T,N,S,n,rr,r5,vector2,col2)
% Hz(k0,T,N,S,n,rr,r5,m,epsilon_r,epsilon_fi,Ufi,Vz,R1,vector2,neff2,col2)

legend('E_\Phi','H_r','H_z')
% 画出纤芯以及高低折射率的边界
 y=-1:3*10^-3:1;
 plot(rcore,y,'-b')
 hold on
 axis([h/2 Thick -1 1])
 title('HE_1_1模的场分布')
 xlabel('radius')
 ylabel('归一化幅度')
toc