tm_fdtd.m

二维TM波的FDTD算法

%设置参数
clear
mu=4*pi*10.^(-7);%真空中的磁导率
epp=(1/(36*pi))*10.^(-9);% 真空介电常数
T=15*10.^(-12);
t0=3*T;
c=3*10.^8;
f=1/(2*T);%频率
lameda=c/f;%波长
dd=lameda/15;%网格分辨率
dt=dd/(2*c);%时间步长
%吸收边界条件迭代方程中的系数
a=(c*dt-dd)/(c*dt+dd);
aa=2*dd/(c*dt+dd);
aaa=((c*dt).^2)/(2*dd*(c*dt+dd));
%
%角点迭代方程中的系数
bb=(c*dt-sqrt(2)*dd)/(c*dt+sqrt(2)*dd);
%
nt=400;
nx=200;
ny=nx;
%---------------------------------------------------初始化--------------------------------------------------------%
for i=1:nx
    for j=1:ny;
        Hx(i,j)=0;
        Hy(i,j)=0;
        Ez(i,j)=0;
        cz(i,j)=0;% n时刻的Ez
        dz(i,j)=0;% n-1时刻的Ez
        hx(i,j)=0;% n-1/2时刻的Hx
        hy(i,j)=0;% n-1/2时刻的Hy
    end
end
%
%------------------------------------------------电场和磁场迭代---------------------------------------------------%
for t=1:nt  %时间步数
    i=1:nx;
    j=1:ny;
    mesh(i,j,Ez);%输出波形
    axis([1,nx,1,ny,-0.05,0.05]);
    pause(dt);
%磁场迭代
    for i=1:nx
        for j=1:ny-1 
           Hx(i,j)=hx(i,j)+(dt/(mu*dd))*(cz(i,j)-cz(i,j+1)); 
        end
    end
    for i=1:nx-1
        for j=1:ny 
            Hy(i,j)=hy(i,j)+(dt/(mu*dd))*(cz(i+1,j)-cz(i,j));
        end
    end
 %
    pulse=exp(-(((t*dt-t0)/T).^2));%激励源高斯脉冲
 %电场迭代
    for i=2:nx
        for j=2:ny-1
            if (i==nx/2&j==ny/2)
            Ez(i,j)=cz(i,j)+(dt/(epp*dd))*(Hy(i,j)-Hy(i-1,j)+Hx(i,j-1)-Hx(i,j))+pulse;
            else
            Ez(i,j)=cz(i,j)+(dt/(epp*dd))*(Hy(i,j)-Hy(i-1,j)+Hx(i,j-1)-Hx(i,j));
            end
        end
    end
 %
 %---------------------------------------------------设置Mur的二阶吸收边界条件----------------------------------------%
 %边界x=1
    for j=2:ny-1
        Ez(1,j)=-dz(2,j)+a*(Ez(2,j)+dz(1,j))+aa*(cz(1,j)+cz(2,j))+aaa*(cz(1,j+1)-2*cz(1,j)+cz(1,j-1)+cz(2,j+1)-2*cz(2,j)+cz(2,j-1));
    end
 %边界x=nx
    for j=2:ny-1
        Ez(nx,j)=-dz(nx-1,j)+a*(Ez(nx-1,j)+dz(nx,j))+aa*(cz(nx,j)+cz(nx-1,j))+aaa*(cz(nx,j+1)-2*cz(nx,j)+cz(nx,j-1)+cz(nx-1,j+1)-2*cz(nx-1,j)+cz(nx-1,j-1));
    end
  %边界y=1  
    for i=2:nx-1
        Ez(i,1)=-dz(i,2)+a*(Ez(i,2)+dz(i,1))+aa*(cz(i,1)+cz(i,2))+aaa*(cz(i+1,1)-2*cz(i,1)+cz(i-1,1)+cz(i+1,2)-2*cz(i,2)+cz(i-1,2));
    end
  %边界y=ny
    for i=2:nx-1
        Ez(i,ny)=-dz(i,ny-1)+a*(Ez(i,ny-1)+dz(i,ny))+aa*(cz(i,ny)+cz(i,ny-1))+aaa*(cz(i+1,ny)-2*cz(i,ny)+cz(i-1,ny)+cz(i+1,ny-1)-2*cz(i,ny-1)+cz(i-1,ny-1));
    end
  %---------------------------------------------------------角点处理------------------------------------------------------------%
    Ez(1,1)=cz(2,2)+bb*(Ez(2,2)-cz(1,1));%角点x=1,y=1
    Ez(1,ny)=cz(2,ny-1)+bb*(Ez(2,ny-1)-cz(1,ny));%角点x=1,y=ny
    Ez(nx,1)=cz(nx-1,2)+bb*(Ez(nx-1,2)-cz(nx,1));%角点x=nx,y=1
    Ez(nx,ny)=cz(nx-1,ny-1)+bb*(Ez(nx-1,ny-1)-cz(nx,ny));%角点x=nx,y=ny
  %
    dz=cz;
    cz=Ez;
    hx=Hx;
    hy=Hy;
end