FDTD算法的Matlab程序

1、精选优质文档-倾情为你奉上%* 5= T$h;O  %    3-D FDTD code with PEC boundaries Iw?*y.z|  %* H!'4A&  % r dCs  %    Program author: Susan C. Hagness B_?7+  %                    Department of Electrical and Compute

2、r Engineering &hrMpD6z6i  %                    University of Wisconsin-Madison T/|nOu 5  %                    1415 Engineering Drive 41P0)o  %          

3、0;         Madison, WI 53706-1691 9MH;=88q  %                    608-265-5739 :s zkh?  %                    |UA=? Xl  % )_;l%&  %    Date of this

4、version:  February 2000 !$Ivro  % yYSmmgrX0  %    This MATLAB M-file implements the finite-difference time-domain ("F$r$9S  %    solution of Maxwell's curl equations over a three-dimensional $gN%X/n"1  %    Cartesian space lattice

5、comprised of uniform cubic grid cells. %hN(79:g  % vq.o ;q /  %    To illustrate the algorithm, an air-filled rectangular cavity < Np Mv!g  %    resonator is modeled.  The length, width, and height of the pojhfn  %    cavity are 10.0 cm (x-di

6、rection), 4.8 cm (y-direction), and VcXr!4 M  %    2.0 cm (z-direction), respectively. DU81(  % #pgD-0_  %    The computational domain is truncated using PEC boundary >Gy+-  %    conditions: yUDoOVC0  %          ex(i,

7、j,k)=0 on the j=1, j=jb, k=1, and k=kb planes kyN)  %          ey(i,j,k)=0 on the i=1, i=ib, k=1, and k=kb planes 7Hv 6>z#m  %          ez(i,j,k)=0 on the i=1, i=ib, j=1, and j=jb planes <9"<,  %    These PEC bo

8、undaries form the outer lossless walls of the cavity. A6z.MdYZ  % a&qdp  %    The cavity is excited by an additive current source oriented x 9 a  %    along the z-direction.  The source waveform is a differentiated Qj1%'wWG  %    Gaussia

9、n pulse given by #K!jh)y  %          J(t)=-J0*(t-t0)*exp(-(t-t0)2/tau2), dMjeF  %    where tau=50 ps.  The FWHM spectral bandwidth of this zero-dc- eYoc(bG(+  %    content pulse is approximately 7 GHz. The grid resolution Gs|a$ V|o

10、60;%    (dx = 2 mm) was chosen to provide at least 10 samples per j1;0kb?  %    wavelength up through 15 GHz. jqjj2 9  % _X<V , p  %    To execute this M-file, type "fdtd3D" at the MATLAB prompt. g(sR ?  %    This M-file displa

11、ys the FDTD-computed Ez fields at every other Wt!;Y,1 s  %    time step, and records those frames in a movie matrix, M, which i WGGnGS  %    is played at the end of the simulation using the "movie" command. WTt /y'6  % j+/EG*/  %* cMy?&  

12、;_3f/lG?&-  clear EOqV5$+  czWUD  %* ov!L8 9u  %    Fundamental constants Tgq,tR  %* <_bGV  d +eb!fi  cc=2.e8;            %speed of light in free space Di_K  muz=4.0*pi*1.0e-7;          %p

13、ermeability of free space <lPHeO<  epsz=1.0/(cc*cc*muz);      %permittivity of free space c/3M>+M  ?nR$>a  %* $:DhK  %    Grid parameters .7#04_aP  %* jRjQDK_"ka  (<G) +*  ie=50;      %number of grid cells

14、 in x-direction pQ+4+7ID  je=24;      %number of grid cells in y-direction L30>?l  ke=10;      %number of grid cells in z-direction Aa%ks+1  s0'Xihsw6  ib=ie+1; gB#$"mq,  jb=je+1; D#&N?<  kb=ke+1; iW <B1'dp  Y

15、JrZ  is=26;      %location of z-directed current source +xtRY"  js=13;      %location of z-directed current source LaiUf_W#X  X<"W  kobs=5; <El6?ml  47XQZ-4  dx=0.002;          %space increment of cu

16、bic lattice jZ3  dt=dx/(2.0*cc);    %time step N XzgI  SuW_6  nmax=500;          %total number of time steps # 1,(I  ghiFI<)VY  cF aYI  %* /Et:',D  %    Differentiated Gaussian pulse excitation >ACMO  %* d

17、YJWQ;j.|  F/PN1#T  rtau=50.0e-12; |k> _ jO  tau=rtau/dt; y_8 8I:O  ndelay=3*tau; /W/ =OPe  srcconst=-dt*3.0e+11; mg*,_3q33  6546"sU  %* %Sfew/"R0  %    Material parameters qI" PI!s  %* C jsy1gA  FU Ip  eps=1.0; 3xhxE

18、 sig=0.0; 5R/!e(m  r+;op_  %* :v&TQ  %    Updating coefficients iu2%S)w  %* SS1-UbL  q<yp6Q3  ca=(1.0-(dt*sig)/(2.0*epsz*eps)/(1.0+(dt*sig)/(2.0*epsz*eps); AF>J8V  cb=(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps); ?k?7GN  da=1.0; +dcB

19、h Dq  db=dt/muz/dx; #TM+Vd$  >f JY  %* -uXf?sTV  %    Field arrays ra/5D  %* $014/IB  fb.jXSR  ex=zeros(ie,jb,kb); eW, E)x:  ey=zeros(ib,je,kb); bfcQ( m5  ez=zeros(ib,jb,ke); ul$k xc=N  hx=zeros(ib,je,ke); W%5o87  hy=zeros(ie,jb,ke)

20、; q#B=PZ'NA  hz=zeros(ie,je,kb); Z;Q2tT /F  OSreS5bg  %* +Ch2Lod  %    Movie initialization C.qN Bl*  %* .n$c+  O0m_  tview(:,:)=ez(:,:,kobs); elF#$  sview(:,:)=ez(:,js,:); B xAyjA6  |5il5UP  subplot('position',0.15 0.45 0.7 0.4

21、5),pcolor(tview'); +_EE  shading flat; mH:8_=(.  *is(-1.0 1.0); 6roq 1=  colorbar; ,GeW_!Q  axis image; Ms<v81z5T  title('Ez(i,j,k=5), time step = 0'); Jb/VITqN4  xlabel('i coordinate'); pTB1I3=.u  ylabel('j coordinate'); c1tM(& &

22、#160;5J-slNNCQ  subplot('position',0.15 0.10 0.7 0.25),pcolor(sview'); Mz&h>  shading flat; :3N6Ej  *is(-1.0 1.0); ajCe&+  colorbar; %E8HLTEvl  axis image; E '|  title('Ez(i,j=13,k), time step = 0'); qWU59:d  xlabel('i coordin

23、ate'); 3qQU-;|  ylabel('k coordinate'); zr 4JTS  9=J 3T66U  rect=get(gcf,'Position'); qEH  rect(1:2)=0 0; # uyAC$  )j QrD  M=moviein(nmax/2,gcf,rect); DUmp6  X2zIFm  %* _gfec4o  %    BEGIN TIME-STEPPING LOOP Z= LLL  %* u

24、F ;8B"  P_Pco  for n=1:nmax rM4RibS  5Z(q|nn7P  %* uG+eF  %    Update electric fields pf#R  %* + jN)$Y3Ya  x f)P  ex(1:ie,2:je,2:ke)=ca*ex(1:ie,2:je,2:ke)+. >J,IxRGi                    cb*(hz(

25、1:ie,2:je,2:ke)-hz(1:ie,1:je-1,2:ke)+. Pk8(2fAYk                        hy(1:ie,2:je,1:ke-1)-hy(1:ie,2:je,2:ke); ()fYhk|W  >?> !#1  ey(2:ie,1:je,2:ke)=ca*ey(2:ie,1:je,2:ke)+. C?rb(m            &

26、#160;       cb*(hx(2:ie,1:je,2:ke)-hx(2:ie,1:je,1:ke-1)+. 1 y7$"N8Xo                        hz(1:ie-1,1:je,2:ke)-hz(2:ie,1:je,2:ke); d 8z9_C-  IQB%v5  ez(2:ie,2:je,1:ke)=ca*ez(2:ie,2:je,1:ke)+. fM<g+X  

27、                  cb*(hx(2:ie,1:je-1,1:ke)-hx(2:ie,2:je,1:ke)+. r+n hm"9                        hy(2:ie,2:je,1:ke)-hy(1:ie-1,2:je,1:ke); <h_7Dn  rgu7g  ez(is,js,1:ke)=ez(is,js,1

28、:ke)+. pJETM                srcconst*(n-ndelay)*exp(-(n-ndelay)2/tau2); -PH qD  j"<F?kQ  %* V7>,  %    Update magnetic fields o(5 ( bJ  %* o n?8l?iQ  ?;Ge/QU5  hx(2:ie,1:je,1:ke)=hx(2:ie,1:je,1:ke)+. ,:/3'

29、L                    db*(ey(2:ie,1:je,2:kb)-ey(2:ie,1:je,1:ke)+. .Ue1'v*,                        ez(2:ie,1:je,1:ke)-ez(2:ie,2:jb,1:ke); 2o-Ie/"d  m_1BB$lyP2  hy(1:ie,2:je,

30、1:ke)=hy(1:ie,2:je,1:ke)+. AqPQeNgz                    db*(ex(1:ie,2:je,1:ke)-ex(1:ie,2:je,2:kb)+. u:D,;)                        ez(2:ib,2:je,1:ke)-ez(1:ie,2:je,1:ke); Sf*b6lcC  WT>2

31、eMK  hz(1:ie,1:je,2:ke)=hz(1:ie,1:je,2:ke)+. QWV12t$v                    db*(ex(1:ie,2:jb,2:ke)-ex(1:ie,1:je,2:ke)+. z/Yrf                        ey(1:ie,1:je,2:ke)-ey(2:ib,1:je,2:ke);

32、=%V(n7=  F"Y.'my8  %* 6vQCghI  %    Visualize fields gK8=A0c  %* A7qKY-4B  HYf  if mod(n,2)=0; $Dm2>:Dmt  plRBfw>N  timestep=int2str(n); BB694  tview(:,:)=ez(:,:,kobs); :d ts>  sview(:,:)=ez(:,js,:); *+qlam4N  >vPDF+u