FDTD算法源代码

%***********************************************************************
%     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（ijk）=0 on the j=1 j=jb k=1 and k=kb planes
%          ey（ijk）=0 on the i=1 i=ib k=1 and k=kb planes
%          ez（ijk）=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