用fdtd方法分析二维tm波在正弦波源激励情况下，正方形mur吸收边界的性能.txt

用FDTD方法分析二维TM波在正弦波源激励情况下

程序1：用FDTD方法分析二维TM波在正弦波源激励情况下，正方形Mur吸收边界的性能
clear;clc;
N=100;
c=3e8;
f=1.5*10^(9);
lambda=c/f
eb=8.85e-12;
u=4*pi*10^-7;
Ez1=zeros(N+1,N+1);
Hx1=zeros(N+1,N);
Hy1=zeros(N,N+1);
Ez2=zeros(N+1,N+1);
Hx2=zeros(N+1,N);
Hy2=zeros(N,N+1);
deltaL=lambda/10;
deltaT=deltaL/(2*c);
 
for  n=1:200 
    ca=1;
    cb=deltaT/eb;
    cp=1;
    cq=deltaT/u;%迭代系数
    It=sin(2*pi*f*n*deltaT);
    %非边界场值迭代
    clear i j;
    for i=1:(N+1)
        for j=1:N
            Hx2(i,j)=cp*Hx1(i,j)-cq*(Ez1(i,j+1)-Ez1(i,j))/deltaL;
        end
    end
    clear i j;
    for i=1:N
        for j=1:(N+1)
            Hy2(i,j)=cp*Hy1(i,j)+cq*(Ez1(i+1,j)-Ez1(i,j))/deltaL;
        end
    end
    clear i j;
    for i=2:N
        for j=2:N
Ez2(i,j)=ca*Ez1(i,j)+cb*((Hy2(i,j)-Hy2(i-1,j))/deltaL-(Hx2(i,j)-Hx2(i,j-1))/deltaL);
        end
    end
    clear i j;
    i=50;j=50;
    Ez2(i,j)=Ez1(i,j)+It;%激励加入
    %边界非角点场值迭代
%左边界
    clear i j;
    i=1;
    for j=2:N
Ez2(i,j)=Ez1(i+1,j)+(c*deltaT-deltaL)*(Ez2(i+1,j)-Ez1(i,j))/(c*deltaT+deltaL);
    end
%右边界
    clear i j;
    i=N+1;
    for j=2:N
Ez2(i,j)=Ez1(i-1,j)+(c*deltaT-deltaL)*(Ez2(i-1,j)-Ez1(i,j))/(c*deltaT+deltaL);
    end
%上边界
    clear i j;
    j=N+1;
    for i=2:N
Ez2(i,j)=Ez1(i,j-1)+(c*deltaT-deltaL)*(Ez2(i,j-1)-Ez1(i,j))/(c*deltaT+deltaL);
    end
%下边界
    clear i j;
    j=1;
    for i=2:N
Ez2(i,j)=Ez1(i,j+1)+(c*deltaT-deltaL)*(Ez2(i,j+1)-Ez1(i,j))/(c*deltaT+deltaL);
    end
%角点的迭代
    clear i j;
Ez2(1,1)=Ez1(2,2)+(c*deltaT-sqrt(2)*deltaL)*(Ez2(2,2)-Ez1(1,1))/(c*deltaT+sqrt(2)*deltaL);
Ez2(N+1,N+1)=Ez1(N,N)+(c*deltaT-sqrt(2)*deltaL)*(Ez2(N,N)-Ez1(N+1,N+1))/(c*deltaT+sqrt(2)*deltaL);
Ez2(1,N+1)=Ez1(2,N)+(c*deltaT-sqrt(2)*deltaL)*(Ez2(2,N)-Ez1(1,N+1))/(c*deltaT+sqrt(2)*deltaL);
Ez2(N+1,1)=Ez1(N,2)+(c*deltaT-sqrt(2)*deltaL)*(Ez2(N,2)-Ez1(N+1,1))/(c*deltaT+sqrt(2)*deltaL);
%一次迭代结束，新得到的值变为下次迭代的旧值
    clear i j;
    for i=1:(N+1);
        for j=1:(N+1)
            Ez1(i,j)=Ez2(i,j);
        end
    end
    clear i j;
    for i=1:(N+1)
        for j=1:N
            Hx1(i,j)=Hx2(i,j);
        end
    end
     clear i j;
    for i=1:N
        for j=1:(N+1)
            Hy1(i,j)=Hy2(i,j);
        end
    end
%一次迭代结束
end
clear i j;
i=1:(N+1);
j=1:(N+1);
figure(1);
contour(i,j,Ez2(i,j));xlabel('x');ylabel('y');