2-d fdtd te code.txt

使用有限时间差分算法来计算2维目标散射的例子

***********************************************************************
     2-D FDTD TE code with PML absorbing boundary conditions
***********************************************************************

     Program author: Susan C. Hagness
                     Department of Electrical and Computer Engineering
                     University of Wisconsin-Madison
                     1415 Engineering Drive
                     Madison, WI 53706-1691
                     608-265-5739
                    [e-mail]

     Date of this version:  February 2000

     This MATLAB M-file implements the finite-difference time-domain
     solution of Maxwell's curl equations over a two-dimensional
     Cartesian space lattice comprised of uniform square grid cells.

     To illustrate the algorithm, a 6-cm-diameter metal cylindrical
     scatterer in free space is modeled. The source excitation is
     a Gaussian pulse with a carrier frequency of 5 GHz.

     The grid resolution (dx = 3 mm) was chosen to provide 20 samples
     per wavelength at the center frequency of the pulse (which in turn
     provides approximately 10 samples per wavelength at the high end
     of the excitation spectrum, around 10 GHz).

     The computational domain is truncated using the perfectly matched
     layer (PML) absorbing boundary conditions.  The formulation used
     in this code is based on the original split-field Berenger PML. The
     PML regions are labeled as shown in the following diagram:

            ----------------------------------------------
           |  |                BACK PML                |  |
            ----------------------------------------------
           |L |                                       /| R|
           |E |                                (ib,jb) | I|
           |F |                                        | G|
           |T |                                        | H|
           |  |                MAIN GRID               | T|
           |P |                                        |  |
           |M |                                        | P|
           |L | (1,1)                                  | M|
           |  |/                                       | L|
            ----------------------------------------------
           |  |                FRONT PML               |  |
            ----------------------------------------------

     To execute this M-file, type "fdtd2D" at the MATLAB prompt.
     This M-file displays the FDTD-computed Ex, Ey, and Hz fields at
     every 4th time step, and records those frames in a movie matrix,
     M, which is played at the end of the simulation using the "movie"
     command.

***********************************************************************

clear

***********************************************************************
     Fundamental constants
***********************************************************************

cc=2.99792458e8;            %speed of light in free space
muz=4.0*pi*1.0e-7;          %permeability of free space
epsz=1.0/(cc*cc*muz);       %permittivity of free space

freq=5.0e+9;                %center frequency of source excitation
lambda=cc/freq;             %center wavelength of source excitation
omega=2.0*pi*freq;

***********************************************************************
     Grid parameters
***********************************************************************

ie=100;           %number of grid cells in x-direction
je=50;            %number of grid cells in y-direction

ib=ie+1;
jb=je+1;

is=15;            %location of z-directed hard source
js=je/2;          %location of z-directed hard source

dx=3.0e-3;        %space increment of square lattice
dt=dx/(2.0*cc);   %time step

nmax=300;         %total number of time steps

iebc=8;           %thickness of left and right PML region
jebc=8;           %thickness of front and back PML region
rmax=0.00001;
orderbc=2;
ibbc=iebc+1;
jbbc=jebc+1;
iefbc=ie+2*iebc;
jefbc=je+2*jebc;
ibfbc=iefbc+1;
jbfbc=jefbc+1;

***********************************************************************
     Material parameters
***********************************************************************
media=2;

eps=[1.0 1.0];
sig=[0.0 1.0e+7];
mur=[1.0 1.0];
sim=[0.0 0.0];

***********************************************************************
     Wave excitation
***********************************************************************

rtau=160.0e-12;
tau=rtau/dt;
delay=3*tau;

source=zeros(1,nmax);
for n=1:7.0*tau
  source(n)=sin(omega*(n-delay)*dt)*exp(-((n-delay)^2/tau^2));
end

***********************************************************************
     Field arrays
***********************************************************************

ex=zeros(ie,jb);           %fields in main grid
ey=zeros(ib,je);
hz=zeros(ie,je);

exbcf=zeros(iefbc,jebc);   %fields in front PML region
eybcf=zeros(ibfbc,jebc);
hzxbcf=zeros(iefbc,jebc);
hzybcf=zeros(iefbc,jebc);

exbcb=zeros(iefbc,jbbc);   %fields in back PML region
eybcb=zeros(ibfbc,jebc);
hzxbcb=zeros(iefbc,jebc);
hzybcb=zeros(iefbc,jebc);

exbcl=zeros(iebc,jb);      %fields in left PML region
eybcl=zeros(iebc,je);
hzxbcl=zeros(iebc,je);
hzybcl=zeros(iebc,je);

exbcr=zeros(iebc,jb);      %fields in right PML region
eybcr=zeros(ibbc,je);
hzxbcr=zeros(iebc,je);
hzybcr=zeros(iebc,je);

***********************************************************************
     Updating coefficients
***********************************************************************

for i=1:media
  eaf  =dt*sig(i)/(2.0*epsz*eps(i));
  ca(i)=(1.0-eaf)/(1.0+eaf);
  cb(i)=dt/epsz/eps(i)/dx/(1.0+eaf);
  haf  =dt*sim(i)/(2.0*muz*mur(i));
  da(i)=(1.0-haf)/(1.0+haf);
  db(i)=dt/muz/mur(i)/dx/(1.0+haf);
end

***********************************************************************
     Geometry specification (main grid)
***********************************************************************

     Initialize entire main grid to free space

caex(1:ie,1:jb)=ca(1);
cbex(1:ie,1:jb)=cb(1);

caey(1:ib,1:je)=ca(1);
cbey(1:ib,1:je)=cb(1);

dahz(1:ie,1:je)=da(1);
dbhz(1:ie,1:je)=db(1);

     Add metal cylinder

diam=20;          % diameter of cylinder: 6 cm
rad=diam/2.0;     % radius of cylinder: 3 cm
icenter=4*ie/5;   % i-coordinate of cylinder's center
jcenter=je/2;     % j-coordinate of cylinder's center

for i=1:ie
for j=1:je
  dist2=(i+0.5-icenter)^2 + (j-jcenter)^2;
  if dist2 <= rad^2
     caex(i,j)=ca(2);
     cbex(i,j)=cb(2);
  end
  dist2=(i-icenter)^2 + (j+0.5-jcenter)^2;
  if dist2 <= rad^2
     caey(i,j)=ca(2);
     cbey(i,j)=cb(2);
  end
end
end

***********************************************************************
     Fill the PML regions
***********************************************************************

delbc=iebc*dx;
sigmam=-log(rmax/100.0)*epsz*cc*(orderbc+1)/(2*delbc);
bcfactor=eps(1)*sigmam/(dx*(delbc^orderbc)*(orderbc+1));

     FRONT region

caexbcf(1:iefbc,1)=1.0;
cbexbcf(1:iefbc,1)=0.0;
for j=2:jebc
  y1=(jebc-j+1.5)*dx;
  y2=(jebc-j+0.5)*dx;
  sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
  ca1=exp(-sigmay*dt/(epsz*eps(1)));
  cb1=(1.0-ca1)/(sigmay*dx);
  caexbcf(1:iefbc,j)=ca1;
  cbexbcf(1:iefbc,j)=cb1;
end
sigmay = bcfactor*(0.5*dx)^(orderbc+1);
ca1=exp(-sigmay*dt/(epsz*eps(1)));
cb1=(1-ca1)/(sigmay*dx);
caex(1:ie,1)=ca1;
cbex(1:ie,1)=cb1;
caexbcl(1:iebc,1)=ca1;
cbexbcl(1:iebc,1)=cb1;
caexbcr(1:iebc,1)=ca1;
cbexbcr(1:iebc,1)=cb1;

for j=1:jebc
  y1=(jebc-j+1)*dx;
  y2=(jebc-j)*dx;
  sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
  sigmays=sigmay*(muz/(epsz*eps(1)));
  da1=exp(-sigmays*dt/muz);
  db1=(1-da1)/(sigmays*dx);
  dahzybcf(1:iefbc,j)=da1;
  dbhzybcf(1:iefbc,j)=db1;
  caeybcf(1:ibfbc,j)=ca(1);
  cbeybcf(1:ibfbc,j)=cb(1);
  dahzxbcf(1:iefbc,j)=da(1);
  dbhzxbcf(1:iefbc,j)=db(1);
end

     BACK region

caexbcb(1:iefbc,jbbc)=1.0;
cbexbcb(1:iefbc,jbbc)=0.0;
for j=2:jebc
  y1=(j-0.5)*dx;
  y2=(j-1.5)*dx;
  sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
  ca1=exp(-sigmay*dt/(epsz*eps(1)));
  cb1=(1-ca1)/(sigmay*dx);
  caexbcb(1:iefbc,j)=ca1;
  cbexbcb(1:iefbc,j)=cb1;
end
sigmay = bcfactor*(0.5*dx)^(orderbc+1);
ca1=exp(-sigmay*dt/(epsz*eps(1)));
cb1=(1-ca1)/(sigmay*dx);
caex(1:ie,jb)=ca1;
cbex(1:ie,jb)=cb1;
caexbcl(1:iebc,jb)=ca1;
cbexbcl(1:iebc,jb)=cb1;
caexbcr(1:iebc,jb)=ca1;
cbexbcr(1:iebc,jb)=cb1;

for j=1:jebc
  y1=j*dx;
  y2=(j-1)*dx;
  sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
  sigmays=sigmay*(muz/(epsz*eps(1)));
  da1=exp(-sigmays*dt/muz);
  db1=(1-da1)/(sigmays*dx);
  dahzybcb(1:iefbc,j)=da1;
  dbhzybcb(1:iefbc,j)=db1;
  caeybcb(1:ibfbc,j)=ca(1);
  cbeybcb(1:ibfbc,j)=cb(1);
  dahzxbcb(1:iefbc,j)=da(1);
  dbhzxbcb(1:iefbc,j)=db(1);
end

     LEFT region

caeybcl(1,1:je)=1.0;
cbeybcl(1,1:je)=0.0;
for i=2:iebc
  x1=(iebc-i+1.5)*dx;
  x2=(iebc-i+0.5)*dx;
  sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
  ca1=exp(-sigmax*dt/(epsz*eps(1)));
  cb1=(1-ca1)/(sigmax*dx);
  caeybcl(i,1:je)=ca1;
  cbeybcl(i,1:je)=cb1;
  caeybcf(i,1:jebc)=ca1;
  cbeybcf(i,1:jebc)=cb1;
  caeybcb(i,1:jebc)=ca1;
  cbeybcb(i,1:jebc)=cb1;
end
sigmax=bcfactor*(0.5*dx)^(orderbc+1);
ca1=exp(-sigmax*dt/(epsz*eps(1)));
cb1=(1-ca1)/(sigmax*dx);
caey(1,1:je)=ca1;
cbey(1,1:je)=cb1;
caeybcf(iebc+1,1:jebc)=ca1;
cbeybcf(iebc+1,1:jebc)=cb1;
caeybcb(iebc+1,1:jebc)=ca1;
cbeybcb(iebc+1,1:jebc)=cb1;

for i=1:iebc
  x1=(iebc-i+1)*dx;
  x2=(iebc-i)*dx;
  sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
  sigmaxs=sigmax*(muz/(epsz*eps(1)));
  da1=exp(-sigmaxs*dt/muz);
  db1=(1-da1)/(sigmaxs*dx);
  dahzxbcl(i,1:je)=da1;
  dbhzxbcl(i,1:je)=db1;