3D FDTD（matlab）源码

离线gwzhao

方恨少

UID ：17098

注册：2008-08-24
登录：2019-01-09
发帖：1374
等级：荣誉管理员

0楼发表于: 2008-10-23 14:40:30

— 本帖被 tensor 从资料库移动到本区(2009-10-28) —

[post]%***********************************************************************
%     1-D FDTD code with simple radiation 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
%                     hagness@engr.wisc.edu
%
%     Date of this version:  February 2000
%
%     This MATLAB M-file implements the finite-difference time-domain
%     solution of Maxwell's curl equations over a one-dimensional space
%     lattice comprised of uniform grid cells.
%
%     To illustrate the algorithm, a sinusoidal wave (1GHz) propagating
%     in a nonpermeable lossy medium (epsr=1.0, sigma=5.0e-3 S/m) is
%     modeled.  The simplified finite difference system for nonpermeable
%     media (discussed in Section 3.6.6 of the text) is implemented.
%
%     The grid resolution (dx = 1.5 cm) is chosen to provide 20
%     samples per wavelength.  The Courant factor S=c*dt/dx is set to
%     the stability limit: S=1.  In 1-D, this is the "magic time step."
%
%     The computational domain is truncated using the simplest radiation
%     boundary condition for wave propagation in free space:
%
%                      Ez(imax,n+1) = Ez(imax-1,n)
%
%     To execute this M-file, type "fdtd1D" at the MATLAB prompt.
%     This M-file displays the FDTD-computed Ez and Hy fields at every
%     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=1.0e+9;                %frequency of source excitation
lambda=cc/freq;             %wavelength of source excitation
omega=2.0*pi*freq;

%***********************************************************************
%     Grid parameters
%***********************************************************************

ie=200;                     %number of grid cells in x-direction

ib=ie+1;

dx=lambda/20.0;             %space increment of 1-D lattice
dt=dx/cc;                   %time step
omegadt=omega*dt;

nmax=round(12.0e-9/dt);     %total number of time steps

%***********************************************************************
%     Material parameters
%***********************************************************************

eps=1.0;
sig=5.0e-3;

%***********************************************************************
%     Updating coefficients for space region with nonpermeable media
%***********************************************************************

scfact=dt/muz/dx;

ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps));
cb=scfact*(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps));

%***********************************************************************
%     Field arrays
%***********************************************************************

ez(1:ib)=0.0;
hy(1:ie)=0.0;

%***********************************************************************
%     Movie initialization
%***********************************************************************

x=linspace(dx,ie*dx,ie);

subplot(2,1,1),plot(x,ez(1:ie)/scfact,'r'),axis([0 3 -1 1]);
ylabel('EZ');

subplot(2,1,2),plot(x,hy,'b'),axis([0 3 -3.0e-3 3.0e-3]);
xlabel('x (meters)');ylabel('HY');

rect=get(gcf,'Position');
rect(1:2)=[0 0];

M=moviein(nmax/2,gcf,rect);

%***********************************************************************
%     BEGIN TIME-STEPPING LOOP
%***********************************************************************

for n=1:nmax

%***********************************************************************
%     Update electric fields
%***********************************************************************

ez(1)=scfact*sin(omegadt*n);

rbc=ez(ie);
ez(2:ie)=ca*ez(2:ie)+cb*(hy(2:ie)-hy(1:ie-1));
ez(ib)=rbc;

%***********************************************************************
%     Update magnetic fields
%***********************************************************************

hy(1:ie)=hy(1:ie)+ez(2:ib)-ez(1:ie);

%***********************************************************************
%     Visualize fields
%*********** ..

未注册仅能浏览部分内容，查看全部内容及附件请先登录或注册

共3条评分，rf币+11 ，技术分+1

chenjiabao1989	rf币 +10	十分期待您分享下一贴！	2014-09-29
casey	rf币 +1	有自己的观点	2008-10-23
casey	技术分 +1	当时加错了，补上	2008-10-23

逆流而上

离线gwzhao

方恨少

UID ：17098

注册：2008-08-24
登录：2019-01-09
发帖：1374
等级：荣誉管理员

1楼发表于: 2008-10-23 14:40:53

%***********************************************************************
%     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
%                     hagness@engr.wisc.edu
%
%     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;
  dahzxbcf(i,1:jebc)=da1;
  dbhzxbcf(i,1:jebc)=db1;
  dahzxbcb(i,1:jebc)=da1;
  dbhzxbcb(i,1:jebc)=db1;
  caexbcl(i,2:je)=ca(1);
  cbexbcl(i,2:je)=cb(1);
  dahzybcl(i,1:je)=da(1);
  dbhzybcl(i,1:je)=db(1);
end

%     RIGHT region

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

for i=1:iebc
  x1=i*dx;
  x2=(i-1)*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);
  dahzxbcr(i,1:je) = da1;
  dbhzxbcr(i,1:je) = db1;
  dahzxbcf(i+ie+iebc,1:jebc)=da1;
  dbhzxbcf(i+ie+iebc,1:jebc)=db1;
  dahzxbcb(i+ie+iebc,1:jebc)=da1;
  dbhzxbcb(i+ie+iebc,1:jebc)=db1;
  caexbcr(i,2:je)=ca(1);
  cbexbcr(i,2:je)=cb(1);
  dahzybcr(i,1:je)=da(1);
  dbhzybcr(i,1:je)=db(1);
end

%***********************************************************************
%     Movie initialization
%***********************************************************************

subplot(3,1,1),pcolor(ex');
shading flat;
caxis([-80.0 80.0]);
axis([1 ie 1 jb]);
colorbar;
axis image;
axis off;
title(['Ex at time step = 0']);

subplot(3,1,2),pcolor(ey');
shading flat;
caxis([-80.0 80.0]);
axis([1 ib 1 je]);
colorbar;
axis image;
axis off;
title(['Ey at time step = 0']);

subplot(3,1,3),pcolor(hz');
shading flat;
caxis([-0.2 0.2]);
axis([1 ie 1 je]);
colorbar;
axis image;
axis off;
title(['Hz at time step = 0']);

rect=get(gcf,'Position');
rect(1:2)=[0 0];

M=moviein(nmax/4,gcf,rect);

%***********************************************************************
%     BEGIN TIME-STEPPING LOOP
%***********************************************************************

for n=1:nmax

%***********************************************************************
%     Update electric fields (EX and EY) in main grid
%***********************************************************************

ex(:,2:je)=caex(:,2:je).*ex(:,2:je)+...
           cbex(:,2:je).*(hz(:,2:je)-hz(:,1:je-1));

ey(2:ie,:)=caey(2:ie,:).*ey(2:ie,:)+...
           cbey(2:ie,:).*(hz(1:ie-1,:)-hz(2:ie,:));

%***********************************************************************
%     Update EX in PML regions
%***********************************************************************

%     FRONT

exbcf(:,2:jebc)=caexbcf(:,2:jebc).*exbcf(:,2:jebc)-...
  cbexbcf(:,2:jebc).*(hzxbcf(:,1:jebc-1)+hzybcf(:,1:jebc-1)-...
                      hzxbcf(:,2:jebc)-hzybcf(:,2:jebc));
ex(1:ie,1)=caex(1:ie,1).*ex(1:ie,1)-...
  cbex(1:ie,1).*(hzxbcf(ibbc:iebc+ie,jebc)+...
                hzybcf(ibbc:iebc+ie,jebc)-hz(1:ie,1));

%     BACK

exbcb(:,2:jebc-1)=caexbcb(:,2:jebc-1).*exbcb(:,2:jebc-1)-...
  cbexbcb(:,2:jebc-1).*(hzxbcb(:,1:jebc-2)+hzybcb(:,1:jebc-2)-...
                        hzxbcb(:,2:jebc-1)-hzybcb(:,2:jebc-1));
ex(1:ie,jb)=caex(1:ie,jb).*ex(1:ie,jb)-...
  cbex(1:ie,jb).*(hz(1:ie,jb-1)-hzxbcb(ibbc:iebc+ie,1)-...
                 hzybcb(ibbc:iebc+ie,1));

%     LEFT

exbcl(:,2:je)=caexbcl(:,2:je).*exbcl(:,2:je)-...
  cbexbcl(:,2:je).*(hzxbcl(:,1:je-1)+hzybcl(:,1:je-1)-...
                    hzxbcl(:,2:je)-hzybcl(:,2:je));
exbcl(:,1)=caexbcl(:,1).*exbcl(:,1)-...
  cbexbcl(:,1).*(hzxbcf(1:iebc,jebc)+hzybcf(1:iebc,jebc)-...
                 hzxbcl(:,1)-hzybcl(:,1));
exbcl(:,jb)=caexbcl(:,jb).*exbcl(:,jb)-...
  cbexbcl(:,jb).*(hzxbcl(:,je)+hzybcl(:,je)-...
                  hzxbcb(1:iebc,1)-hzybcb(1:iebc,1));

%     RIGHT

exbcr(:,2:je)=caexbcr(:,2:je).*exbcr(:,2:je)-...
  cbexbcr(:,2:je).*(hzxbcr(:,1:je-1)+hzybcr(:,1:je-1)-...
                    hzxbcr(:,2:je)-hzybcr(:,2:je));
exbcr(:,1)=caexbcr(:,1).*exbcr(:,1)-...
  cbexbcr(:,1).*(hzxbcf(1+iebc+ie:iefbc,jebc)+...
                 hzybcf(1+iebc+ie:iefbc,jebc)-...
                 hzxbcr(:,1)-hzybcr(:,1));
exbcr(:,jb)=caexbcr(:,jb).*exbcr(:,jb)-...
  cbexbcr(:,jb).*(hzxbcr(:,je)+hzybcr(:,je)-...
                  hzxbcb(1+iebc+ie:iefbc,1)-...
                  hzybcb(1+iebc+ie:iefbc,1));

%***********************************************************************
%     Update EY in PML regions
%***********************************************************************

%     FRONT

eybcf(2:iefbc,:)=caeybcf(2:iefbc,:).*eybcf(2:iefbc,:)-...
  cbeybcf(2:iefbc,:).*(hzxbcf(2:iefbc,:)+hzybcf(2:iefbc,:)-...
                       hzxbcf(1:iefbc-1,:)-hzybcf(1:iefbc-1,:));

%     BACK

eybcb(2:iefbc,:)=caeybcb(2:iefbc,:).*eybcb(2:iefbc,:)-...
  cbeybcb(2:iefbc,:).*(hzxbcb(2:iefbc,:)+hzybcb(2:iefbc,:)-...
                       hzxbcb(1:iefbc-1,:)-hzybcb(1:iefbc-1,:));

%     LEFT

eybcl(2:iebc,:)=caeybcl(2:iebc,:).*eybcl(2:iebc,:)-...
  cbeybcl(2:iebc,:).*(hzxbcl(2:iebc,:)+hzybcl(2:iebc,:)-...
                      hzxbcl(1:iebc-1,:)-hzybcl(1:iebc-1,:));
ey(1,:)=caey(1,:).*ey(1,:)-...
  cbey(1,:).*(hz(1,:)-hzxbcl(iebc,:)-hzybcl(iebc,:));

%     RIGHT

eybcr(2:iebc,:)=caeybcr(2:iebc,:).*eybcr(2:iebc,:)-...
  cbeybcr(2:iebc,:).*(hzxbcr(2:iebc,:)+hzybcr(2:iebc,:)-...
                      hzxbcr(1:iebc-1,:)-hzybcr(1:iebc-1,:));
ey(ib,:)=caey(ib,:).*ey(ib,:)-...
  cbey(ib,:).*(hzxbcr(1,:)+hzybcr(1,:)- hz(ie,:));

%***********************************************************************
%     Update magnetic fields (HZ) in main grid
%***********************************************************************

hz(1:ie,1:je)=dahz(1:ie,1:je).*hz(1:ie,1:je)+...
              dbhz(1:ie,1:je).*(ex(1:ie,2:jb)-ex(1:ie,1:je)+...
                                ey(1:ie,1:je)-ey(2:ib,1:je));

hz(is,js)=source(n);

%***********************************************************************
%     Update HZX in PML regions
%***********************************************************************

%     FRONT

hzxbcf(1:iefbc,:)=dahzxbcf(1:iefbc,:).*hzxbcf(1:iefbc,:)-...
  dbhzxbcf(1:iefbc,:).*(eybcf(2:ibfbc,:)-eybcf(1:iefbc,:));

%     BACK

hzxbcb(1:iefbc,:)=dahzxbcb(1:iefbc,:).*hzxbcb(1:iefbc,:)-...
  dbhzxbcb(1:iefbc,:).*(eybcb(2:ibfbc,:)-eybcb(1:iefbc,:));

%     LEFT

hzxbcl(1:iebc-1,:)=dahzxbcl(1:iebc-1,:).*hzxbcl(1:iebc-1,:)-...
  dbhzxbcl(1:iebc-1,:).*(eybcl(2:iebc,:)-eybcl(1:iebc-1,:));
hzxbcl(iebc,:)=dahzxbcl(iebc,:).*hzxbcl(iebc,:)-...
  dbhzxbcl(iebc,:).*(ey(1,:)-eybcl(iebc,:));

%     RIGHT

hzxbcr(2:iebc,:)=dahzxbcr(2:iebc,:).*hzxbcr(2:iebc,:)-...
  dbhzxbcr(2:iebc,:).*(eybcr(3:ibbc,:)-eybcr(2:iebc,:));
hzxbcr(1,:)=dahzxbcr(1,:).*hzxbcr(1,:)-...
  dbhzxbcr(1,:).*(eybcr(2,:)-ey(ib,:));

%***********************************************************************
%     Update HZY in PML regions
%***********************************************************************

%     FRONT

hzybcf(:,1:jebc-1)=dahzybcf(:,1:jebc-1).*hzybcf(:,1:jebc-1)-...
  dbhzybcf(:,1:jebc-1).*(exbcf(:,1:jebc-1)-exbcf(:,2:jebc));
hzybcf(1:iebc,jebc)=dahzybcf(1:iebc,jebc).*hzybcf(1:iebc,jebc)-...
  dbhzybcf(1:iebc,jebc).*(exbcf(1:iebc,jebc)-exbcl(1:iebc,1));
hzybcf(iebc+1:iebc+ie,jebc)=...
  dahzybcf(iebc+1:iebc+ie,jebc).*hzybcf(iebc+1:iebc+ie,jebc)-...
  dbhzybcf(iebc+1:iebc+ie,jebc).*(exbcf(iebc+1:iebc+ie,jebc)-...
                                  ex(1:ie,1));
hzybcf(iebc+ie+1:iefbc,jebc)=...
  dahzybcf(iebc+ie+1:iefbc,jebc).*hzybcf(iebc+ie+1:iefbc,jebc)-...
  dbhzybcf(iebc+ie+1:iefbc,jebc).*(exbcf(iebc+ie+1:iefbc,jebc)-...
                                   exbcr(1:iebc,1));

%     BACK

hzybcb(1:iefbc,2:jebc)=dahzybcb(1:iefbc,2:jebc).*hzybcb(1:iefbc,2:jebc)-...
  dbhzybcb(1:iefbc,2:jebc).*(exbcb(1:iefbc,2:jebc)-exbcb(1:iefbc,3:jbbc));
hzybcb(1:iebc,1)=dahzybcb(1:iebc,1).*hzybcb(1:iebc,1)-...
  dbhzybcb(1:iebc,1).*(exbcl(1:iebc,jb)-exbcb(1:iebc,2));
hzybcb(iebc+1:iebc+ie,1)=...
  dahzybcb(iebc+1:iebc+ie,1).*hzybcb(iebc+1:iebc+ie,1)-...
  dbhzybcb(iebc+1:iebc+ie,1).*(ex(1:ie,jb)-exbcb(iebc+1:iebc+ie,2));
hzybcb(iebc+ie+1:iefbc,1)=...
  dahzybcb(iebc+ie+1:iefbc,1).*hzybcb(iebc+ie+1:iefbc,1)-...
  dbhzybcb(iebc+ie+1:iefbc,1).*(exbcr(1:iebc,jb)-...
                                exbcb(iebc+ie+1:iefbc,2));

%     LEFT

hzybcl(:,1:je)=dahzybcl(:,1:je).*hzybcl(:,1:je)-...
  dbhzybcl(:,1:je).*(exbcl(:,1:je)-exbcl(:,2:jb));

%     RIGHT

hzybcr(:,1:je)=dahzybcr(:,1:je).*hzybcr(:,1:je)-...
  dbhzybcr(:,1:je).*(exbcr(:,1:je)-exbcr(:,2:jb));

%***********************************************************************
%     Visualize fields
%***********************************************************************

if mod(n,4)==0;

timestep=int2str(n);

subplot(3,1,1),pcolor(ex');
shading flat;
caxis([-80.0 80.0]);
axis([1 ie 1 jb]);
colorbar;
axis image;
axis off;
title(['Ex at time step = ',timestep]);

subplot(3,1,2),pcolor(ey');
shading flat;
caxis([-80.0 80.0]);
axis([1 ib 1 je]);
colorbar;
axis image;
axis off;
title(['Ey at time step = ',timestep]);

subplot(3,1,3),pcolor(hz');
shading flat;
caxis([-0.2 0.2]);
axis([1 ie 1 je]);
colorbar;
axis image;
axis off;
title(['Hz at time step = ',timestep]);

nn=n/4;
M(:,nn)=getframe(gcf,rect);

end;

%***********************************************************************
%     END TIME-STEPPING LOOP
%***********************************************************************

end

movie(gcf,M,0,10,rect);

共1条评分

cem-uestc

rf币 +5

优秀资料+RF币

2008-10-23

逆流而上

离线gwzhao

方恨少

UID ：17098

注册：2008-08-24
登录：2019-01-09
发帖：1374
等级：荣誉管理员

2楼发表于: 2008-10-23 14:41:19

%***********************************************************************
%     3-D FDTD code with PEC boundaries
%***********************************************************************
%
%     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
%                     hagness@engr.wisc.edu
%
%     Date of this version:  February 2000
%
%     This MATLAB M-file implements the finite-difference time-domain
%     solution of Maxwell's curl equations over a three-dimensional
%     Cartesian space lattice comprised of uniform cubic grid cells.
%
%     To illustrate the algorithm, an air-filled rectangular cavity
%     resonator is modeled.  The length, width, and height of the
%     cavity are 10.0 cm (x-direction), 4.8 cm (y-direction), and
%     2.0 cm (z-direction), respectively.
%
%     The computational domain is truncated using PEC boundary
%     conditions:
%          ex(i,j,k)=0 on the j=1, j=jb, k=1, and k=kb planes
%          ey(i,j,k)=0 on the i=1, i=ib, k=1, and k=kb planes
%          ez(i,j,k)=0 on the i=1, i=ib, j=1, and j=jb planes
%     These PEC boundaries form the outer lossless walls of the cavity.
%
%     The cavity is excited by an additive current source oriented
%     along the z-direction.  The source waveform is a differentiated
%     Gaussian pulse given by
%          J(t)=-J0*(t-t0)*exp(-(t-t0)^2/tau^2),
%     where tau=50 ps.  The FWHM spectral bandwidth of this zero-dc-
%     content pulse is approximately 7 GHz. The grid resolution
%     (dx = 2 mm) was chosen to provide at least 10 samples per
%     wavelength up through 15 GHz.
%
%     To execute this M-file, type "fdtd3D" at the MATLAB prompt.
%     This M-file displays the FDTD-computed Ez fields at every other
%     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

%***********************************************************************
%     Grid parameters
%***********************************************************************

ie=50;       %number of grid cells in x-direction
je=24;       %number of grid cells in y-direction
ke=10;       %number of grid cells in z-direction

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

is=26;       %location of z-directed current source
js=13;       %location of z-directed current source

kobs=5;

dx=0.002;          %space increment of cubic lattice
dt=dx/(2.0*cc);    %time step

nmax=500;          %total number of time steps

%***********************************************************************
%     Differentiated Gaussian pulse excitation
%***********************************************************************

rtau=50.0e-12;
tau=rtau/dt;
ndelay=3*tau;
srcconst=-dt*3.0e+11;

%***********************************************************************
%     Material parameters
%***********************************************************************

eps=1.0;
sig=0.0;

%***********************************************************************
%     Updating coefficients
%***********************************************************************

ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps));
cb=(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps));
da=1.0;
db=dt/muz/dx;

%***********************************************************************
%     Field arrays
%***********************************************************************

ex=zeros(ie,jb,kb);
ey=zeros(ib,je,kb);
ez=zeros(ib,jb,ke);
hx=zeros(ib,je,ke);
hy=zeros(ie,jb,ke);
hz=zeros(ie,je,kb);

%***********************************************************************
%     Movie initialization
%***********************************************************************

tview(:,:)=ez(:,:,kobs);
sview(:,:)=ez(:,js,:);

subplot('position',[0.15 0.45 0.7 0.45]),pcolor(tview');
shading flat;
caxis([-1.0 1.0]);
colorbar;
axis image;
title(['Ez(i,j,k=5), time step = 0']);
xlabel('i coordinate');
ylabel('j coordinate');

subplot('position',[0.15 0.10 0.7 0.25]),pcolor(sview');
shading flat;
caxis([-1.0 1.0]);
colorbar;
axis image;
title(['Ez(i,j=13,k), time step = 0']);
xlabel('i coordinate');
ylabel('k coordinate');

rect=get(gcf,'Position');
rect(1:2)=[0 0];

M=moviein(nmax/2,gcf,rect);

%***********************************************************************
%     BEGIN TIME-STEPPING LOOP
%***********************************************************************

for n=1:nmax

%***********************************************************************
%     Update electric fields
%***********************************************************************

ex(1:ie,2:je,2:ke)=ca*ex(1:ie,2:je,2:ke)+...
                   cb*(hz(1:ie,2:je,2:ke)-hz(1:ie,1:je-1,2:ke)+...
                       hy(1:ie,2:je,1:ke-1)-hy(1:ie,2:je,2:ke));

ey(2:ie,1:je,2:ke)=ca*ey(2:ie,1:je,2:ke)+...
                   cb*(hx(2:ie,1:je,2:ke)-hx(2:ie,1:je,1:ke-1)+...
                       hz(1:ie-1,1:je,2:ke)-hz(2:ie,1:je,2:ke));

ez(2:ie,2:je,1:ke)=ca*ez(2:ie,2:je,1:ke)+...
                   cb*(hx(2:ie,1:je-1,1:ke)-hx(2:ie,2:je,1:ke)+...
                       hy(2:ie,2:je,1:ke)-hy(1:ie-1,2:je,1:ke));

ez(is,js,1:ke)=ez(is,js,1:ke)+...
               srcconst*(n-ndelay)*exp(-((n-ndelay)^2/tau^2));

%***********************************************************************
%     Update magnetic fields
%***********************************************************************

hx(2:ie,1:je,1:ke)=hx(2:ie,1:je,1:ke)+...
                   db*(ey(2:ie,1:je,2:kb)-ey(2:ie,1:je,1:ke)+...
                       ez(2:ie,1:je,1:ke)-ez(2:ie,2:jb,1:ke));

hy(1:ie,2:je,1:ke)=hy(1:ie,2:je,1:ke)+...
                   db*(ex(1:ie,2:je,1:ke)-ex(1:ie,2:je,2:kb)+...
                       ez(2:ib,2:je,1:ke)-ez(1:ie,2:je,1:ke));

hz(1:ie,1:je,2:ke)=hz(1:ie,1:je,2:ke)+...
                   db*(ex(1:ie,2:jb,2:ke)-ex(1:ie,1:je,2:ke)+...
                       ey(1:ie,1:je,2:ke)-ey(2:ib,1:je,2:ke));

%***********************************************************************
%     Visualize fields
%***********************************************************************

if mod(n,2)==0;

timestep=int2str(n);
tview(:,:)=ez(:,:,kobs);
sview(:,:)=ez(:,js,:);

subplot('position',[0.15 0.45 0.7 0.45]),pcolor(tview');
shading flat;
caxis([-1.0 1.0]);
colorbar;
axis image;
title(['Ez(i,j,k=5), time step = ',timestep]);
xlabel('i coordinate');
ylabel('j coordinate');

subplot('position',[0.15 0.10 0.7 0.25]),pcolor(sview');
shading flat;
caxis([-1.0 1.0]);
colorbar;
axis image;
title(['Ez(i,j=13,k), time step = ',timestep]);
xlabel('i coordinate');
ylabel('k coordinate');

nn=n/2;
M(:,nn)=getframe(gcf,rect);

end;

%***********************************************************************
%     END TIME-STEPPING LOOP
%***********************************************************************

end

movie(gcf,M,0,10,rect);

共1条评分

cem-uestc

rf币 +5

优秀资料+RF币

2008-10-23

逆流而上

离线cem-uestc

UID ：9061