fdtd.m

一个参照葛德彪书编写的3D FDTD源文件

obj=ones(51,51,51);                           %存放空间媒质编号，包括散射体，迭代次数
Ntime=100;
freq=1e9;                                        %频率上限
Nmedia=2;media=zeros(Nmedia,4);                  %可存Nmedia种介质的电磁参数表
media=[1,1,0,0;1,1,3.72e7,0];                    %[epslon_r,miu_r,sigma,sigma_m]

Tin=1/freq;                                      %入射波时间常数
Vin=1;t0=1e-9;
Phi=0;                                           %平面波入射角
WL=40;                                           %数值波长，wavel有40个网格
Ninwave=100;
%---------不变参数-------------
miu0=4e-7*pi;epslon0=8.854187818e-12;            %真空电磁参数
vc=299792458;wavel=vc/freq;                      %光速和波长设置
Z0=sqrt(miu0/epslon0);                           %波阻抗
%------------------------------
%------导入模型到odj离散空间----------
%readmodel;
%[xs,ys,zs]=size(mmd);
%obj(101-floor(xs/2):101-floor(xs/2)+xs-1,101-floor(ys/2):101-floor(ys/2)+ys-1,101-floor(zs/2):101-floor(zs/2)+zs-1)=mmd+1;
obj(24:28,24:28,24:28)=obj(24:28,24:28,24:28)+1;
%----------导入模型结束---------------

Izmin=11;Izmax=41;%总场区边界
Ismin=6;Ismax=46;%散射场区边界
Fzmin=9;Fzmax=43;%远场外推边界
TEM_flag=1;        %TE（＝1） or TM（＝2）
Ez=zeros(51,51,51);Ex=Ez;Ey=Ez;Hx=Ez;Hy=Ez;Hz=Ez;%存放电磁场的空间
Ein=zeros(1,Ninwave);Hin=Ein;%存放1-D FDTD入射电磁场的空间
CA=zeros(Nmedia,Nmedia,Nmedia,Nmedia);CB=CA;%电场迭代系数
CP=zeros(Nmedia,Nmedia);CQ=CP;%磁场迭代系数
delta=wavel/WL;
dt=delta/(2*vc);
%---申请吸收边界暂存空间(参考葛德彪P63)----
%6个面，每个面2个切向吸收E分量，每个分量需4层暂存层空间(1:Pn, 2:Pn-1 ,3:Qn, 4:Qn-1 )
Ab1x=zeros(51,4,51);Ab1z=Ab1x;      %1面x,z分量,4层
Ab2x=zeros(51,4,51);Ab2z=Ab2x;      %2面x,z分量,4层
Ab3x=zeros(51,51,4);Ab3y=Ab3x;      %3面x,y分量,4层
Ab4x=zeros(51,51,4);Ab4y=Ab4x;      %4面x,y分量,4层
Ab5z=zeros(4,51,51);Ab5y=Ab5z;      %5面z,y分量,4层
Ab6z=zeros(4,51,51);Ab6y=Ab6z;      %6面z,y分量,4层
%---吸收边界暂存空间 end----

%计算迭代系数
for i=1:Nmedia
    for j=1:Nmedia
        for k=1:Nmedia
            for l=1:Nmedia
                epslon_eff=0.25*(media(i,1)+media(j,1)+media(k,1)+media(l,1));
                sigma_eff=0.25*(media(i,3)+media(j,3)+media(k,3)+media(l,3));
                CA(i,j,k,l)=(1-sigma_eff*dt/(2*epslon0*epslon_eff))/(1+sigma_eff*dt/(2*epslon0*epslon_eff));
                CB(i,j,k,l)=(dt/(epslon0*delta))/(epslon_eff+sigma_eff*dt/(2*epslon0));
            end
        end
        miu_eff=0.5*(media(i,2)+media(j,2));
        sigmam_eff=0.5*(media(i,4)+media(j,4));
        CP(i,j)=(1-sigmam_eff*dt/(2*miu0*miu_eff))/(1+sigmam_eff*dt/(2*miu0*miu_eff));
        CQ(i,j)=(dt/miu0/delta)/(miu_eff+sigmam_eff*dt/(2*miu0));
    end
end
FE=dt/epslon0/delta;FH=dt/miu0/delta; %入射波1-D FDTD迭代系数
EBin=zeros(1,4);%1-D FDTD吸收边界变量
%====================================prepare work end=====================================================
%
%=========================================Main loop=======================================================
for Nt=1:Ntime
    Nt
    A(Nt)=Ez(26,9,26)
    %---------creat add_in wave using 1-D FDTD----------------
    for i=1:Ninwave-1
        Hin(i)=Hin(i)-FH*(Ein(i+1)-Ein(i));
    end
    for i=2:Ninwave
        Ein(i)=Ein(i)-FE*(Hin(i)-Hin(i-1));
    end
    Ein(5)=Vin*exp(-((Nt*dt-t0)/Tin)^2); %add sourse in 10th element
    %1-D FDTD absolute
    Ein(1)=EBin(4);EBin(4)=EBin(3);EBin(3)=Ein(2);
    Ein(Ninwave)=EBin(2);EBin(2)=EBin(1);EBin(1)=Ein(Ninwave-1);
    %-----Ein的函数表达-----
    %for i=1:Ninwave
    %   Ein(i)=Vin*exp(-((delta*i-vc*dt*(Nt-300))/(Tin*vc))^2);
    %end
    %-------------inwave 1-D FDTD end-------------------------
    %
    %-----start to caculate incident wave E and H feild on each connect plane-----
    %           --------
    %         /  3,4  y/|
    %        /___x____/ |  看的见得在前
    %      5-|z  1,2  |6/        
    %        |________|/
    %申请空间存放表面 切向电场分量和磁场分量,算出个表面的E，除以Z0就是H，乘以角度就是分量
    Ein34=zeros(51,51,1);%3,4面，由于入射波在Z向均匀，所以3，4面切向场相同共用一个
    Ein1=zeros(51,1,51);Ein2=Ein1;
    Ein5=zeros(1,51,51);Ein6=Ein5;
    for i=Izmin:Izmax
        for j=Izmin:Izmax
            %--------3 4 面---------
            Tty=i*sin(Phi)+j*cos(Phi);  Ttz=floor(Tty);  Tyu=Tty-floor(Tty);          
            Ein34(i,j,1)=Ein(Ttz)+Tyu*(Ein(Ttz+1)-Ein(Ttz));
            %---------1 面----------
            Tty=i*sin(Phi)+Izmin*cos(Phi);  Ttz=floor(Tty);  Tyu=Tty-floor(Tty);          
            Ein1(i,1,j)=Ein(Ttz)+Tyu*(Ein(Ttz+1)-Ein(Ttz));
            %---------2 面----------
            Tty=i*sin(Phi)+Izmax*cos(Phi);  Ttz=floor(Tty);  Tyu=Tty-floor(Tty);          
            Ein2(i,1,j)=Ein(Ttz)+Tyu*(Ein(Ttz+1)-Ein(Ttz));
            %---------5 面----------
            Tty=Izmin*sin(Phi)+j*cos(Phi);  Ttz=floor(Tty);  Tyu=Tty-floor(Tty);          
            Ein5(1,j,i)=Ein(Ttz)+Tyu*(Ein(Ttz+1)-Ein(Ttz));
            %---------6 面----------
            Tty=Izmax*sin(Phi)+j*cos(Phi);  Ttz=floor(Tty);  Tyu=Tty-floor(Tty);          
            Ein6(1,j,i)=Ein(Ttz)+Tyu*(Ein(Ttz+1)-Ein(Ttz));
        end
    end
    %------end of caculate incident wave E and H feild on each connect plane------
    %
    %--------------------------------------------------start E FDTD----------------------------------------------------
    for i=Ismin+1:Ismax-1
        for j=Ismin+1:Ismax-1
            for k=Ismin+1:Ismax-1
                T1=Hy(i,j,k)-Hy(i-1,j,k)-Hx(i,j,k)+Hx(i,j-1,k);
                Ez(i,j,k)=CA( obj(i,j,k),obj(i-1,j,k),obj(i,j,k),obj(i,j-1,k) ) * Ez(i,j,k) + CB( obj(i,j,k),obj(i-1,j,k),obj(i,j,k),obj(i,j-1,k) ) * T1;
                T1=Hz(i,j,k)-Hz(i,j-1,k)-Hy(i,j,k)+Hy(i,j,k-1);
                Ex(i,j,k)=CA( obj(i,j,k),obj(i-1,j,k),obj(i,j,k),obj(i,j-1,k) ) * Ex(i,j,k) + CB( obj(i,j,k),obj(i-1,j,k),obj(i,j,k),obj(i,j-1,k) ) * T1;
                T1=Hx(i,j,k)-Hx(i,j,k-1)-Hz(i,j,k)+Hz(i-1,j,k);
                Ey(i,j,k)=CA( obj(i,j,k),obj(i-1,j,k),obj(i,j,k),obj(i,j-1,k) ) * Ey(i,j,k) + CB( obj(i,j,k),obj(i-1,j,k),obj(i,j,k),obj(i,j-1,k) ) * T1;
            end
        end
    end
    %-----连接边界面引入Ez分量-----
    Ez(Izmin+1:Izmax-1,Izmin,Izmin+1:Izmax-1)=Ez(Izmin+1:Izmax-1,Izmin,Izmin+1:Izmax-1)+Ein1(Izmin+1:Izmax-1,1,Izmin+1:Izmax-1)/Z0*FE*cos(Phi);%面1
    Ez(Izmin+1:Izmax-1,Izmax,Izmin+1:Izmax-1)=Ez(Izmin+1:Izmax-1,Izmax,Izmin+1:Izmax-1)-Ein2(Izmin+1:Izmax-1,1,Izmin+1:Izmax-1)/Z0*FE*cos(Phi);%面2
    %----面3----
    Ex(Izmin+1:Izmax-1,Izmin+1:Izmax-1,Izmin)=Ex(Izmin+1:Izmax-1,Izmin+1:Izmax-1,Izmin)+Ein34(Izmin+1:Izmax-1,Izmin+1:Izmax-1,1)/Z0*FE*sin(Phi);
    Ey(Izmin+1:Izmax-1,Izmin+1:Izmax-1,Izmin)=Ey(Izmin+1:Izmax-1,Izmin+1:Izmax-1,Izmin)-Ein34(Izmin+1:Izmax-1,Izmin+1:Izmax-1,1)/Z0*FE*cos(Phi);
    %--面3 end--
    %----面4----
    Ex(Izmin+1:Izmax-1,Izmin+1:Izmax-1,Izmax)=Ex(Izmin+1:Izmax-1,Izmin+1:Izmax-1,Izmax)-Ein34(Izmin+1:Izmax-1,Izmin+1:Izmax-1,1)/Z0*FE*sin(Phi);
    Ey(Izmin+1:Izmax-1,Izmin+1:Izmax-1,Izmax)=Ey(Izmin+1:Izmax-1,Izmin+1:Izmax-1,Izmax)+Ein34(Izmin+1:Izmax-1,Izmin+1:Izmax-1,1)/Z0*FE*cos(Phi);
    %--面4 end--
    Ez(Izmin,Izmin+1:Izmax-1,Izmin+1:Izmax-1)=Ez(Izmin,Izmin+1:Izmax-1,Izmin+1:Izmax-1)-Ein5(1,Izmin+1:Izmax-1,Izmin+1:Izmax-1)/Z0*FE*sin(Phi);%面5
    Ez(Izmax,Izmin+1:Izmax-1,Izmin+1:Izmax-1)=Ez(Izmax,Izmin+1:Izmax-1,Izmin+1:Izmax-1)+Ein6(1,Izmin+1:Izmax-1,Izmin+1:Izmax-1)/Z0*FE*sin(Phi);%面6
    %----楞边引入入射波--Hx=H*cos;Hy=H*sin----
    %--//Z Ez分量----
    Ez(Izmin,Izmin,Izmin+1:Izmax-1)=Ez(Izmin,Izmin,Izmin+1:Izmax-1)-FE*(sin(Phi)-cos(Phi))*Ein1(Izmin,1,Izmin+1:Izmax-1)/Z0;       %1,5面交线
    Ez(Izmax,Izmin,Izmin+1:Izmax-1)=Ez(Izmax,Izmin,Izmin+1:Izmax-1)+FE*(sin(Phi)+cos(Phi))*Ein1(Izmax,1,Izmin+1:Izmax-1)/Z0;       %1,6面交线
    Ez(Izmax,Izmax,Izmin+1:Izmax-1)=Ez(Izmax,Izmax,Izmin+1:Izmax-1)+FE*(sin(Phi)-cos(Phi))*Ein2(Izmax,1,Izmin+1:Izmax-1)/Z0;       %2,6面交线
    Ez(Izmin,Izmax,Izmin+1:Izmax-1)=Ez(Izmin,Izmax,Izmin+1:Izmax-1)-FE*(sin(Phi)+cos(Phi))*Ein2(Izmin,1,Izmin+1:Izmax-1)/Z0;       %2,5面交线
    %--//X Ex分量----
    Ex(Izmin+1:Izmax-1,Izmin,Izmin)=Ex(Izmin+1:Izmax-1,Izmin,Izmin)+FE*sin(Phi)*Ein34(Izmin+1:Izmax-1,Izmin,1)/Z0;       %1,3面交线
    Ex(Izmin+1:Izmax-1,Izmax,Izmin)=Ex(Izmin+1:Izmax-1,Izmax,Izmin)+FE*sin(Phi)*Ein34(Izmin+1:Izmax-1,Izmax,1)/Z0;       %2,3面交线
    Ex(Izmin+1:Izmax-1,Izmax,Izmax)=Ex(Izmin+1:Izmax-1,Izmax,Izmax)-FE*sin(Phi)*Ein34(Izmin+1:Izmax-1,Izmax,1)/Z0;       %2,4面交线
    Ex(Izmin+1:Izmax-1,Izmin,Izmax)=Ex(Izmin+1:Izmax-1,Izmin,Izmax)-FE*sin(Phi)*Ein34(Izmin+1:Izmax-1,Izmin,1)/Z0;       %1,4面交线
    %--//Y Ey分量----
    Ey(Izmin,Izmin+1:Izmax-1,Izmin)=Ey(Izmin,Izmin+1:Izmax-1,Izmin)-FE*cos(Phi)*Ein34(Izmin,Izmin+1:Izmax-1,1)/Z0;       %5,3面交线
    Ey(Izmax,Izmin+1:Izmax-1,Izmin)=Ey(Izmax,Izmin+1:Izmax-1,Izmin)+FE*cos(Phi)*Ein34(Izmax,Izmin+1:Izmax-1,1)/Z0;       %6,3面交线
    Ey(Izmax,Izmin+1:Izmax-1,Izmax)=Ey(Izmax,Izmin+1:Izmax-1,Izmax)+FE*cos(Phi)*Ein34(Izmax,Izmin+1:Izmax-1,1)/Z0;       %6,4面交线
    Ey(Izmin,Izmin+1:Izmax-1,Izmax)=Ey(Izmin,Izmin+1:Izmax-1,Izmax)-FE*cos(Phi)*Ein34(Izmin,Izmin+1:Izmax-1,1)/Z0;       %5,4面交线
    %-------------------------------------------------吸收边界-------------------------------------------------------
    %-------------------------------------6个面-----------------------------------------
    %---------------面1(y=Ismin,吸收Ez、Ex)----------------------
    for i=Ismin+2:Ismax-2
        for j=Ismin+2:Ismax-2
            %--Ez--
            T1=Ez(i,Ismin+1,j)+Ab1z(i,2,j);
            T2=Ab1z(i,3,j)+Ab1z(i,1,j);
            T3=Ab1z(i+1,1,j)+Ab1z(i-1,1,j)+Ab1z(i,1,j+1)+Ab1z(i,1,j-1) + Ab1z(i+1,3,j)+Ab1z(i-1,3,j)+Ab1z(i,3,j+1)+Ab1z(i,3,j-1);
            Ez(i,Ismin,j)=-Ab1z(i,4,j)-T1/3+T2+T3/12;
            %--Ex--
            T1=Ex(i,Ismin+1,j)+Ab1x(i,2,j);
            T2=Ab1x(i,3,j)+Ab1x(i,1,j);
            T3=Ab1x(i+1,1,j)+Ab1x(i-1,1,j)+Ab1x(i,1,j+1)+Ab1x(i,1,j-1) + Ab1x(i+1,3,j)+Ab1x(i-1,3,j)+Ab1x(i,3,j+1)+Ab1x(i,3,j-1);
            Ex(i,Ismin,j)=-Ab1x(i,4,j)-T1/3+T2+T3/12;
        end
    end
    Ez(Ismin+1:Ismax-1,Ismin,Ismin+1)=Ab1z(Ismin+1:Ismax-1,3,Ismin+1)-(Ez(Ismin+1:Ismax-1,Ismin+1,Ismin+1)-Ab1z(Ismin+1:Ismax-1,1,Ismin+1))/3;
    Ex(Ismin+1:Ismax-1,Ismin,Ismin+1)=Ab1x(Ismin+1:Ismax-1,3,Ismin+1)-(Ex(Ismin+1:Ismax-1,Ismin+1,Ismin+1)-Ab1x(Ismin+1:Ismax-1,1,Ismin+1))/3;
    Ez(Ismin+1:Ismax-1,Ismin,Ismax-1)=Ab1z(Ismin+1:Ismax-1,3,Ismax-1)-(Ez(Ismin+1:Ismax-1,Ismin+1,Ismax-1)-Ab1z(Ismin+1:Ismax-1,1,Ismax-1))/3;