bl_fdtd1_pm_mur.m

标量FDTD模拟平面光波导内部TE波的传播

% 标量FDTD分析平面光波导（只存在低阶模，沿z轴传播，纵轴为y轴）mur吸收边界条件(把精确吸收边界条件中的根号部分进行taylor展开，然后取前两项)
% 折射率只在x方向有变化，光沿z方向传播，
% 波导的形状及结构均与y无关,三维可变成二维问题来处理，平面波导中可存在TE模或TM模，以TE模为例分析
% TE波:Ey,Hx,Hz
% 利用标量FDTD法分析，中心差分近似，利用导模条件取m=0只有基模TE0
% 平面光源，球面光源，高斯光源，sin*sin,cos*sin分别激励
clc
clear

n1=1.563;  % 芯层折射率
n2=1.550;  % 包层和衬底折射率

c=3*10^8;  % 真空光速 
lam=1.30*10^(-6); % 波长 
epsz=1/(4*pi*9*10^9);   %  真空介电常数
u0=4*pi*10^(-7); % 真空磁导率

dels=lam/(20*n1); % 保证计算精度，ds=lma/(20*n1)  
vmax=1/(sqrt(u0*epsz)*min(n1,n2));
delt=dels/(2*vmax);    %  稳定条件

N=fix(1*10^(-6)/dels)+26;  %  空间步长(i:50)
NN=fix(1.0397*10^(-5)/dels);  %  (k:250)
T=1000;    % 时间步长

ey1=zeros(N,NN);  %  初始化
ey2=zeros(N,NN);
ey3=zeros(N,NN);
for i=1:N
    for j=1:NN
        if j>=13&j<=37
            n(i,j)=n2;
        else
            n(i,j)=n1;
        end
    end
end

for t=1:T
    for i=1:N
        if i>=13&i<=37
%             ey2(i,1)=sin(2*pi*c/lam*(t+1)*delt);% 芯层端面上入射的平面谐振光源
            ey2(i,1)=exp(-(2*(i-25)*dels/(1*10^(-6)/2))^2)*sin(2*pi*c/lam*(t+1)*delt);%高斯光源
% %             ey2(i,1)=sin(pi*i/48)*sin(2*pi*c/lam*(t+1)*delt); % sin*sin光源
%             ey2(i,1)=cos(pi*i/48)*sin(2*pi*c/lam*(t+1)*delt); % cos*sin光源
%            % 球面光源??
%             if (t+1)<2*sqrt((i-25)^2+81)/1.55
%                 ey2(i,1)=0;
%             else
%                 ey2(i,1)=6/sqrt((i-25)^2)*cos(0.05*pi*(2*sqrt((i-25)^2+81)...
%                      /1.563-(t+1)));
%                    ey2(i,1)=6/sqrt((i-25)^2)*cos(0.05*pi*(2*sqrt((i-25)^2+81)...
%                      /1.563-(t+1)));

%              end
         end 
    end
    for k=2:NN-1 
        for i=2:N-1
            if i>=13&i<=37
                p1(i,k)=delt^2/(epsz*u0*(dels^2)*n(i,j)^2);
                ey3(i,k)=p1(i,k)*(ey2(i-1,k)+ey2(i+1,k)+ey2(i,k-1)+ey2(i,k+1)...
                    -4*ey2(i,k))+2*ey2(i,k)-ey1(i,k);
            else
                p2(i,k)=delt^2/(epsz*u0*(dels^2)*n(i,j)^2);
                ey3(i,k)=p2(i,k)*(ey2(i-1,k)+ey2(i+1,k)+ey2(i,k-1)+ey2(i,k+1)...
                    -4*ey2(i,k))+2*ey2(i,k)-ey1(i,k);
                
            end
        end
    end
 % 边界点(右,上,下边界)mur吸收边界条件二阶近似
   km=NN;im=N;
   for i=2:im-1
       ey3(i,km)=-ey1(i,km-1)-1/3*(ey3(i,km-1)+ey1(i,km))+4/3*(ey2(i,km)...
            +ey2(i,km-1))+1/12*(ey2(i+1,km)-2*ey2(i,km)+ey2(i-1,km)...
            +ey2(i+1,km-1)-2*ey2(i,km-1)+ey2(i-1,km-1));
    end
    
    for k=2:km-1
        ey3(im,k)=-ey1(im-1,k)-1/3*(ey3(im-1,k)+ey1(im,k))+4/3*(ey2(im,k)...
            +ey2(im-1,k))+1/12*(ey2(im,k+1)-2*ey2(im,k)+ey2(im,k-1)...
            +ey2(im-1,k+1)-2*ey2(im-1,k)+ey2(im-1,k-1));
        ey3(1,k)=-ey1(2,k)-1/3*(ey3(2,k)+ey1(1,k))+4/3*(ey2(1,k)...
            +ey2(2,k))+1/12*(ey2(1,k+1)-2*ey2(1,k)+ey2(1,k-1)...
            +ey2(2,k+1)-2*ey2(2,k)+ey2(2,k-1));

    end
    % 角点(右边界上下两个角点）
    a1=48/(48^2+248^2)^(1/2);
    b1=248/(48^2+248^2)^(1/2);
    En1=(1-a1)*(1-b1)*ey2(im,km)+(1-a1)*b1*ey2(im,km-1)+a1*(1-b1)*ey2(im-1,km)...
        +a1*b1*ey2(im-1,km-1);
    ey3(im,km)=0.9960333*En1;
    En2=(1-a1)*(1-b1)*ey1(1,km)+(1-a1)*b1*ey1(1,km-1)+a1*(1-b1)*ey1(1,km)...
        +a1*b1*ey1(2,km-1);
    ey3(1,km)=0.9960333*En1;
    ey1=ey2;
    ey2=ey3;
    mesh(ey3);
    axis([1 NN 1 N -1 1]);
    drawnow;
end
          figure(3)
          hold on
%           axis([13 37 0 1]);
          plot(1:50,abs(ey3(1:50,3)))
%           for i=1:50;
% if i<=13
% ey3(i,3)=0;
% end
% if i>=37
% ey3(i,3)=0;
% end
% end
%  plot(1:50,abs(ey3(1:50,3)))
% 结论：
% 不同光源作用与同一个平面波导，波导模式与激励的类型无关，
% 只由波导的结构参数与光源的频率(或波长)决定。
% 达到稳定后五种情况下的稳态基模模场分布相同。