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

%* =a)6y % 3-D FDTD code with PEC boundaries rd5xlk %* G)hH?_U#T % &5yh % Program author: Susan C. Hagness UD2%8L % Madison, WI 53706-1691 (w4.! % 608-265-5739 B m:JOx % v* %x % Date of this version: February 2000 ypF 8V 8 % d50IAap6J % This MATLAB M-file implements the finite-difference time-domain 1!.V* % solution of Maxwells curl equations over a three-dimensional $qtU % Cartesian space lattice comprised of uniform cubic grid cells. ;eEGY % !| f%UO % To illustrate the algorithm, an air-filled rectangular cavity I0l.KiBm % resonator is modeled. The length, width, and height of the ay:P.5) % cavity are 10.0 cm (x-direction), 4.8 cm (y-direction), and Tc3ihLvG % 2.0 cm (z-direction), respectively. b*FU*)4. % Xf:E1J % The computational domain is truncated using PEC boundary p, 3A:i % conditions: GwMUIevO_ % ex(i,j,k)=0 on the j=1, j=jb, k=1, and k=kb planes O+c*|w % ey(i,j,k)=0 on the i=1, i=ib, k=1, and k=kb planes hHfe6P | % ez(i,j,k)=0 on the i=1, i=ib, j=1, and j=jb planes =p,4=wo % These PEC boundaries form the outer lossless walls of the cavity. rY3_NG% % ,;lv( % The cavity is excited by an additive current source oriented (:&;sI % along the z-direction. The source waveform is a differentiated q#Yg0w % Gaussian pulse given by HC9vc,Fp % J(t)=-J0*(t-t0)*exp(-(t-t0)2/tau2), zH=/.31Q % where tau=50 ps. The FWHM spectral bandwidth of this zero-dc- uS: A4tN % content pulse is approximately 7 GHz. The grid resolution F1?CqN M % (dx = 2 mm) was chosen to provide at least 10 samples per z&C8aQ % wavelength up through 15 GHz. a!c/5)v( % l#2r.q$| % To execute this M-file, type fdtd3D at the MATLAB prompt. Xmf % This M-file displays the FDTD-computed Ez fields at every other 3JlC/v#0 % time step, and records those frames in a movie matrix, M, which ;H7EB % is played at the end of the simulation using the movie command. la4L4o % Qag|nLoT %* a-T*F FSW3 clear %-h7Z3YcN &q , !:L %* Bl kSWW/ % Fundamental constants E0sbU11 %* tpv?(DDU pd=7; cc=2.99792458e8; %speed of light in free space O?U muz=4.0*pi*1.0e-7; %permeability of free space aj&CJ epsz=1.0/(cc*cc*muz); %permittivity of free space MB :GY? Hm.XHO0L %* jFgZXp % Grid parameters Vj3yK %* =Ov9Kf V78Mq:7d ie=50; %number of grid cells in x-direction Xfj)gPt je=24; %number of grid cells in y-direction OlIT|bzkb ke=10; %number of grid cells in z-direction H!p!sn yq1c*| ib=ie+1; BeAk 21xb jb=je+1; 8a7YHULv|d T=.-Cl1A nmax=500; %total number of time steps 5|5=Y/ WBa /IM S3QaYqv %* hc6u % Differentiated Gaussian pulse excitation hw*1gm %* =Ajj-r& =88t*dH(, rtau=50.0e-12; +3) 1Z tau=rtau/dt; moaodmtx ndelay=3*tau; &lUNy L srcconst=-dt*3.0e+11; ;&1A RvR:e| %* Wtf % Material parameters ex!XB$X %* UDM yyVd %6Rn 4J eps=1.0; h07eE g sig=0.0; qfa3k8et 4T%cTH:.9N %* 2uEhOi0I % Updating coefficients G0pBR_5z$ %* VyY.r# M7vjmt? ca=(1.0-(dt*sig)/(2.0*epsz*eps)/(1.0+(dt*sig)/(2.0*epsz*eps); ZtR&wk cb=(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps); R0Y+#$8h da=1.0; rn H#u+ db=dt/muz/dx; wFb1ae -bm,:Iy! %* B=dseeGTo % Field arrays lJ&y&NO %* L+t&1cW HCazwX ex=zeros(ie,jb,kb); K$ AB Fvc ey=zeros(ib,je,kb); |xeE3,8 ez=zeros(ib,jb,ke); mW +tV1XjG hx=zeros(ib,je,ke); UJn/s;$.e hy=zeros(ie,jb,ke); SX_4= hz=zeros(ie,je,kb); c&m9)rzP o%qkqK1 %* Q1ToxV % Movie initialization Bm iU(Z %* v|&s4x?D /?Jylj tview(:,:)=ez(:,:,kobs); mf2Mx=oy sview(:,:)=ez(:,js,:); LBzsIv ONP0E subplot(position,0.15 0.45 0.7 0.45),pcolor(tview); 2kJ!En7 shading flat; M2pR#9 xlabel(i coordinate); .YC;zn ylabel(j coordinate); IgS u8,!) 4i subplot(position,0.15 0.10 0.7 0.25),pcolor(sview); T.2ZBG | shading flat; mG/(LX8 caxis(-1.0 1.0); !7c+Hm colorbar; M)bHB(v axis image; s|gp title(Ez(i,j=13,k), time step = 0); *2R.u xlabel(i coordinate); !oICHHH ylabel(k coordinate); 9=2g2Q 738Z/) rect=get(gcf,Position); yZb+=UM rect(1:2)=0 0; RM#fX)= O% VA) %* jH0SBKC % BEGIN TIME-STEPPING LOOP Wt4 %* cyHU!Z*Zq 5Edo%Hd6 for n=1:nmax N UFIjWh %* sNf +lga0 % Update electric fields Y!it!9 %* 1xIFvXru IW mHp ex(1:ie,2:je,2:ke)=ca*ex(1:ie,2:je,2:ke)+. Q%Jc( cb*(hz(1:ie,2:je,2:ke)-hz(1:ie,1:je-1,2:ke)+. d9| T=R hy(1:ie,2:je,1:ke-1)-hy(1:ie,2:je,2:ke); 97(nWt 2 r5 ey(2:ie,1:je,2:ke)=ca*ey(2:ie,1:je,2:ke)+. | 4 .#4 cb*(hx(2:ie,1:je,2:ke)-hx(2:ie,1:je,1:ke-1)+. n8F7 hz(1:ie-1,1:je,2:ke)-hz(2:ie,1:je,2:ke); Xa)7bp &%bRPUl ez(2:ie,2:je,1:ke)=ca*ez(2:ie,2:je,1:ke)+. N:d D*QZ cb*(hx(2:ie,1:je-1,1:ke)-hx(2:ie,2:je,1:ke)+. dbg|V oNf hy(2:ie,2:je,1:ke)-hy(1:ie-1,2:je,1:ke); _m lo +HQP8 ez(is,js,1:ke)=ez(is,js,1:ke)+. +B(5z4 srcconst*(n-ndelay)*exp(-(n-ndelay)2/tau2); cK/PQsMP 9*f2b.Aj %* #%;FFu % Update magnetic fields 2n _T2 %* 1k*n1t): x_,2 W hx(2:ie,1:je,1:ke)=hx(2:ie,1:je,1:ke)+. _: / A db*(ey(2:ie,1:je,2:kb)-ey(2:ie,1:je,1:ke)+. e=0lRj ez(2:ie,1:je,1:ke)-ez(2:ie,2:jb,1:ke); u XaL K*g, s+ hy(1:ie,2:je,1:ke)=hy(1:ie,2:je,1:ke)+. CaOgt/U db*(ex(1:ie,2:je,1:ke)-ex(1:ie,2:je,2:kb)+. Gr3* ez(2:ib,2:je,1:ke)-ez(1:ie,2:je,1:ke); o9OPoof:. Vp5V m hz(1:ie,1:je,2:ke)=hz(1:ie,1:je,2:ke)+. Ah(%35& db*(ex(1:ie,2:jb,2:ke)-ex(1:ie,1:je,2:ke)+. C5.;7& ey(1:ie,1:je,2:ke)-ey(2:ib,1:je,2:ke); AHg gwT bCJ&6O4 if mod(n,2)=0; 7cOg(6N pw(U: w title(Ez(i,j,k=5), time step = ,timestep); tl,xp xlabel(i coordinate); x :Bw; ylabel(j coordinate); ,NIcl r6B+3J subplot(position,0.15 0.10 0.7 0.25),pcolor(sview); Z+h7 0,| shading flat; PE6u8ZAb caxis(-1.0 1.0); -tH&X colorbar; Rcxa:k axis image; Rp%+Xz title(Ez(i,j=13,k), time step = ,ti