時域有限差分法算法的Matlab源程序

1、%* color=#FFFFFF本文來自:/color  %    3-D FDTD code with PEC boundaries color=#FFFFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  % color=#FFFFFF本文來自:/color  %    Program author: Susan C. Hagness color=#FFFFFF本文來自:/color  %          

2、0;         Department of Electrical and Computer Engineering color=#FFFFFF本文來自:/color  %                    University of Wisconsin-Madison color=#FFFFFF本文來自:/color  %              

3、0;     1415 Engineering Drive color=#FFFFFF本文來自:/color  %                    Madison, WI 53706-1691 color=#FFFFFF本文來自:/color  %                    608-265-5739 color=#FFFFFF本文來自:/c

4、olor  %                    color=#FFFFFF本文來自:/color  % color=#FFFFFF本文來自:/color  %    Date of this version:  February 2000 color=#FFFFFF本文來自:/color  % color=#FFFFFF本文來自:/color  %  

5、0; This MATLAB M-file implements the finite-difference time-domain color=#FFFFFF本文來自:/color  %    solution of Maxwell's curl equations over a three-dimensional color=#FFFFFF本文來自:/color  %    Cartesian space lattice comprised of uniform cubic grid cells. color=#FFFFF

6、F本文來自:/color  % color=#FFFFFF本文來自:/color  %    To illustrate the algorithm, an air-filled rectangular cavity color=#FFFFFF本文來自:/color  %    resonator is modeled.  The length, width, and height of the color=#FFFFFF本文來自:/color  %    cavity are 10.

7、0 cm (x-direction), 4.8 cm (y-direction), and color=#FFFFFF本文來自:/color  %    2.0 cm (z-direction), respectively. color=#FFFFFF本文來自:/color  % color=#FFFFFF本文來自:/color  %    The computational domain is truncated using PEC boundary color=#FFFFFF本文來自:/color  %

8、0;   conditions: color=#FFFFFF本文來自:/color  %          ex(i,j,k)=0 on the j=1, j=jb, k=1, and k=kb planes color=#FFFFFF本文來自:/color  %          ey(i,j,k)=0 on the i=1, i=ib, k=1, and k=kb planes color=#FFFFFF本文來自:/color  % 

9、60;       ez(i,j,k)=0 on the i=1, i=ib, j=1, and j=jb planes color=#FFFFFF本文來自:/color  %    These PEC boundaries form the outer lossless walls of the cavity. color=#FFFFFF本文來自:/color  % color=#FFFFFF本文來自:/color  %    The cavity is excited by an addit

10、ive current source oriented color=#FFFFFF本文來自:/color  %    along the z-direction.  The source waveform is a differentiated color=#FFFFFF本文來自:/color  %    Gaussian pulse given by color=#FFFFFF本文來自:/color  %          J(t)=-J0*(t-t0)*exp(

11、-(t-t0)2/tau2), color=#FFFFFF本文來自:/color  %    where tau=50 ps.  The FWHM spectral bandwidth of this zero-dc- color=#FFFFFF本文來自:/color  %    content pulse is approximately 7 GHz. The grid resolution color=#FFFFFF本文來自:/color  %    (dx = 2 mm) was chos

12、en to provide at least 10 samples per color=#FFFFFF本文來自:/color  %    wavelength up through 15 GHz. color=#FFFFFF本文來自:/color  % color=#FFFFFF本文來自:/color  %    To execute this M-file, type "fdtd3D" at the MATLAB prompt. color=#FFFFFF本文來自:/color  % 

13、;   This M-file displays the FDTD-computed Ez fields at every other color=#FFFFFF本文來自:/color  %    time step, and records those frames in a movie matrix, M, which color=#FFFFFF本文來自:/color  %    is played at the end of the simulation using the "movie" comm

14、and. color=#FFFFFF本文來自:/color  % color=#FFFFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  clear color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  %    Fundamental constants color=#FFFFFF本文來自:/color &#

15、160;%* color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  cc=2.99792458e8;            %speed of light in free space color=#FFFFFF本文來自:/color  muz=4.0*pi*1.0e-7;          %permeability of free space color=#FFFFFF本文來自:/color  

16、epsz=1.0/(cc*cc*muz);      %permittivity of free space color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  %    Grid parameters color=#FFFFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  ie=50;&

17、#160;     %number of grid cells in x-direction color=#FFFFFF本文來自:/color  je=24;      %number of grid cells in y-direction color=#FFFFFF本文來自:/color  ke=10;    %number of grid cells in z-direction color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color 

18、0;ib=ie+1; color=#FFFFFF本文來自:/color  jb=je+1; color=#FFFFFF本文來自:/color  kb=ke+1; color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  is=26;      %location of z-directed current source color=#FFFFFF本文來自:/color  js=13;      %location of z-directed

19、 current source color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  kobs=5; color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  dx=0.002;          %space increment of cubic lattice color=#FFFFFF本文來自:/color  dt=dx/(2.0*cc);    %time step col

20、or=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  nmax=500;        %total number of time steps color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  %    Differentiated Gaussian pulse exci

21、tation color=#FFFFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  rtau=50.0e-12; color=#FFFFFF本文來自:/color  tau=rtau/dt; color=#FFFFFF本文來自:/color  ndelay=3*tau; color=#FFFFFF本文來自:/color  srcconst=-dt*3.0e+11; color=#FFFFFF本文來自:/color  color=#FF

22、FFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  %    Material parameters color=#FFFFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  eps=1.0; color=#FFFFFF本文來自:/color  sig=0.0; color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  %*

23、 color=#FFFFFF本文來自:/color  %    Updating coefficients color=#FFFFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  ca=(1.0-(dt*sig)/(2.0*epsz*eps)/(1.0+(dt*sig)/(2.0*epsz*eps); color=#FFFFFF本文來自:/color  cb=(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*eps

24、z*eps); color=#FFFFFF本文來自:/color  da=1.0; color=#FFFFFF本文來自:/color  db=dt/muz/dx; color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  %    Field arrays color=#FFFFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/

25、color  ex=zeros(ie,jb,kb); color=#FFFFFF本文來自:/color  ey=zeros(ib,je,kb); color=#FFFFFF本文來自:/color  ez=zeros(ib,jb,ke); color=#FFFFFF本文來自:/color  hx=zeros(ib,je,ke); color=#FFFFFF本文來自:/color  hy=zeros(ie,jb,ke); color=#FFFFFF本文來自:/color  hz=zeros(ie,je,kb); color=#FFFFFF

26、本文來自:/color  color=#FFFFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  %    Movie initialization color=#FFFFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  tview(:,:)=ez(:,:,kobs); color=#FFFFFF本文來自:/color  sview(:,:)=ez(:,js,:); color=

27、#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  subplot('position',0.15 0.45 0.7 0.45),pcolor(tview'); color=#FFFFFF本文來自:/color  shading flat; color=#FFFFFF本文來自:/color  caxis(-1.0 1.0); color=#FFFFFF本文來自:/color  colorbar; color=#FFFFFF本文來自:/color  axis image;

28、color=#FFFFFF本文來自:/color  title('Ez(i,j,k=5), time step = 0'); color=#FFFFFF本文來自:/color  xlabel('i coordinate'); color=#FFFFFF本文來自:/color  ylabel('j coordinate'); color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  subplot('position',0.15 0.10

29、 0.7 0.25),pcolor(sview'); color=#FFFFFF本文來自:/color  shading flat; color=#FFFFFF本文來自:/color  caxis(-1.0 1.0); color=#FFFFFF本文來自:/color  colorbar; color=#FFFFFF本文來自:/color  axis image; color=#FFFFFF本文來自:/color  title('Ez(i,j=13,k), time step = 0'); color=#FFFFFF本文

30、來自:/color  xlabel('i coordinate'); color=#FFFFFF本文來自:/color  ylabel('k coordinate'); color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  rect=get(gcf,'Position'); color=#FFFFFF本文來自:/color  rect(1:2)=0 0; color=#FFFFFF本文來自:/color  color=#FFFFFF本文來

31、自:/color  M=moviein(nmax/2,gcf,rect); color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  %    BEGIN TIME-STEPPING LOOP color=#FFFFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  for n=1:nmax color=#FFFFF

32、F本文來自:/color  color=#FFFFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  %    Update electric fields color=#FFFFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  ex(1:ie,2:je,2:ke)=ca*ex(1:ie,2:je,2:ke)+. color=#FFFFFF本文來自:/color    &

33、#160;               cb*(hz(1:ie,2:je,2:ke)-hz(1:ie,1:je-1,2:ke)+. color=#FFFFFF本文來自:/color                        hy(1:ie,2:je,1:ke-1)-hy(1:ie,2:je,2:ke); color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/

34、color  ey(2:ie,1:je,2:ke)=ca*ey(2:ie,1:je,2:ke)+. color=#FFFFFF本文來自:/color                    cb*(hx(2:ie,1:je,2:ke)-hx(2:ie,1:je,1:ke-1)+. color=#FFFFFF本文來自:/color                      

35、0; hz(1:ie-1,1:je,2:ke)-hz(2:ie,1:je,2:ke); color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  ez(2:ie,2:je,1:ke)=ca*ez(2:ie,2:je,1:ke)+. color=#FFFFFF本文來自:/color                    cb*(hx(2:ie,1:je-1,1:ke)-hx(2:ie,2:je,1:ke)+. color=#FFFF

36、FF本文來自:/color                        hy(2:ie,2:je,1:ke)-hy(1:ie-1,2:je,1:ke); color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  ez(is,js,1:ke)=ez(is,js,1:ke)+. color=#FFFFFF本文來自:/color            &#

37、160;   srcconst*(n-ndelay)*exp(-(n-ndelay)2/tau2); color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  %    Update magnetic fields color=#FFFFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  hx(2:ie,1:je,1

38、:ke)=hx(2:ie,1:je,1:ke)+. color=#FFFFFF本文來自:/color                    db*(ey(2:ie,1:je,2:kb)-ey(2:ie,1:je,1:ke)+. color=#FFFFFF本文來自:/color                        ez(2:ie,1:je,1:ke)-ez(2:ie,2

39、:jb,1:ke); color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  hy(1:ie,2:je,1:ke)=hy(1:ie,2:je,1:ke)+. color=#FFFFFF本文來自:/color                    db*(ex(1:ie,2:je,1:ke)-ex(1:ie,2:je,2:kb)+. color=#FFFFFF本文來自:/color      

40、0;                 ez(2:ib,2:je,1:ke)-ez(1:ie,2:je,1:ke); color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  hz(1:ie,1:je,2:ke)=hz(1:ie,1:je,2:ke)+. color=#FFFFFF本文來自:/color                    db*(ex

41、(1:ie,2:jb,2:ke)-ex(1:ie,1:je,2:ke)+. color=#FFFFFF本文來自:/color                        ey(1:ie,1:je,2:ke)-ey(2:ib,1:je,2:ke); color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  %    Visu

42、alize fields color=#FFFFFF本文來自:/color  %* color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  if mod(n,2)=0; color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  timestep=int2str(n); color=#FFFFFF本文來自:/color  tview(:,:)=ez(:,:,kobs); color=#FFFFFF本文來自:/color  sview(:,:)=ez(:,js,:); color=#FFFFFF本文來自:/color  color=#FFFFFF本文來自:/color  subplot('position',0.15 0.45 0.7 0.45),pcolor(tview'); col