matlab模拟te波abc边界条件的程序.txt

MATLAB在ABC边界上的FDTD程序仿真程序。

clear;clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
IE=100; 
JE=100;% 网格数量
L=400; % 迭代次数
f=1.6e9; % 频率
dx=.01; % 网格尺寸
dt=dx/2/3e8; % 时间间隔
ic=3;%平面波源位置

e0=8.8e-12; % 真空介电常数
m0=pi*4e-7; % 真空磁导率
epsilon=1; % 相对介电常数
sigma=0; % 电导率

ia=7;
ib=IE-ia+1;
ja=7;
jb=JE-ja+1;

for i=1:IE 
    for j=1:JE
        ga(i,j)=1/(epsilon+sigma*dt/e0); % 求和参量
        gb(i,j)=sigma*dt/e0; % 求和参量
    end
end
ez_inc=zeros(1,IE);
hx_inc=zeros(1,JE);

ez_low1=0;
ez_low2=0;
ez_low3=0;
ez_low4=0;

dz=zeros(IE,JE); % z方向电荷密度
ez=zeros(IE,JE); % z方向电场
iz=zeros(IE,JE); % z方向电场求和参量
hx=zeros(IE,JE); % x方向磁场
hy=zeros(IE,JE); % y方向磁场
ihx=zeros(IE,JE); % x方向磁场参量
ihy=zeros(IE,JE); % y方向磁场参量

%%%%%%%%%%%%%%%%%%%%%%%%%迭带求解电磁场%%%%%%%%%%%%%%%%%
for T=1:L
    for j=2:JE
        ez_inc(j)=ez_inc(j)+0.5*(hx_inc(j-1)-hx_inc(j));
    end
    
    ez_inc(1)=ez_low2;
    ez_low2=ez_low1;
    ez_low1=ez_inc(2);
    ez_inc(JE)=ez_low3;
    ez_low3=ez_low4;
    ez_low4=ez_inc(JE-1);  
%%%%%%%%%%%%%% 电荷密度dz%%%%%%%%%%%%%%%%%%%%%%%%%
    for i=2:IE-1 
        for j=2:JE
            dz(i,j)=dz(i,j)+0.5*( hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1));
        end
    end  
    
    source=sin(2*pi*f*dt*T); 
    ez_inc(ic)=source; 
    
    for i=ia:ib
        dz(i,ja)=dz(i,ja)+0.5*hx_inc(ja-1);
        dz(i,jb)=dz(i,jb)-0.5*hx_inc(jb);
    end
%%%%%%%%%%%%%% 电场ez%%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:IE
        for j=1:JE
            ez(i,j)=ga(i,j)*( dz(i,j)-iz(i,j) );
            iz(i,j)=iz(i,j)+gb(i,j)*ez(i,j) ;
        end
    end
%%%%%%%%%%%%%%%%%%%%边界电场%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:IE
        ez(1,j)=0;
        ez(IE,j)=0;
    end
    for i=1:JE
        ez(i,1)=0;
        ez(i,JE)=0;
    end
    
    for j=1:JE-1
        hx_inc(j)=hx_inc(j)+0.5*(ez_inc(j)-ez_inc(j+1));
    end
%%%%%%%%%%%%%%%%%%%%%%%%磁场Hx%%%%%%%%%%%%%%%%%%%%%%
    for i=1:IE
        for j=1:JE-1
            hx(i,j)=hx(i,j)+0.5*(ez(i,j)-ez(i,j+1));
        end
    end

    for i=ia:ib
        hx(i,ja-1)=hx(i,ja-1)+0.5*ez_inc(ja);
        hx(i,jb)=hx(i,jb)-0.5*ez_inc(jb);
    end
%%%%%%%%%%%%%%%%%磁场Hy%%%%%%%%%%%%%%%%%%%
    for i=1:IE-1 
        for j=1:JE
            hy(i,j)=hy(i,j)+0.5*(ez(i+1,j)-ez(i,j));
        end
    end

    for j=ja:jb;
        hy(ia-1,j)=hy(ia-1,j)-0.5*ez_inc(j);
        hy(ib,j)=hy(ib,j)+0.5*ez_inc(j);
    end
%%%%%%%%%%%%%%%%%%绘图%%%%%%%%%%%%%%%%%%
    surf(ez)
    axis([0 100 0 100 -2 2])
    shading interp;
    %view(0,90);
    drawnow;
end