fdtd.rar_FDTD matlab_FDTD matlab 电磁_fdtd_fdtd 脉冲_高斯脉冲

% 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/tau^2), % 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:\MATLAB7\work\cavity\figures') % mkdir 'C:\MATLAB7\work\cavity\figures' % 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 SimRunSt