upml_2d_tmmode.m

二维TM模电磁波FDTD仿真

%这里为仅有Hx,Hy,Ez时的极化情况
clear 
% 边界条件参数
% iebc=input('请输入UPML的层数： iebc=')           %PML thickness
% m=input('please input the order of UPML')        %order of PML 
iebc=8;
m=4;
%some physical constant
cc=2.99792458e8;            %speed of light in free space
muz=4.0*pi*1.0e-7;
epsz=1.0/(cc^2*muz);    

%simulator parameters
 ix =80;       %X cell number
 jy =80;      %Y Grid number
 ixb=ix+1;
 jyb=jy+1;
 ie=ix+2*iebc;         %for UPML surround calculate whole space grid
 je=jy+2*iebc;
 ib=ie+1;               %for some uses we need calculat X_size+1 .etc
 jb=je+1;         
 nmax=3000;         %number of time steps


%分配内存
%系数矩阵
DCP=zeros(ib,jb);
DCQ=zeros(ib,jb);
ECP=zeros(ib,jb);
ECQ=zeros(ib,jb);

BxCA=zeros(ib,jb);
BxCB=zeros(ib,jb);
HxCA=zeros(ib,jb);
HxCB1=zeros(ib,jb);
HxCB2=zeros(ib,jb);

ByCA=zeros(ib,jb);
ByCB=zeros(ib,jb);
HyCA=zeros(ib,jb);
HyCB1=zeros(ib,jb);
HyCB2=zeros(ib,jb);
%场分量
Bx=zeros(ib,jb);
By=zeros(ib,jb);
Hx=zeros(ib,jb);
Hy=zeros(ib,jb);
Dz=zeros(ib,jb);
Ez=zeros(ib,jb);
%相对介电常数矩阵
M_epsr=ones(ib,jb);
% IOtest1;
%建立入射脉冲源
freq  = 193*10^12;      
lambda=cc/freq;
omega=2.0*pi*freq;
% N     = 20 ;                %cell numbers per wavelength
% ds    =lambda/N ;  
dx=0.02e-6;        %space increment in the x_direction
dy=0.02e-6 ;                   %dy
                     
dt=dx/(2*cc);
rtau=1/freq; 
tau=rtau/dt; 
delay=3*tau; 

source=zeros(1,nmax); 
for n=1:7.0*tau 
  source(n)=sin(omega*(n-delay)*dt)*exp(-((n-delay)^2/tau^2)); 
end 
% for n=1:nmax
%   source(n)=sin(omega*(n)*dt); 
% end 
%以上假定各种参数都已经准备好,各种参数都采用完备的存储单元分配方式
%left side
 %Ez left 
e=epsz;
Consx=1/dx/((iebc*dx)^m)/(m+1);
Consy=1/dy/((iebc*dy)^m)/(m+1);
rmax =10^-8;
Omax=-log(rmax/100.0)*epsz*cc*(m+1)/(2*dx*iebc);
for ii=1:iebc
    for jj=iebc+1:iebc+jy
    eff=M_epsr(ii,jj);
    x1=(iebc-ii)*dx;
    x2=(iebc-ii+1)*dx;
    Ox=Omax*Consx*(x2^(m+1)-x1^(m+1));     %应用半网格点处的Ox值和未平均时的eps作为eff
    Oy=0;
    DCP(ii,jj)=(2*e-Oy*dt)/(2*e+Oy*dt);
    DCQ(ii,jj)=2*e*dt/(2*e+Oy*dt);
    ECP(ii,jj)=(2*e-Ox*dt)/(2*e+Ox*dt);
    ECQ(ii,jj)=2*e/(2*e+Ox*dt)/eff/e;
    end
end

 %Hy left
for ii=2:iebc % <--+++++++++++++++++++++ for ii=2:iebc
    for jj=iebc+1:iebc+jy
    eff= (M_epsr(ii,jj)+M_epsr(ii+1,jj))/2;   %应用x方向上平均后的eps作为eff    <--++++++++++++ eff= (M_epsr(ii,jj)+M_epsr(ii-1,jj))/2;
    x1=(iebc-ii+0.5)*dx; % <--+++++++++++++++++  x1=(iebc-ii+0.5)*dx;
    x2=(iebc-ii+1.5)*dx; % <--+++++++++++++++++  x2=(iebc-ii+1.5)*dx;
    Oy=0;
    Ox=Omax*Consx*(x2^(m+1)-x1^(m+1));      %应用整网格点处的Ox值
    ByCA(ii,jj)=1;
    ByCB(ii,jj)=dt;
    HyCA(ii,jj)=(2*e-Ox*dt)/(2*e+Ox*dt);
    HyCB1(ii,jj)=(2*e+Oy*dt)/(2*e+Ox*dt)/muz;
    HyCB2(ii,jj)=(2*e-Oy*dt)/(2*e+Ox*dt)/muz;
    end
end

 %Hx left
%对front 和left的交接处做特殊处理
for ii=1:iebc
    for jj=iebc+1:iebc+jy
    eff=(M_epsr(ii,jj)+M_epsr(ii,jj+1))/2 ;        %Y方向平均后的eps
%     Omax=(m+1)/sqrt(eff)/150/pi/dx;
    x1=(iebc-ii)*dx; % <--++++++++++++++++
    x2=(iebc-ii+1)*dx; % <--+++++++++++++++++
    Ox=Omax*Consx*(x2^(m+1)-x1^(m+1));          %半网格处的Ox
    if jj==iebc+1
        
        Oy=Omax*Consy*(0.5*dy)^(m+1);
    else
        Oy=0;
    end
    BxCA(ii,jj)=1;
    BxCB(ii,jj)=dt;
    HxCA(ii,jj)=(2*e-Oy*dt)/(2*e+Oy*dt);
    HxCB1(ii,jj)=(2*e+Ox*dt)/(2*e+Oy*dt)/muz;
    HxCB2(ii,jj)=(2*e-Ox*dt)/(2*e+Oy*dt)/muz;
    end
end
%modified to this line 9:42 pm 4.10
%front
 %Ez front
 for ii=1:2*iebc+ix
     for jj=1:iebc
         eff=M_epsr(ii,jj);   %未平均化的eps

         if (ii<=iebc)
           
             x1  =(iebc-ii)*dx;
             x2  =(iebc-ii+1)*dx;
             Ox  =Omax*Consx*(x2^(m+1)-x1^(m+1));   %半网格处的值
         else
             if ii>iebc+ix
                 x1=(ii-iebc-ix)*dx;
                 x2=(ii-iebc-ix-1)*dx;
                 Ox  =Omax*Consx*(x1^(m+1)-x2^(m+1)); 
             else
                 Ox=0;
             end
         end
%        
         y1=(iebc-jj)*dy;
         y2=(iebc-jj+1)*dy;
         Oy =Omax*Consy*(y2^(m+1)-y1^(m+1));   %半网格处的值
         DCP(ii,jj)=(2*e-Oy*dt)/(2*e+Oy*dt);
         DCQ(ii,jj)=2*e*dt/(2*e+Oy*dt);
         ECP(ii,jj)=(2*e-Ox*dt)/(2*e+Ox*dt);
         ECQ(ii,jj)=2*e/(2*e+Ox*dt)/eff/e;
         
     end
 end
 
 %Hy front
 %对左右边界上的Hy赋0
 for ii=2:2*iebc+ix
     for jj=1:iebc
         eff=(M_epsr(ii,jj)+M_epsr(ii-1,jj))/2;     %应用X方向平均后的eps

         y1=(iebc-jj)*dy;
         y2=(iebc-jj+1)*dy;
         Oy =Omax*Consy*(y2^(m+1)-y1^(m+1));    %应用半网格处的Oy
         %分为三个区域,角区域应用整网格处的Ox并对交接处特殊处理Ox

         if ii==iebc+1
         Ox= Omax*Consx*(0.5*dx)^(m+1);     %边界上特殊处理的Ox
         else
             if ii==iebc+ix+1
             Ox= Omax*Consx*(0.5*dx)^(m+1);      %边界上特殊处理的Ox
             else
                 if  ii<iebc+1
                     x1=(iebc-ii+0.5)*dx;
                     x2=(iebc-ii+1.5)*dx;
                     Ox=Omax*Consx*(x2^(m+1)-x1^(m+1));                         %用整网格上的Ox
                 else
                     if ii>iebc+ix+1
                         x1=(ii-iebc-ix-0.5)*dx;
                         x2=(ii-iebc-ix-1.5)*dx;
                         Ox=Omax*Consx*(x1^(m+1)-x2^(m+1));
                     else
                         Ox=0;
                     end
                 end
             end
         end
         ByCA(ii,jj)=1;
         ByCB(ii,jj)=dt;
         HyCA(ii,jj)=(2*e-Ox*dt)/(2*e+Ox*dt);
         HyCB1(ii,jj)=(2*e+Oy*dt)/(2*e+Ox*dt)/muz;
         HyCB2(ii,jj)=(2*e-Oy*dt)/(2*e+Ox*dt)/muz;
     end
 end
 
 %Hx front
 %对下边界所有的Hx赋0
 for ii=1:2*iebc+ix
     for jj=2:iebc
         eff=(M_epsr(ii,jj)+M_epsr(ii,jj-1))/2 ; %应用y方向平均后的eps
%          Omax=(m+1)/sqrt(eff)/150/pi/dy;
         y1=(iebc-jj+1.5)*dy;
         y2=(iebc-jj+0.5)*dy;
         Oy  =Omax*Consy*(y1^(m+1)-y2^(m+1));  %应用整网格点处的Oy
         
%          Omax=(m+1)/sqrt(eff)/150/pi/dx;
         if  (ii<=iebc)
         x1=(iebc-ii)*dx;
         x2=(iebc-ii+1)*dx;
         Ox=Omax*Consx*(x2^(m+1)-x1^(m+1));   %应用半网格点处的Ox
         else
             if ii>iebc+ix
                 x1=(ii-iebc-ix)*dx;
                 x2=(ii-iebc-ix-1)*dx;
                 Ox=Omax*Consx*(x1^(m+1)-x2^(m+1));   %应用半网格点处的Ox
             else
                 Ox=0;
             end
         end
        BxCA(ii,jj)=1;
        BxCB(ii,jj)=dt;
        HxCA(ii,jj)=(2*e-Oy*dt)/(2*e+Oy*dt);
        HxCB1(ii,jj)=(2*e+Ox*dt)/(2*e+Oy*dt)/muz;
        HxCB2(ii,jj)=(2*e-Ox*dt)/(2*e+Oy*dt)/muz;
     end
 end
 
 %right     
 %Ez right
 for ii=iebc+ix+1:iebc+ix+iebc
    for jj=iebc+1:iebc+jy
    Oy=0;
    eff=M_epsr(ii,jj);    %未平均时的eps作为eff
    x1=(ii-iebc-ix)*dx;
    x2=(ii-iebc-ix-1)*dx;
    Ox=Omax*Consx*(x1^(m+1)-x2^(m+1));     %应用半网格点处的Ox值
    e=epsz;
    DCP(ii,jj)=(2*e-Oy*dt)/(2*e+Oy*dt);
    DCQ(ii,jj)=2*e*dt/(2*e+Oy*dt);
    ECP(ii,jj)=(2*e-Ox*dt)/(2*e+Ox*dt);
    ECQ(ii,jj)=2*e/(2*e+Ox*dt)/eff/e;
    end
end

%Hy right
%将右边界置0
for ii=iebc+ix+1:iebc+ix+iebc
    for jj=iebc+1:iebc+jy
        Oy=0;
        eff=(M_epsr(ii,jj)+M_epsr(ii-1,jj))/2;     %在x方向上作平均
        %对CS 和right的交接处特殊处理
%         Omax=(m+1)/sqrt(eff)/150/pi/dx;
        if ii==iebc+ix+1
            Ox=Omax*Consx*(0.5*dx)^(m+1);
        else
            x1=(ii-iebc-ix-0.5)*dx;
            x2=(ii-iebc-ix-1.5)*dx;
            Ox=Omax*Consx*(x1^(m+1)-x2^(m+1));  %用整网格上的Ox
        end
         ByCA(ii,jj)=1;
         ByCB(ii,jj)=dt;
         HyCA(ii,jj)=(2*e-Ox*dt)/(2*e+Ox*dt);
         HyCB1(ii,jj)=(2*e+Oy*dt)/(2*e+Ox*dt)/muz;
         HyCB2(ii,jj)=(2*e-Oy*dt)/(2*e+Ox*dt)/muz;
    end
end

%Hx right
%对front和right交接处做特殊处理

for ii=iebc+ix+1:iebc+ix+iebc
    for jj=iebc+1:iebc+jy
        eff=(M_epsr(ii,jj)+M_epsr(ii,jj-1))/2;  %在y方向上作平均
%         Omax=(m+1)/sqrt(eff)/150/pi/dx;
        x1=(ii-iebc-ix)*dx;
        x2=(ii-iebc-ix-1)*dx;
        Ox =Omax*Consx*(x1^(m+1)-x2^(m+1));  %半网格点上的值
%         Omax=(m+1)/sqrt(eff)/150/pi/dy;
        if jj==iebc+1
             Oy =Omax*Consy*(0.5*dy)^(m+1);
        else
             Oy =0;
        end
    BxCA(ii,jj)=1;
    BxCB(ii,jj)=dt;
    HxCA(ii,jj)=(2*e-Oy*dt)/(2*e+Oy*dt);
    HxCB1(ii,jj)=(2*e+Ox*dt)/(2*e+Oy*dt)/muz;
    HxCB2(ii,jj)=(2*e-Ox*dt)/(2*e+Oy*dt)/muz;
    end
end