《三维电磁场FDTD程序PEC边界.docx》由会员分享，可在线阅读，更多相关《三维电磁场FDTD程序PEC边界.docx（15页珍藏版）》请在三一办公上搜索。
1、三维电磁场FDTD程序PEC边界三维电磁场FDTD程序PEC边界 % FDTD Main Function Jobs to Workers % %* % 3-D FDTD code with PEC boundaries %* % % This MATLAB M-file implements the finite-difference time-domain % solution of Maxwells curl equations over a three-dimensional % Cartesian space lattice comprised of uniform cubic gr
2、id cells. % % To illustrate the algorithm, an air-filled rectangular cavity % resonator is modeled. The length, width, and height of the % cavity are X cm (x-direction), Y cm (y-direction), and % Z cm (z-direction), respectively. % % The computational domain is truncated using PEC boundary % conditi
3、ons: % ex(i,j,k)=0 on the j=1, j=jb, k=1, and k=kb planes % ey(i,j,k)=0 on the i=1, i=ib, k=1, and k=kb planes % ez(i,j,k)=0 on the i=1, i=ib, j=1, and j=jb planes % These PEC boundaries form the outer lossless walls of the cavity. % % The cavity is excited by an additive current source oriented % a
4、long the z-direction. The source waveform is a differentiated % Gaussian pulse given by % J(t)=-J0*(t-t0)*exp(-(t-t0)2/tau2), % 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 sa
5、mples per % wavelength up through 15 GHz. % % To execute this M-file, type fdtd3D at the MATLAB prompt. % This M-file displays the FDTD-computed Ez fields at every other % time step, and records those frames in a movie matrix, M, which % is played at the end of the simulation using the movie command
6、. % %* function Ex,Ey,Ez=FDTD3D_Main(handles) global SimRunStop % if isdir(C:MATLAB7workcavityfigures) % mkdir C:MATLAB7workcavityfigures % end %* % Grid Partition %* p.ip = get(handles.XdirPar,Value); p.jp = get(handles.YdirPar,Value); p.PN = get(handles.PartNo,Value); %* % Grid Dimensons %* ie = g
7、et(handles.xslider,Value); %number of grid cells in x-direction je = get(handles.yslider,Value); %number of grid cells in y-direction ke = get(handles.zslider,Value); %number of grid cells in z-direction ib=ie+1; jb=je+1; kb=ke+1; %* % All Domains Fields Ini. %* Ex=zeros(ie,jb,kb); Ey=zeros(ib,je,kb
8、); Ez=zeros(ib,jb,ke); Hx=zeros(ib,je,ke); Hy=zeros(ie,jb,ke); Hz=zeros(ie,je,kb); %* % Fundamental constants %* param.cc=2.99792458e8; %speed of light in free space param.muz=4.0*pi*1.0e-7; %permeability of free space param.epsz = 1.0/(param.cc*param.cc*param.muz); %permittivity of free space %* %
9、Grid parameters %* param.is = get(handles.xsource,Value); %location of z-directed current source param.js = get(handles.ysource,Value); %location of z-directed current source param.kobs = floor(ke/2); %Surface of observation param.dx = get(handles.CellSize,Value); %space increment of cubic lattice p
10、aram.dt = param.dx/(2.0*param.cc); %time step param.nmax = get(handles.TimeStep,Value); %total number of time steps %* % Differentiated Gaussian pulse excitation %* param.rtau=get(handles.tau,Value)*100e-12; param.tau=param.rtau/param.dt; param.ndelay=3*param.tau; param.Amp = get(handles.sourceamp,V
11、alue)*10e11; param.srcconst=-param.dt*param.Amp; %* % Material parameters %* param.eps= get(handles.epsilon,Value); param.sig= get(handles.sigma,Value); %* % Updating coefficients %* 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); p
12、aram.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; %* % Calling FDTD Algorithm %* ex=zeros(ib,jb,kb); ey=zeros(ib,jb,kb); ez=zeros(ib,jb,kb); hx=zeros(ib,jb,kb); hy=zeros(ib,jb,kb); hz=zeros(ib,jb
13、,kb); X,Y,Z = meshgrid(1:ib,1:jb,1:kb); % Grid coordinates Psim = 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 n=1:1:param.nmax for m1=1:1:p.PN ex,ey,ez=Efields(para
14、m,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); end Psim(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_FFT pause(0.01); end elseif (p.ip = 0)&(p.jp = 1
15、) 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; for n=1:1:param.nmax for m2=1:1:p.PN ex,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); end Ps
16、im(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); end elseif (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; for n=1:1:param.n
17、max for m2=1:1:p.PN for m1=1:1:p.PN ex,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); end end Psim(n),Panl(n) = Cavity_Power(param,handles,ex,ey,ez,n); field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,
18、Panl,p); pause(0.01); end else param.a = 1; param.b = ie; param.c = 1; param.d = je; m1 = 1;m2=1; for n=1:1:param.nmax ex,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
19、_sim,value); if SimRunStop = 1 h = warndlg(Simulation Run is Stopped by User !,! Warning !); waitfor(h); break; end Psim(n),Panl(n) = Cavity_Power(param,handles,ex,ey,ez,n); if n=2 Panl(n)=Panl(n) + Panl(n-1); end field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p); pause(0.01); end end %文件2 % Cav
20、ity Field Viz. % function = field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p) %* % Cross-Section initialization %* tview = squeeze(ez(:,:,param.kobs); sview = squeeze(ez(:,param.js,:); ax1 = handles.axes1; ax2 = handles.axes2; ax3 = handles.axes3; %* % Visualize fields %* timestep=int2str(n); ez
21、view=squeeze(ez(:,:,param.kobs); exview=squeeze(ex(:,:,param.kobs); eyview=squeeze(ey(:,:,param.kobs); xmin = min(X(:); xmax = max(X(:); ymin = 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; cx cy cz = meshg
22、rid(xrange,yrange,zrange); % sview=squeeze(ez(:,param.js,:); rect = -50 -35 360 210; rectp = -50 -40 350 260; axes(ax1); axis tight set(gca,nextplot,replacechildren); E_total = sqrt(ex.2 + ey.2 + ez.2); etview=squeeze(E_total(:,:,param.kobs); sview = squeeze(E_total(:,param.js,:); surf(etview); shad
23、ing interp; caxis(-1.0 1.0); colorbar; axis image; title(Total E-Field, time step = ,timestep,fontname,Times New Roman,fontsize,12,FontWeight,BOLD); xlabel(i coordinate); ylabel(j coordinate); set(gca,fontname,Times New Roman,fontsize,10); % F1 = getframe(gca,rect); % M1 = frame2im(F1); % filename =
24、 fullfile(C:MATLAB7workcavityfigures,strcat(XY_ETotal,num2str(n),.jpeg); % imwrite(M1,filename,jpeg) axes(ax2); axis tight set(gca,nextplot,replacechildren); surf(sview); shading interp; caxis(-1.0 1.0); colorbar; axis image; title(Ez(i,j=13,k), time step = ,timestep,fontname,Times New Roman,fontsiz
25、e,12,FontWeight,BOLD); xlabel(i coordinate); ylabel(k coordinate); set(gca,fontname,Times New Roman,fontsize,10); % F2 = getframe(gca,rect); % M2 = frame2im(F2); % filename = fullfile(C:MATLAB7workcavityfigures,strcat(XZ_ETotal,num2str(n),.jpeg); % imwrite(M2,filename,jpeg) %* % Cavity Power - Anali
26、tic expression %* axes(ax3); % axis tight % 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) = Normalized Analytic Vs. Simulated Power; str(2) = as function of time-steps; title(str,fontname,Times N
27、ew Roman,fontsize,12,FontWeight,BOLD ); xlabel(Time Step); ylabel(Log(Power); legend(Simulation,Analytic,location,SouthEast); set(gca,fontname,Times New Roman,fontsize,10); % F3 = getframe(gca,rectp); % M3 = frame2im(F3); % filename = fullfile(C:MATLAB7workcavityfigures,strcat(Power,num2str(n),.jpeg
28、); % imwrite(M3,filename,jpeg) %文件3 % Source waveform function function source=waveform(param,handles,n) ax1 = handles.axes1; type = get(handles.sourcetype,value); amp = get(handles.sourceamp,value)*1e4; f = get(handles.Frequency,value)*1e9; switch type case 2 source = param.srcconst*(n-param.ndelay
29、)*exp(-(n-param.ndelay)2/param.tau2); case 1 source = param.srcconst*sin(param.dt*n*2*pi*f); case 3 source = param.srcconst*exp(-(n-param.ndelay)2/param.tau2)*sin(2*pi*f*(n-param.ndelay)*param.dt); otherwise source = param.srcconst*(n-param.ndelay)*exp(-(n-param.ndelay)2/param.tau2); end %文件4 functi
30、on hx,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(
31、a:b,c+1:d,1:ke)+. param.db*(ex(a:b,c+1:d,1:ke)-. ex(a:b,c+1:d,2:kb)+. ez(a+1:b+1,c+1:d,1:ke)-. ez(a:b,c+1:d,1:ke); hz(a:b,c:d,2:ke)=. hz(a:b,c:d,2:ke)+. param.db*(ex(a:b,c+1:d+1,2:ke)-. ex(a:b,c:d,2:ke)+. ey(a:b,c:d,2:ke)-. ey(a+1:b+1,c:d,2:ke); %文件5 function ex,ey,ez=Efields(param,handles,ex,ey,ez,
32、hx,hy,hz,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); if (param.is=a)&(param.is=c)&(param.js=d) source = waveform(param,handles,n); else source = 0; end ex(a:b,c+1:d,2:ke)=. param.ca*ex(a:b,c+1:d,2:ke)+. param.cb*(hz(a:b,c+1:d,2:ke)-. hz(a:b,c:d-1,2:ke)+
33、. hy(a:b,c+1:d,1:ke-1)-. hy(a:b,c+1:d,2:ke); ey(a+1:b,c:d,2:ke)=. param.ca*ey(a+1:b,c:d,2:ke)+. param.cb*(hx(a+1:b,c:d,2:ke)-. hx(a+1:b,c:d,1:ke-1)+. hz(a:b-1,c:d,2:ke)-. hz(a+1:b,c:d,2:ke); ez(a+1:b,c+1:d,1:ke)=. param.ca*ez(a+1:b,c+1:d,1:ke)+. param.cb*(hx(a+1:b,c:d-1,1:ke)-. hx(a+1:b,c+1:d,1:ke)+. hy(a+1:b,c+1:d,1:ke)-. hy(a:b-1,c+1:d,1:ke); ez(param.is,param.js,1:ke)=ez(param.is,param.js,1:ke)+source;
