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

%***********************************************************************%3-DFDTDcodewithPECboundaries%***********************************************************************%%Programauthor:SusanC.Hagness%DepartmentofElectricalandComputerEngineering%UniversityofWisconsin-Madison%1415EngineeringDrive%Madison,[e-mail]%%Dateofthisversion:February2000%%ThisMATLABM-fileimplementsthefinite-differencetime-domain%solutionofMaxwell'scurlequationsoverathree-dimensional%Cartesianspacelatticecomprisedofuniformcubicgridcells.%%Toillustratethealgorithm,anair-filledrectangularcavity%resonatorismodeled.Thelength,width,andheightofthe%cavityare10.0cm(x-direction),4.8cm(y-direction),and%2.0cm(z-direction),respectively.%%ThecomputationaldomainistruncatedusingPECboundary%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-direction.Thesourcewaveformisadifferentiated%Gaussianpulsegivenby%J(t)=-J0*(t-t0)*exp(-(t-t0)^2/tau^2),%wheretau=50ps.TheFWHMspectralbandwidthofthiszero-dc-%contentpulseisapproximately7GHz.Thegridresolution%(dx=2mm)waschosentoprovideatleast10samplesper%wavelengthupthrough15GHz.%%ToexecutethisM-file,typefdtd3DattheMATLABprompt.%ThisM-filedisplaystheFDTD-computedEzfieldsateveryother%timestep,andrecordsthoseframesinamoviematrix,M,which%isplayedattheendofthesimulationusingthemoviecommand.%%***********************************************************************clear%***********************************************************************%Fundamentalconstants%***********************************************************************cc=2.99792458e8;%speedoflightinfreespacemuz=4.0*pi*1.0e-7;%permeabilityoffreespaceepsz=1.0/(cc*cc*muz);%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-directedcurrentsourcekobs=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.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(:,:)=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=0']);xlabel('icoordinate');ylabel('jcoordinate');subplot('position',[0.150.100.70.25]),pcolor(sview');shadingflat;caxis([-1.01.0]);colorbar;axisimage;title(['Ez(i,j=13,k),timestep=0']);xlabel('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,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(-((n-ndelay)^2/tau^2));%***********************************************************************%Updatemagneticfields%*********************************************************************