2D FDTD TE mode with a plane wave source and a PML ABC

2D FDTD TE mode with a plane wave source and a PML abc

%2D FDTD TE mode with a plane wave source and a PML abc
%program edited by: K.B.A., NIP,U.P.,Diliman
%Date: January 2008
%this program borrowed the (2D FDTD)code by susan hagness
%***********************************************************************
%CHANGES: The metallic cylindrical scatterer is removed.
%          The whole computational region becomes free space only.
%          Program implements the TF/SF source formulation
%***********************************************************************
%Note:if you want to include the cylindrical metal, replace the sig
% expression with the commented sig in line 124
%
%***********************************************************************
%     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;

is1=7;              %connecting boundaries
is2=ie-is1-1;          
js1=7;
js2=je-js1-1;

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

media=2;

eps=[1.0 1.0];
sig=[0.0 0.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
ey_inc=zeros(ib,1);
hz_inc=zeros(ie,1);
ey_low_m1=0;
ey_low_m2=0;
ey_high_m1=0;
ey_high_m2=0;

%***********************************************************************
%     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));

%**************************************************************************       
%ex field correction(all from is1 to is2-1)
%ex(js1)_total=ex(js1)_tot+(dt/(epsz*dx))*(hz(js1)_total-hz(js1-1)_scattered)
%hz(js1)_total-hz(js1-1)_scattered is not consistent because
%hz(js1-1)_scattered is scattered not total. so we add hz_inc to hz(js1-1)
ex(is1:is2-1,js1)=ex(is1:is2-1,js1)-((dt/(epsz*dx)).*hz_inc(is1:is2-1));
ex(is1:is2-1,js2)=ex(is1:is2-1,js2)+((dt/(epsz*dx)).*hz_inc(is1:is2-1));
%*************************************************************************

%1D buffer
ey_inc(2:ib-1)=ey_inc(2:ib-1)-(dt/(epsz*dx)).*(hz_inc(2:ie)-hz_inc(1:ie-1));
ey_inc(1)=ey_low_m2;
ey_low_m2=ey_low_m1;
ey_low_m1=ey_inc(2);

ey_inc(ib)=ey_high_m2;
ey_high_m2=ey_high_m1;
ey_high_m1=ey_inc(ib-1);

ey_inc(3)=ey_inc(3)+sin(omega*dt*n);

ey(2:ie,:)=caey(2:ie,:).*ey(2:ie,:)+...
           cbey(2:ie,:).*(hz(1:ie-1,:)-hz(2:ie,:));
       
%***********************************************************************
%ey field correction
ey(is1,js1:js2-1)=ey(is1,js1:js2-1)+(dt/(epsz*dx)).*hz_inc(is1-1);
ey(is2,js1:js2-1)=ey(is2,js1:js2-1)-(dt/(epsz*dx)).*hz_inc(is2);
%***********************************************************************

%***********************************************************************
%     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_inc(1:ie)=hz_inc(1:ie)-(dt/(muz*dx)).*(ey_inc(2:ib)-ey_inc(1:ib-1));

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);
%************************************************************************
%hz field correction
hz(is1-1,js1:js2-1)=hz(is1-1,js1:js2-1)+(dt/(muz*dx)).*ey_inc(is1);
hz(is2,js1:js2-1)=hz(is2,js1:js2-1)-(dt/(muz*dx)).*ey_inc(is2);
%************************************************************************


%***********************************************************************
%     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([-1.3 1.3]);
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([-1.3 1.3]);
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.003 0.003]);
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);


－－－－－－－－－－－－－
附一个国外doctor的观点
I "think" the methods you are referring to are:

Finite-Difference Frequency-Domain ( FDFD)
Finite-Difference Time-Domain (FDTD)
Transmission Line Method (TLM)
Time-Domain Transmission Line Method (TLM)

All four of these methods typically work by fitting the fields and materials to discrete points on a Cartesion grid. The formulation of equations relating the fields through Maxwell's equations is slightly different in each method, but they are all essentially just enforcing Maxwell's equations across the grid.

Time-domain methods have the ability to simulate a device over a very broad frequency range in just a single simulation. A pulse source is typically used so at all points of interest the impulse response is recorded. Fourier transforming the impulse response gives you the spectral response. FDTD has another big advantage in that it does use linear algebra. This mean that if you double the size of your problem, memory and run time only doubles. Other methods, doubling the size of the problem could increase your memory and run time by 4 to 16 times or more! Time-domain methods, however, are poor for modeling highly resonant devices, or devices where the spectral response changes abruptly. In these cases, the frequency-domain methods tend to be more accurate.

Frequency-domain methods tend to be more efficient and accurate for highly resonant devices or when you are only interested in one or two frequency points. They are based on linear algebra so memory and run time increases exponentially with problem size. Typically you convert Maxwell's equations to either A*x=b for scattering problems, or A*x=a*x for eigen-value problems.

I will also say that FDFD is the simplest numerical method I know to implement. It is not as efficient as the finite element method (FEM), but FEM is more tedious to formulate. I will also say that I do not see TLM or TDTLM used nearly as much as FDTD so there is much more literate available for FDTD. Most of this is also applicable to FDFD, but otherwise FDFD can be difficult to find in the literature.

In my Ph.D. dissertation, I describe in good detail a number of numerical methods including:

Finite-Difference Frequency-Domain ( FDFD)
Finite-Difference Time-Domain (FDTD)
Plane Wave Expansion Method (PWEM)
Rigorous-Coupled Wave Analysis (RCWA)
String Method
Level Set Method (LSM)
Fast Marching Method (FMM)

You can download it here:

http://purl.fcla.edu/fcla/etd/CFE0001159