fdtd_matlab 源码

%*********************************************************************** % 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 % [email protected] % % 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=17; 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=100; %*********************************************************************** % 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 %*********************************************************************** 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_tot-upml)=C1; C2ey(:,:,kh_bc:kh_tot-upml)=C2; C3ey(ih_bc:ih_tot-upml,:,:)=C3; C4ey(ih_bc:ih_tot-upml,:,:)=C4; C5ey(:,jh_bc:je_tot-upml,:)=C5; C6ey(:,jh_bc:je_tot-upml,:)=C6; C1ez(ih_bc:ih_tot-upml,:,:)=C1; C2ez(ih_bc:ih_tot-upml,:,:)=C2; C3ez(:,jh_bc:jh_tot-upml,:)=C3; C4ez(:,jh_bc:jh_tot-upml,:)=C4; C5ez(:,:,kh_bc:ke_tot-upml)=C5; C6ez(:,:,kh_bc:ke_tot-upml)=C6; D1hx(:,jh_bc:je_tot-upml,:)=D1; D2hx(:,jh_bc:je_tot-upml,:)=D2; D3hx(:,:,kh_bc:ke_tot-upml)=