%***********************************************************************
% 3-D FDTD code with UPML absorbing boundary conditions
%***********************************************************************
%
% Program author: Keely J. Willis, Graduate Student
% UW Computational Electromagnetics Laboratory
% Director: Susan C. Hagness
% Department of Electrical and Computer Engineering
% University of Wisconsin-Madison
% 1415 Engineering Drive
% Madison, WI 53706-1691
% kjwillis@wisc.edu
%
% Copyright 2005
%
% 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.
%
% The dimensions of the computational domain are 8.2 cm
% (x-direction), 3.4 cm (y-direction), and 3.2 cm (z-direction).
% The grid is terminated with UPML absorbing boundary conditions.
%
% An electric current source comprised of two collinear Jz components
% (realizing a Hertzian dipole) excites a radially propagating wave.
% The current source is located in the center of the grid. 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_UPML" at the MATLAB prompt.
%
% This code has been tested in the following Matlab environments:
% Matlab version 6.1.0.450 Release 12.1 (May 18, 2001)
% Matlab version 6.5.1.199709 Release 13 Service Pack 1 (August 4, 2003)
% Matlab version 7.0.0.19920 R14 (May 6, 2004)
% Matlab version 7.0.1.24704 R14 Service Pack 1 (September 13, 2004)
% Matlab version 7.0.4.365 R14 Service Pack 2 (January 29, 2005)
%
% Note: if you are using Matlab version 6.x, you may wish to make
% one or more of the following modifications to this code:
% --uncomment line numbers 485 and 486
% --comment out line numbers 552 and 561
%
%***********************************************************************
clear
%***********************************************************************
% Fundamental constants
%***********************************************************************
cc=2.99792458e8;
muz=4.0*pi*1.0e-7;
epsz=1.0/(cc*cc*muz);
etaz=sqrt(muz/epsz);
%***********************************************************************
% Material parameters
%***********************************************************************
mur=1.0;
epsr=1.0;
eta=etaz*sqrt(mur/epsr);
%***********************************************************************
% Grid parameters
%
% Each grid size variable name describes the number of sampled points
% for a particular field component in the direction of that component.
% Additionally, the variable names indicate the region of the grid
% for which the dimension is relevant. For example, ie_tot is the
% number of sample points of Ex along the x-axis in the total
% computational grid, and jh_bc is the number of sample points of Hy
% along the y-axis in the y-normal UPML regions.
%
%***********************************************************************
ie=41; % Size of main grid
je=41;
ke=16;
ih=ie+1;
jh=je+1;
kh=ke+1;
upml=10; % Thickness of PML boundaries
ih_bc=upml+1;
jh_bc=upml+1;
kh_bc=upml+1;
ie_tot=ie+2*upml; % Size of total computational domain
je_tot=je+2*upml;
ke_tot=ke+2*upml;
ih_tot=ie_tot+1;
jh_tot=je_tot+1;
kh_tot=ke_tot+1;
is=round(ih_tot/2); % Location of z-directed current source
js=round(jh_tot/2);
ks=round(ke_tot/2);
%***********************************************************************
% Fundamental grid parameters
%***********************************************************************
delta=0.002;
dt=delta*sqrt(epsr*mur)/(2.0*cc);
nmax=200;
%***********************************************************************
% Differentiated Gaussian pulse excitation
%***********************************************************************
rtau=50.0e-12;
tau=rtau/dt;
ndelay=3*tau;
J0=-1.0*epsz;
%***********************************************************************
% Initialize field arrays
%***********************************************************************
ex=zeros(ie_tot,jh_tot,kh_tot);
ey=zeros(ih_tot,je_tot,kh_tot);
ez=zeros(ih_tot,jh_tot,ke_tot);
dx=zeros(ie_tot,jh_tot,kh_tot);
dy=zeros(ih_tot,je_tot,kh_tot);
dz=zeros(ih_tot,jh_tot,ke_tot);
hx=zeros(ih_tot,je_tot,ke_tot);
hy=zeros(ie_tot,jh_tot,ke_tot);
hz=zeros(ie_tot,je_tot,kh_tot);
bx=zeros(ih_tot,je_tot,ke_tot);
by=zeros(ie_tot,jh_tot,ke_tot);
bz=zeros(ie_tot,je_tot,kh_tot);
%***********************************************************************
% Initialize update coefficient arrays
%***********************************************************************
TEST = size(ex);
C1ex=zeros(size(ex));
C2ex=zeros(size(ex));
C3ex=zeros(size(ex));
C4ex=zeros(size(ex));
C5ex=zeros(size(ex));
C6ex=zeros(size(ex));
C1ey=zeros(size(ey));
C2ey=zeros(size(ey));
C3ey=zeros(size(ey));
C4ey=zeros(size(ey));
C5ey=zeros(size(ey));
C6ey=zeros(size(ey));
C1ez=zeros(size(ez));
C2ez=zeros(size(ez));
C3ez=zeros(size(ez));
C4ez=zeros(size(ez));
C5ez=zeros(size(ez));
C6ez=zeros(size(ez));
D1hx=zeros(size(hx));
D2hx=zeros(size(hx));
D3hx=zeros(size(hx));
D4hx=zeros(size(hx));
D5hx=zeros(size(hx));
D6hx=zeros(size(hx));
D1hy=zeros(size(hy));
D2hy=zeros(size(hy));
D3hy=zeros(size(hy));
D4hy=zeros(size(hy));
D5hy=zeros(size(hy));
D6hy=zeros(size(hy));
D1hz=zeros(size(hz));
D2hz=zeros(size(hz));
D3hz=zeros(size(hz));
D4hz=zeros(size(hz));
D5hz=zeros(size(hz));
D6hz=zeros(size(hz));
%***********************************************************************
% Update coefficients, as described in Section 7.8.2.
%
% In order to simplify the update equations used in the time-stepping
% loop, we implement UPML update equations throughout the entire
% grid. In the main grid, the electric-field update coefficients of
% Equations 7.91a-f and the correponding magnetic field update
% coefficients extracted from Equations 7.89 and 7.90 are simplified
% for the main grid (free space) and calculated below.
%
%***********************************************************************
C1=1.0;
C2=dt;
C3=1.0;
C4=1.0/2.0/epsr/epsr/epsz/epsz;
C5=2.0*epsr*epsz;
C6=2.0*epsr*epsz;
D1=1.0;
D2=dt;
D3=1.0;
D4=1.0/2.0/epsr/epsz/mur/muz;
D5=2.0*epsr*epsz;
D6=2.0*epsr*epsz;
%***********************************************************************
% Initialize main grid update coefficients
%***********************************************************************
C1ex(:,jh_bc:jh_tot-upml,:)=C1;
C2ex(:,jh_bc:jh_tot-upml,:)=C2;
C3ex(:,:,kh_bc:kh_tot-upml)=C3;
C4ex(:,:,kh_bc:kh_tot-upml)=C4;
C5ex(ih_bc:ie_tot-upml,:,:)=C5;
C6ex(ih_bc:ie_tot-upml,:,:)=C6;
C1ey(:,:,kh_bc:kh
评论1
最新资源