三維電磁場FDTD程序(3DFDTDmatlab)PEC邊界

1、%三維電磁場FDTD程序(3D-FDTDmatlab)PEC邊界%FDTDMainFunctionJobstoWorkers%*%3-DFDTDcodewithPECboundaries%*%ThisMATLABM-fileimplementsthefinite-differencetime-domain%solutionofMaxwellscurlequationsoverathree-dimensional%Cartesianspacelatticecomprisedofuniformcubicgridcells.%Toillustratethealgorithm,anair-filled

2、rectangularcavity%resonatorismodeled.Thelength,width,andheightofthe%cavityareXcm(x-direction),Ycm(y-direction),and%Zcm(z-direction),respectively.%ThecomputationaldomainistruncatedusingPECboundaryconditions:ex(i,j,k)=0onthej=1,j=jb,k=1,andk=kbplanesey(i,j,k)=0onthei=1,i=ib,k=1,andk=kbplanesez(i,j,k)=

3、0onthei=1,i=ib,j=1,andj=jbplanesThesePECboundariesformtheouterlosslesswallsofthecavity.Thecavityisexcitedbyanadditivecurrentsourceorientedalongthez-direction.ThesourcewaveformisadifferentiatedGaussianpulsegivenbyJ(t)=-J0*(t-tO)*exp(-(t-tO廠2/ta2),wheretau=50ps.TheFWHMspectralbandwidthofthiszero-dc-co

4、ntentpulseisapproximately7GHz.Thegridresolution(dx=2mm)waschosentoprovideatleast1Osamplesperwavelengthupthrough15GHz.ToexecutethisM-file,typefdtd3DattheMATLABprompt.ThisM-filedisplaystheFDTD-computedEzfieldsateveryothertimestep,andrecordsthoseframesinamoviematrix,M,whichisplayedattheendofthesimulati

5、onusingthemoviecommand.%*functionEx,Ey,Ez=FDTD3D_Main(handles)globalSimRunStop%ifisdir(C:MATLAB7workcavityfigures)%mkdirC:MATLAB7workcavityfigures%end%*%GridPartition%*p.ip=get(handles.XdirPar,Value);p.jp=get(handles.YdirPar,Value);p.PN=get(handles.PartNo,Value);%*%GridDimensons%*%numberofgridcellsi

6、n%numberofgridcellsin%numberofgridcellsinie=get(handles.xslider,Value);x-directionje=get(handles.yslider,Value);y-directionke=get(handles.zslider,Value);z-directionib=ie+1;jb=je+1;kb=ke+1;%*%AllDomainsFieldsIni.%*Ex=zeros(ie,jb,kb);Ey=zeros(ib,je,kb);Ez=zeros(ib,jb,ke);Hx=zeros(ib,je,ke);Hy=zeros(ie

7、,jb,ke);Hz=zeros(ie,je,kb);%*%Fundamentalconstants%*param.cc=2.99792458e8;%speedoflightinfreespaceparam.muz=4.0*pi*1.0e-7;%permeabilityoffreespaceparam.epsz=1.0/(param.cc*param.cc*param.muz);%permittivityoffreespace%*%Gridparameters%*param.is=get(handles.xsource,Value);%locationofz-directedcurrentso

8、urceparam.js=get(handles.ysource,Value);%locationofz-directedcurrentsourceparam.kobs=floor(ke/2);%Surfaceofobservationparam.dx=get(handles.CellSize,Value);%spaceincrementofcubiclatticeparam.dt=param.dx/(2.0*param.cc);%timestepparam.nmax=get(handles.TimeStep,Value);%totalnumberoftimesteps%*%Different

9、iatedGaussianpulseexcitation%*param.rtau=get(handles.tau,Value)*100e-12;param.tau=param.rtau/param.dt;param.ndelay=3*param.tau;param.Amp=get(handles.sourceamp,Value)*10e11;param.srcconst=-param.dt*param.Amp;%*%Materialparameters%*param.eps=get(handles.epsilon,Value);param.sig=get(handles.sigma,Value

10、);%*%Updatingcoefficients%*param.ca=(1.0-(param.dt*param.sig)/(2.0*param.epsz*param.eps)/(1.0+(param.dt*param.sig)/(2.0*param.epsz*param.eps);param.cb=(param.dt/param.epsz/param.eps/param.dx)/(1.0+(param.dt*param.sig)/(2.0*param.epsz*param.eps);param.da=1.0;param.db=param.dt/param.muz/param.dx;%*%Ca

11、llingFDTDAlgorithm%*bbbjjj)bbbkex=zeros(ibey=zeros(ibez=zeros(ibhx=zeros(ibhy=zeros(ibhz=zeros(ibX,Y,Z=meshgrid(1:ib,1:jb,1:kb);%GridcoordinatesPsim=zeros(param.nmax,1);Panl=zeros(param.nmax,1);if(p.ip=1)&(p.jp=0)x=ceil(ie/p.PN)param.a=1:x-1:ie-x;param.b=x+1:x-1:ie;param.c=1,1;param.d=je,je;m2=1;for

12、n=1:1:param.nmaxform1=1:1:p.PNex,ey,ez=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p);hx,hy,hz=Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p);endPsim(n),Panl(n)=Cavity_Power(param,handles,ex,ey,ez,n);field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p);Dyn_FFTpau

13、se(0.01);endelseif(p.ip=0)&(p.jp=1)y=ceil(je/p.PN);param.c=1:y-1:je-y;param.d=y+1:y-1:je;param.a=1,1;param.b=ie,ie;m1=1;forn=1:1:param.nmaxform2=1:1:p.PNex,ey,ez=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p);hx,hy,hz=Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p)

14、;endPsim(n),Panl(n)=Cavity_Power(param,handles,ex,ey,ez,n);field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p);pause(0.01);endelseif(p.ip=1)&(p.jp=1)x=ceil(ie/p.PN);param.a=1:x-1:ie-x;param.b=x+1:x-1:ie;y=ceil(je/p.PN);param.c=1:y-1:je-y;param.d=y+1:y-1:je;forn=1:1:param.nmaxform2=1:1:p.PNform1=1:

15、1:p.PNex,ey,ez=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p);hx,hy,hz=Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p);endendPsim(n),Panl(n)=Cavity_Power(param,handles,ex,ey,ez,n);field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p);pause(0.01);endelseparam.a=1;pa

16、ram.b=ie;param.c=1;param.d=je;m1=1;m2=1;forn=1:1:param.nmaxex,ey,ez=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p);hx,hy,hz=Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p);SimRunStop=get(handles.Stop_sim,value);ifSimRunStop=1h=warndlg(SimulationRunisStoppedbyUser!,

17、!Warning!);waitfor(h);break;endPsim(n),Panl(n)=Cavity_Power(param,handles,ex,ey,ez,n);ifn=2Panl(n)=Panl(n)+Panl(n-1);endfield_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p);pause(0.01);endend%文件2%CavityFieldViz.%function=field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p)%*%Cross-Sectioninitializ

18、ation%*tview=squeeze(ez(:,:,param.kobs);sview=squeeze(ez(:,param.js,:);ax1=handles.axes1;ax2=handles.axes2;ax3=handles.axes3;%*%Visualizefields%*timestep=int2str(n);ezview=squeeze(ez(:,:,param.kobs);exview=squeeze(ex(:,:,param.kobs);eyview=squeeze(ey(:,:,param.kobs);xmin=min(X(:);xmax=max(X(:);ymin=

19、min(Y(:);ymax=max(Y(:);zmin=min(Z(:);daspect(2,2,1)xrange=linspace(xmin,xmax,8);yrange=linspace(ymin,ymax,8);zrange=3:4:15;cxcycz=meshgrid(xrange,yrange,zrange);%sview=squeeze(ez(:,param.js,:);rect=-50-35360210;rectp=-50-40350260;axes(ax1);axistightset(gca,nextplot,replacechildren);E_total=sqrt(ex.2

20、+ey.2+ez.2);etview=squeeze(E_total(:,:,param.kobs);sview=squeeze(E_total(:,param.js,:);surf(etview);shadinginterp;caxis(-1.01.0);colorbar;axisimage;title(TotalE-Field,timestep=,timestep,fontname,TimesNewRoman,fontsize,12,FontWeight,BOLD);xlabel(icoordinate);ylabel(jcoordinate);set(gca,fontname,Times

21、NewRoman,fontsize,10);%F1=getframe(gca,rect);%M1=frame2im(F1);%filename=fullfile(C:MATLAB7workcavityfigures,strcat(XY_ETotal,num2str(n),.jpeg);%imwrite(M1,filename,jpeg)axes(ax2);axistightset(gca,nextplot,replacechildren);surf(sview);shadinginterp;caxis(-1.01.0);colorbar;axisimage;title(Ez(i,j=13,k)

22、,timestep=,timestep,fontname,TimesNewRoman,fontsize,12,FontWeight,BOLD);xlabel(icoordinate);ylabel(kcoordinate);set(gca,fontname,TimesNewRoman,fontsize,10);%F2=getframe(gca,rect);%M2=frame2im(F2);%filename=fullfile(C:MATLAB7workcavityfigures,strcat(XZ_ETotal,num2str(n),.jpeg);%imwrite(M2,filename,jp

23、eg)%*%CavityPower-Analiticexpression%*axes(ax3);%axistight%set(gca,nextplot,replacechildren);t=1:1:param.nmax;Psim=1e3*Psim./max(Psim);Panl=1e3*Panl./max(Panl);semilogy(t,Psim,b.-,t,Panl,rx-);str(1)=NormalizedAnalyticVs.SimulatedPower;str(2)=asfunctionoftime-steps;title(str,fontname,TimesNewRoman,fo

24、ntsize,12,FontWeight,BOLD);xlabel(TimeStep);ylabel(Log(Power);legend(Simulation,Analytic,location,SouthEast);set(gca,fontname,TimesNewRoman,fontsize,10);%F3=getframe(gca,rectp);%M3=frame2im(F3);%filename=fullfile(C:MATLAB7workcavityfigures,strcat(Power,num2str(n),.jpeg);%imwrite(M3,filename,jpeg)%文件

25、3%Sourcewaveformfunctionfunctionsource=waveform(param,handles,n)ax1=handles.axes1;type=get(handles.sourcetype,value);amp=get(handles.sourceamp,value)*1e4;f=get(handles.Frequency,value)*1e9;switchtypecase2source=param.srcconst*(n-param.ndelay)*exp(-(n-param.ndelay廠2/param.tau2);case1source=param.srcc

26、onst*sin(param.dt*n*2*pi*f);case3source=param.srcconst*exp(-(n-param.ndelay廠2/param.ta2)*sin(2*pi*f*(n-param.ndelay)*param.dt);otherwisesource=param.srcconst*(n-param.ndelay)*exp(-(n-param.ndelay廠2/param.tau2);end%文件4functionhx,hy,hz=Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,x,y,p)a=param.a(x);b=param.b(x);c=param.c(y);d=param.d(y);hx(a+1:b,c:d,1:ke)=.hx(a+1:b,c:d,1:ke)+.param.db*(ey(a+1:b,c:d,2:kb)-.ey(a+1:b,c:d,1:ke)+.ez(a+1:b,c:d,1:ke)-.ez(a+1:b,c+1:d+1,1:ke);hy(a:b,c+1:d,1:ke)=.hy(a:b,c+1:d,1:ke)+.param.db*(ex(a:b,c+1:d,1:ke)-.ex(a:b,c+1:d