%***********************************************************************
%?????3-D?FDTD?code?with?PML?absorbing?boundary?conditions
%***********************************************************************
%??????????????????????Lijijun?
%??????????????????????Department?of?Physics??Ynagtze?University
%??????????????????????lijijun@163.net
%???????????????????
%????Date?of?this?version:?August?2006??????????????????
%
%***********************************************************************
function?fdtd3D

clear
clc
%***********************************************************************
%?????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
etaz=sqrt（muz/epsz）;

freq=1e12;
lambda=cc/freq;
omega=2.0*pi*freq;????

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

ie=101;???????????%number?of?grid?cells?in?x-direction
if?mod（ie2）==0
???ie=ie+1;
end???????%Force?it?to?be?a?odd?integer!
je=131;????????????%number?of?grid?cells?in?y-direction
if?mod（je2）==0
???je=je+1;
end???????%Force?it?to?be?a?odd?integer!
ke=51;?????????????%number?of?grid?cells?in?z-direction
if?mod（ke2）==0
???ke=ke+1;
end???????%Force?it?to?be?a?odd?integer!

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

is=（ie+1）/2;????????????%location?of??hard?source?in?x?axis
js=（je+1）/2;??????????%location?of??hard?source?in?y?axis
ks=（ke+1）/2;??????????%location?of??hard?source?in?z?axis

ds=lambda/20;????????%space?increment?of?square?lattice
dt=ds/（2.0*cc）;???%time?step

nmax=100;%?????????%total?number?of?time?steps

iebc=8;???????????%thickness?of?left?and?right?PML?region
jebc=8;???????????%thickness?of?front?and?back?PML?region
kebc=8;???????????%thickness?of?bottom?and?top?PML?region
ibbc=iebc+1;
jbbc=jebc+1;
kbbc=kebc+1;
iefbc=ie+2*iebc;
jefbc=je+2*jebc;
kefbc=ke+2*kebc;
ibfbc=iefbc+1;
jbfbc=jefbc+1;
kbfbc=kefbc+1;

rmax=1.0e-7;
orderbc=2;

%***********************************************************************
%?????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];

ds?=?ds/sqrt（eps（1）*mur（1））;

%***********************************************************************
%?????Wave?excitation
%***********************************************************************

source=zeros（1nmax）;

for?n=1:nmax
????source（n）=sin（omega*n*dt）;
end

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

ex=zeros（iejbkb）;???????????%fields?in?main?grid?
ey=zero