FDTD(时域有限差分法)算法的Matlab源程序

1、%*%3-DFDTDcodewithPECboundaries%*%Programauthor:SusanC.Hagness%DepartmentofElectricalandComputerEngineering%UniversityofWisconsin-Madison%1415EngineeringDrive%Madison,WI53706-1691%608-265-5739*.edu*.edu%Dateofthisversion:February2000%ThisMATLABM-fileimplementsthefinite-differencetime-domain%solution

2、ofMaxwellscurlequationsoverathree-dimensional%Cartesianspacelatticecomprisedofuniformcubicgridcells.%Toillustratethealgorithm,anair-filledrectangularcavity%resonatorismodeled.Thelength,width,andheightofthe%cavityare10.0cm(x-direction),4.8cm(y-direction),and%2.0cm(z-direction),respectively.%Thecomput

3、ationaldomainistruncatedusingPECboundary%conditions:%ex(i,j,k)=0onthej=1,j=jb,k=1,andk=kbplanes%ey(i,j,k)=0onthei=1,i=ib,k=1,andk=kbplanes%ez(i,j,k)=0onthei=1,i=ib,j=1,andj=jbplanes%ThesePECboundariesformtheouterlosslesswallsofthecavity.%Thecavityisexcitedbyanadditivecurrentsourceoriented%alongthez-

4、direction.Thesourcewaveformisadifferentiated%Gaussianpulsegivenby%J(t)=-J0*(t-t0)*exp(-(t-t0)A2/tauA2),%wheretau=50ps.TheFWHMspectralbandwidthofthiszero-dc-%contentpulseisapproximately7GHz.Thegridresolution%(dx=2mm)waschosentoprovideatleast10samplesper%wavelengthupthrough15GHz.%ToexecutethisM-file,t

5、ypefdtd3DattheMATLABprompt.%ThisM-filedisplaystheFDTD-computedEzfieldsateveryother%timestep,andrecordsthoseframesinamoviematrix,M,which%isplayedattheendofthesimulationusingthemoviecommand.%*clear%*%Fundamentalconstantscc=2.99792458e8;muz=4.0*pi*1.0e-7;cc=2.99792458e8;muz=4.0*pi*1.0e-7;epsz=1.0/(cc*c

6、c*muz);%speedoflightinfreespace%permeabilityoffreespace%permittivityoffreespace%*%Gridparameters%*ie=50;%numberofgridcellsinx-directionje=24;%numberofgridcellsiny-directionke=10;%numberofgridcellsinz-directionib=ie+1;jb=je+1;kb=ke+1;is=26;%locationofz-directedcurrentsourcejs=13;%locationofz-directed

7、currentsourcekobs=5;dx=0.002;%spaceincrementofcubiclatticedt=dx/(2.0*cc);%timestepnmax=500;%totalnumberoftimesteps%*%DifferentiatedGaussianpulseexcitation%*rtau=50.0e-12;tau=rtau/dt;ndelay=3*tau;srcconst=-dt*3.0e+11;%*%Materialparameters%*eps=1.0;sig=0.0;%*%Updatingcoefficients%*ca=(1.0-(dt*sig)/(2.

8、0*epsz*eps)/(1.0+(dt*sig)/(2.0*epsz*eps);cb=(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps);da=1.0;db=dt/muz/dx;%*%Fieldarrays%*ex=zeros(ie,jb,kb);ey=zeros(ib,je,kb);ez=zeros(ib,jb,ke);hx=zeros(ib,je,ke);hy=zeros(ie,jb,ke);hz=zeros(ie,je,kb);%*%Movieinitialization%*tview(:,:)=ez(:,:,kobs);sview(:,:)=e

9、z(:,js,:);subplot(position,0.150.450.70.45),pcolor(tview);shadingflat;caxis(-1.01.0);colorbar;axisimage;title(Ez(i,j,k=5),timestep=0);xlabel(icoordinate);ylabel(jcoordinate);subplot(position,0.25),pcolor(sview);shadingflat;caxis(-1.01.0);colorbar;axisimage;title(Ez(i,j=13,k),timestep=0);x

10、label(icoordinate);ylabel(kcoordinate);rect=get(gcf,Position);rect(1:2)=00;M=moviein(nmax/2,gcf,rect);%*%BEGINTIME-STEPPINGLOOP%*forn=1:nmax%*%Updateelectricfields%*ex(1:ie,2:je,2:ke)=ca*ex(1:ie,2:je,2:ke)+.cb*(hz(1:ie,2:je,2:ke)-hz(1:ie,1:je-1,2:ke)+.hy(1:ie,2:je,1:ke-1)-hy(1:ie,2:je,2:ke);ey(2:ie,

11、1:je,2:ke)=ca*ey(2:ie,1:je,2:ke)+.cb*(hx(2:ie,1:je,2:ke)-hx(2:ie,1:je,1:ke-1)+.hz(1:ie-1,1:je,2:ke)-hz(2:ie,1:je,2:ke);ez(2:ie,2:je,1:ke)=ca*ez(2:ie,2:je,1:ke)+.cb*(hx(2:ie,1:je-1,1:ke)-hx(2:ie,2:je,1:ke)+.hy(2:ie,2:je,1:ke)-hy(1:ie-1,2:je,1:ke);ez(is,js,1:ke)=ez(is,js,1:ke)+.srcconst*(n-ndelay)*exp

12、(-(n-ndelayF2/tauA2);%*%Updatemagneticfields%*hx(2:ie,1:je,1:ke)=hx(2:ie,1:je,1:ke)+.db*(ey(2:ie,1:je,2:kb)-ey(2:ie,1:je,1:ke)+.ez(2:ie,1:je,1:ke)-ez(2:ie,2:jb,1:ke);hy(1:ie,2:je,1:ke)=hy(1:ie,2:je,1:ke)+.db*(ex(1:ie,2:je,1:ke)-ex(1:ie,2:je,2:kb)+.ez(2:ib,2:je,1:ke)-ez(1:ie,2:je,1:ke);hz(1:ie,1:je,2

13、:ke)=hz(1:ie,1:je,2:ke)+.db*(ex(1:ie,2:jb,2:ke)-ex(1:ie,1:je,2:ke)+.ey(1:ie,1:je,2:ke)-ey(2:ib,1:je,2:ke);%*%Visualizefields%*ifmod(n,2)=0;timestep=int2str(n);tview(:,:)=ez(:,:,kobs);sview(:,:)=ez(:,js,:);subplot(position,0.150.450.70.45),pcolor(tview);shadingflat;caxis(-1.01.0);colorbar;axisimage;title(Ez(i,j,k=5),timestep=,timestep);xlabel(icoordinate);ylabel(jcoordinate);subplot(position,0.25),pcolor(sview);sha