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

%* 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 % Department of Electrical and Computer Engineering color=#FFFFFF本文来自:/color % University of Wisconsin-Madison color=#FFFFFF本文来自:/color % 1415 Engineering Drive color=#FFFFFF本文来自:/color % Madison, WI 53706-1691 color=#FFFFFF本文来自:/color % 608-265-5739 color=#FFFFFF本文来自:/color % color=#FFFFFF本文来自:/color % color=#FFFFFF本文来自:/color % Date of this version: February 2000 color=#FFFFFF本文来自:/color % color=#FFFFFF本文来自:/color % This MATLAB M-file implements the finite-difference time-domain color=#FFFFFF本文来自:/color % solution of Maxwells curl equations over a three-dimensional color=#FFFFFF本文来自:/color % Cartesian space lattice comprised of uniform cubic grid cells. color=#FFFFFF本文来自:/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.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 % 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 % 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 additive 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(-(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 chosen 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 % 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 command. 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 %* 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 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; %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 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 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 color=#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 excitation 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=#FFFFFF本文来自:/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 %* 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*epsz*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本文来自:/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本文来自:/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=#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; 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 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本文来自:/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本文来自:/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=#FFFFFF本文来自:/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 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本文来自:/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 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=#FFFFFF本文来自:/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 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: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: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 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(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 % Visualize 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); color=#FFFFFF本文来自:/color shading flat; color=#FFFFFF本文来自:/color caxis(-1.0 1.0); color=#FFFFFF本文来自:/color colorbar; color=#FFFFFF本文来自:/color axis image; color=#FF