三维电磁场FDTD程序PEC边界.docx

三维电磁场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 Maxwell's curl equations over a three-dimensional % Cartesian space lattice comprised of uniform cubic grid 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 % conditions: % 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 % along 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 samples 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. % %* 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 = get(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); 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 %* % 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 param.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,'Value')*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); 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; %* % 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,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(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 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) 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 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); 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.nmax 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,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_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 % Cavity 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); ezview=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 = meshgrid(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'); shading 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 = 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','fontsize',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 - Analitic 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 New 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'); % 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)*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 function 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(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,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<=b)|(param.js>=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)+. 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;