xg.f90

fortran 编写的三维fdtd计算程序

program xg
!变量声明
integer Ntime,wge,obj(101,101,101),Nmedia,Ninwave,wgec,Nt,Ttz,i,j,k,l; 
integer Izmin,Izmax,Ismin,Ismax,Fzmin,Fzmax,TEM_flag;
double precision freq,Tin,media(2,4),Vin,t0,Phi,Tty,Tyu,WL,miu0,epslon0,vc,wavel,Z0,delta,dt,FE,FH,EBin(4),  A(300)  ,pi;
double precision Ex(101,101,101),Ey(101,101,101),Ez(101,101,101),Hx(101,101,101),Hy(101,101,101),Hz(101,101,101);
double precision Ein(2*300), Hin(2*300),   CA(2,2,2,2),CB(2,2,2,2),CP(2,2),CQ(2,2);
double precision Ab1x(101,4,101),Ab1z(101,4,101),Ab2x(101,4,101),Ab2z(101,4,101);
double precision Ab3x(101,101,4),Ab3y(101,101,4),Ab4x(101,101,4),Ab4y(101,101,4);
double precision Ab5z(4,101,101),Ab5y(4,101,101),Ab6z(4,101,101),Ab6y(4,101,101);
double precision Ein34(101,101,1),Hin34(101,101,1),Ein1(101,1,101),Ein2(101,1,101),Hin1(101,1,101),Hin2(101,1,101),Ein5(1,101,101),Ein6(1,101,101),Hin5(1,101,101),Hin6(1,101,101);
!变量声明结束

Ntime=300;

wge=101;
obj=1;                 !存放空间媒质编号，包括散射体，迭代次数
freq=1.0e9; Tin=0.15/freq;   !频率上限,入射波时间常数
Nmedia=2;
media(1,:)=[1.0,1.0,0.0,0.0];media(2,:)=[1.0,1.0,3.72E7,0.0]; !可存Nmedia种介质的[epslon_r,miu_r,sigma,sigma_m]
Vin=1;t0=0.6e-9;
Phi=0;WL=40;            !平面波入射角和数值波长，wavel有40个网格
Ninwave=2*wge+20;       !一阶FDTD长度
!---------不变参数-------------
pi=3.14159265358979;
miu0=4e-7*pi;epslon0=8.854187818e-12;            !真空电磁参数
vc=299792458;wavel=vc/freq;                      !光速和波长设置
Z0=sqrt(miu0/epslon0);                           !波阻抗
wgec=real(wge)/2+0.5;                            !中心网格坐标
!------------------------------
!------导入模型到odj离散空间----------
!readmodel;
![xs,ys,zs]=size(mmd);
!obj(wge-floor(xs/2):wge-floor(xs/2)+xs-1,wge-floor(ys/2):wge-floor(ys/2)+ys-1,wge-floor(zs/2):wge-floor(zs/2)+zs-1)=mmd+1;
obj(wgec-10:wgec+10,wgec-10:wgec+10,wgec-10:wgec+10)=obj(wgec-10:wgec+10,wgec-10:wgec+10,wgec-10:wgec+10)+1;
!----------导入模型结束---------------
Izmin=25;Izmax=wge-25;!总场区边界
Ismin=5;Ismax=wge-5;!散射场区边界
Fzmin=9;Fzmax=wge-9;!远场外推边界
TEM_flag=1;        !TE（＝1） or TM（＝2）

Ez=0;Ex=Ez;Ey=Ez;Hx=Ez;Hy=Ez;Hz=Ez;!存放电磁场的空间
Ein=0;Hin=Ein;!存放1-D FDTD入射电磁场的空间
CA=0;CB=CA;!电场迭代系数
CP=0;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=0; Ab1z=Ab1x;      !1面x,z分量,4层
Ab2x=0; Ab2z=Ab2x;      !2面x,z分量,4层
Ab3x=0; Ab3y=Ab3x;      !3面x,y分量,4层
Ab4x=0; Ab4y=Ab4x;      !4面x,y分量,4层
Ab5z=0; Ab5y=Ab5z;      !5面z,y分量,4层
Ab6z=0; Ab6y=Ab6z;      !6面z,y分量,4层
!---吸收边界暂存空间 end----

!计算迭代系数
do i=1,Nmedia
    do j=1,Nmedia
        do k=1,Nmedia
            do 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 do
        end do
        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 do
end do
FE=dt/epslon0/delta;  FH=dt/miu0/delta;           !入射波1-D FDTD迭代系数
EBin=0;                                           !1-D FDTD吸收边界变量

!申请空间存放表面入射切向电场分量和磁场分量,算出个表面的E，除以Z0就是H，乘以角度就是分量
Ein34=0;  Hin34=Ein34;!3,4面，由于入射波在Z向均匀，所以3，4面切向场相同共用一个
Ein1=0; Ein2=Ein1; Hin1=Ein1; Hin2=Ein1;
Ein5=0; Ein6=Ein5; Hin5=Ein5; Hin6=Ein5;
!====================================prepare work end=====================================================
!
open (1,file='a1.txt')
!=========================================Main loop=======================================================
do Nt=1,Ntime
    A(Nt)=Ez(wgec,Izmin-5,wgec);
    !---------creat add_in wave using 1-D FDTD----------------
    do i=1,Ninwave-1
        Hin(i)=Hin(i)-FH*(Ein(i+1)-Ein(i));
    end do
    do i=2,Ninwave-1
        Ein(i)=Ein(i)-FE*(Hin(i)-Hin(i-1));
    end do
    Ein(20)=Vin*exp(-((Nt*dt-t0)/Tin)**2); !sin(10*Nt*dt/Tin);!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-----
    !                 2
    !                /
    !          _________
    !         /   4   y/|
    !        /___x____/ |  
    !        |        |6|
    !      5-|z   1   | /        
    !        |________|/
    !            |
    !            3
    !申请空间存放表面 切向电场分量和磁场分量,算出个表面的E&H，乘以角度就是分量
    !Ein34=zeros(wge,wge,1);!3,4面，由于入射波在Z向均匀，所以3，4面切向场相同共用一个
    !Ein1=zeros(wge,1,wge);Ein2=Ein1;
    !Ein5=zeros(1,wge,wge);Ein6=Ein5;
    do i=Izmin,Izmax
        do j=Izmin,Izmax
            !--------3 4 面---------
            Tty=i*sin(Phi)+j*cos(Phi);  Ttz=int(Tty);  Tyu=Tty-Ttz;          
            Ein34(i,j,1)=Ein(Ttz)+Tyu*(Ein(Ttz+1)-Ein(Ttz));
            Hin34(i,j,1)=Hin(Ttz)+Tyu*(Hin(Ttz+1)-Hin(Ttz));
            !---------1 面----------
            Tty=i*sin(Phi)+Izmin*cos(Phi);  Ttz=int(Tty);  Tyu=Tty-Ttz;          
            Ein1(i,1,j)=Ein(Ttz)+Tyu*(Ein(Ttz+1)-Ein(Ttz));
            Tty=i*sin(Phi)+(Izmin-0.5)*cos(Phi);  Ttz=int(Tty-0.5)+1;  Tyu=Ttz+0.5-Tty;
            Hin1(i,1,j)=Hin(Ttz)+Tyu*(Hin(Ttz-1)-Hin(Ttz));
            !---------2 面----------
            Tty=i*sin(Phi)+Izmax*cos(Phi);  Ttz=int(Tty);  Tyu=Tty-Ttz;          
            Ein2(i,1,j)=Ein(Ttz)+Tyu*(Ein(Ttz+1)-Ein(Ttz));
            Tty=i*sin(Phi)+(Izmax+0.5)*cos(Phi);  Ttz=int(Tty-0.5)+1;  Tyu=Ttz+0.5-Tty;
            Hin2(i,1,j)=Hin(Ttz)+Tyu*(Hin(Ttz-1)-Hin(Ttz));
            !---------5 面----------
            Tty=Izmin*sin(Phi)+j*cos(Phi);  Ttz=int(Tty);  Tyu=Tty-Ttz;          
            Ein5(1,j,i)=Ein(Ttz)+Tyu*(Ein(Ttz+1)-Ein(Ttz));
            Tty=(Izmin-0.5)*sin(Phi)+j*cos(Phi);  Ttz=int(Tty-0.5)+1;  Tyu=Ttz+0.5-Tty;
            Hin5(1,j,i)=Hin(Ttz)+Tyu*(Hin(Ttz-1)-Hin(Ttz));
            !---------6 面----------
            Tty=Izmax*sin(Phi)+j*cos(Phi);  Ttz=int(Tty);  Tyu=Tty-Ttz;          
            Ein6(1,j,i)=Ein(Ttz)+Tyu*(Ein(Ttz+1)-Ein(Ttz));
            Tty=(Izmax+0.5)*sin(Phi)+j*cos(Phi);  Ttz=int(Tty-0.5)+1;  Tyu=Ttz+0.5-Tty;
            Hin6(1,j,i)=Hin(Ttz)+Tyu*(Hin(Ttz-1)-Hin(Ttz));
        end do
    end do
    !------end of caculate incident wave E and H feild on each connect plane------
    !
    !--------------------------------------------------start E FDTD----------------------------------------------------
    do i=Ismin+1,Ismax-1
        do j=Ismin+1,Ismax-1
            do 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 do
        end do
    end do
    !-----连接边界面引入Ez分量-----
    Ez(Izmin+1:Izmax-1,Izmin,Izmin+1:Izmax-1)=Ez(Izmin+1:Izmax-1,Izmin,Izmin+1:Izmax-1)+FE*cos(Phi)*Hin1(Izmin+1:Izmax-1,1,Izmin+1:Izmax-1);!面1
    Ez(Izmin+1:Izmax-1,Izmax,Izmin+1:Izmax-1)=Ez(Izmin+1:Izmax-1,Izmax,Izmin+1:Izmax-1)-FE*cos(Phi)*Hin2(Izmin+1:Izmax-1,1,Izmin+1:Izmax-1);!面2
    !----面3----
    Ex(Izmin+1:Izmax-1,Izmin+1:Izmax-1,Izmin)=Ex(Izmin+1:Izmax-1,Izmin+1:Izmax-1,Izmin)+FE*sin(Phi)*Hin34(Izmin+1:Izmax-1,Izmin+1:Izmax-1,1);
    Ey(Izmin+1:Izmax-1,Izmin+1:Izmax-1,Izmin)=Ey(Izmin+1:Izmax-1,Izmin+1:Izmax-1,Izmin)-FE*cos(Phi)*Hin34(Izmin+1:Izmax-1,Izmin+1:Izmax-1,1);
    !--面3 end--
    !----面4----
    Ex(Izmin+1:Izmax-1,Izmin+1:Izmax-1,Izmax)=Ex(Izmin+1:Izmax-1,Izmin+1:Izmax-1,Izmax)-FE*sin(Phi)*Hin34(Izmin+1:Izmax-1,Izmin+1:Izmax-1,1);
    Ey(Izmin+1:Izmax-1,Izmin+1:Izmax-1,Izmax)=Ey(Izmin+1:Izmax-1,Izmin+1:Izmax-1,Izmax)+FE*cos(Phi)*Hin34(Izmin+1:Izmax-1,Izmin+1:Izmax-1,1);
    !--面4 end--
    Ez(Izmin,Izmin+1:Izmax-1,Izmin+1:Izmax-1)=Ez(Izmin,Izmin+1:Izmax-1,Izmin+1:Izmax-1)-FE*sin(Phi)*Hin5(1,Izmin+1:Izmax-1,Izmin+1:Izmax-1);!面5
    Ez(Izmax,Izmin+1:Izmax-1,Izmin+1:Izmax-1)=Ez(Izmax,Izmin+1:Izmax-1,Izmin+1:Izmax-1)+FE*sin(Phi)*Hin6(1,Izmin+1:Izmax-1,Izmin+1:Izmax-1);!面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)*Hin5(1,Izmin,Izmin+1:Izmax-1)+FE*cos(Phi)*Hin1(Izmin,1,Izmin+1:Izmax-1);!1,5面交线
    Ez(Izmax,Izmin,Izmin+1:Izmax-1)=Ez(Izmax,Izmin,Izmin+1:Izmax-1)+FE*sin(Phi)*Hin6(1,Izmin,Izmin+1:Izmax-1)+FE*cos(Phi)*Hin1(Izmax,1,Izmin+1:Izmax-1);!1,6面交线
    Ez(Izmax,Izmax,Izmin+1:Izmax-1)=Ez(Izmax,Izmax,Izmin+1:Izmax-1)+FE*sin(Phi)*Hin6(1,Izmax,Izmin+1:Izmax-1)-FE*cos(Phi)*Hin2(Izmax,1,Izmin+1:Izmax-1);!2,6面交线
    Ez(Izmin,Izmax,Izmin+1:Izmax-1)=Ez(Izmin,Izmax,Izmin+1:Izmax-1)-FE*sin(Phi)*Hin5(1,Izmax,Izmin+1:Izmax-1)-FE*cos(Phi)*Hin2(Izmin,1,Izmin+1:Izmax-1);!2,5面交线
    !--//X Ex分量----
    Ex(Izmin+1:Izmax-1,Izmin,Izmin)=Ex(Izmin+1:Izmax-1,Izmin,Izmin)+FE*sin(Phi)*Hin34(Izmin+1:Izmax-1,Izmin,1);       !1,3面交线
    Ex(Izmin+1:Izmax-1,Izmax,Izmin)=Ex(Izmin+1:Izmax-1,Izmax,Izmin)+FE*sin(Phi)*Hin34(Izmin+1:Izmax-1,Izmax,1);       !2,3面交线
    Ex(Izmin+1:Izmax-1,Izmax,Izmax)=Ex(Izmin+1:Izmax-1,Izmax,Izmax)-FE*sin(Phi)*Hin34(Izmin+1:Izmax-1,Izmax,1);       !2,4面交线
    Ex(Izmin+1:Izmax-1,Izmin,Izmax)=Ex(Izmin+1:Izmax-1,Izmin,Izmax)-FE*sin(Phi)*Hin34(Izmin+1:Izmax-1,Izmin,1);       !1,4面交线
    !--//Y Ey分量----
    Ey(Izmin,Izmin+1:Izmax-1,Izmin)=Ey(Izmin,Izmin+1:Izmax-1,Izmin)-FE*cos(Phi)*Hin34(Izmin,Izmin+1:Izmax-1,1);       !5,3面交线
    Ey(Izmax,Izmin+1:Izmax-1,Izmin)=Ey(Izmax,Izmin+1:Izmax-1,Izmin)+FE*cos(Phi)*Hin34(Izmax,Izmin+1:Izmax-1,1);       !6,3面交线
    Ey(Izmax,Izmin+1:Izmax-1,Izmax)=Ey(Izmax,Izmin+1:Izmax-1,Izmax)+FE*cos(Phi)*Hin34(Izmax,Izmin+1:Izmax-1,1);       !6,4面交线
    Ey(Izmin,Izmin+1:Izmax-1,Izmax)=Ey(Izmin,Izmin+1:Izmax-1,Izmax)-FE*cos(Phi)*Hin34(Izmin,Izmin+1:Izmax-1,1);       !5,4面交线
    !-------------------------------------------------吸收边界-------------------------------------------------------
    !-------------------------------------6个面-----------------------------------------
    !---------------面1(y=Ismin,吸收Ez、Ex)----------------------
    do i=Ismin+2,Ismax-2
        do j=Ismin+2,Ismax-2