FDTD matlab 程序

%*********************************************************************** % 2-D FDTD TE code with PML absorbing boundary conditions %*********************************************************************** % % Program author: Susan C. Hagness % Department of Electrical and Computer Engineering % University of Wisconsin-Madison % 1415 Engineering Drive % Madison, WI 53706-1691 % 608-265-5739 % hagness@engr.wisc.edu % % Date of this version: February 2000 % % This MATLAB M-file implements the finite-difference time-domain % solution of Maxwell's curl equations over a two-dimensional % Cartesian space lattice comprised of uniform square grid cells. % % To illustrate the algorithm, a 6-cm-diameter metal cylindrical % scatterer in free space is modeled. The source excitation is % a Gaussian pulse with a carrier frequency of 5 GHz. % % The grid resolution (dx = 3 mm) was chosen to provide 20 samples % per wavelength at the center frequency of the pulse (which in turn % provides approximately 10 samples per wavelength at the high end % of the excitation spectrum, around 10 GHz). % % The computational domain is truncated using the perfectly matched % layer (PML) absorbing boundary conditions. The formulation used % in this code is based on the original split-field Berenger PML. The % PML regions are labeled as shown in the following diagram: % % ---------------------------------------------- % | | BACK PML | | % ---------------------------------------------- % |L | /| R| % |E | (ib,jb) | I| % |F | | G| % |T | | H| % | | MAIN GRID | T| % |P | | | % |M | | P| % |L | (1,1) | M| % | |/ | L| % ---------------------------------------------- % | | FRONT PML | | % ---------------------------------------------- % % To execute this M-file, type "fdtd2D" at the MATLAB prompt. % This M-file displays the FDTD-computed Ex, Ey, and Hz fields at % every 4th time step, and records those frames in a movie matrix, % M, which is played at the end of the simulation using the "movie" % command. % %*********************************************************************** clear %*********************************************************************** % Fundamental constants %*********************************************************************** cc=2.99792458e8; %speed of light in free space muz=4.0*pi*1.0e-7; %permeability of free space epsz=1.0/(cc*cc*muz); %permittivity of free space freq=5.0e+9; %center frequency of source excitation lambda=cc/freq; %center wavelength of source excitation omega=2.0*pi*freq; %*********************************************************************** % Grid parameters %*********************************************************************** ie=100; %number of grid cells in x-direction je=50; %number of grid cells in y-direction ib=ie+1; jb=je+1; is=15; %location of z-directed hard source js=je/2; %location of z-directed hard source dx=3.0e-3; %space increment of square lattice dt=dx/(2.0*cc); %time step nmax=300; %total number of time steps iebc=8; %thickness of left and right PML region jebc=8; %thickness of front and back PML region rmax=0.00001; orderbc=2; ibbc=iebc+1; jbbc=jebc+1; iefbc=ie+2*iebc; jefbc=je+2*jebc; ibfbc=iefbc+1; jbfbc=jefbc+1; %*********************************************************************** % Material parameters %*********************************************************************** media=2; eps=[1.0 1.0]; sig=[0.0 1.0e+7]; mur=[1.0 1.0]; sim=[0.0 0.0]; %*********************************************************************** % Wave excitation %*********************************************************************** rtau=160.0e-12; tau=rtau/dt; delay=3*tau; source=zeros(1,nmax); for n=1:7.0*tau source(n)=sin(omega*(n-delay)*dt)*exp(-((n-delay)^2/tau^2)); end %*********************************************************************** % Field arrays %*********************************************************************** ex=zeros(ie,jb); %fields in main grid ey=zeros(ib,je); hz=zeros(ie,je); exbcf=zeros(iefbc,jebc); %fields in front PML region eybcf=zeros(ibfbc,jebc); hzxbcf=zeros(iefbc,jebc); hzybcf=zeros(iefbc,jebc); exbcb=zeros(iefbc,jbbc); %fields in back PML region eybcb=zeros(ibfbc,jebc); hzxbcb=zeros(iefbc,jebc); hzybcb=zeros(iefbc,jebc); exbcl=zeros(iebc,jb); %fields in left PML region eybcl=zeros(iebc,je); hzxbcl=zeros(iebc,je); hzybcl=zeros(iebc,je); exbcr=zeros(iebc,jb); %fields in right PML region eybcr=zeros(ibbc,je); hzxbcr=zeros(iebc,je); hzybcr=zeros(iebc,je); %*********************************************************************** % Updating coefficients %*********************************************************************** for i=1:media eaf =dt*sig(i)/(2.0*epsz*eps(i)); ca(i)=(1.0-eaf)/(1.0+eaf); cb(i)=dt/epsz/eps(i)/dx/(1.0+eaf); haf =dt*sim(i)/(2.0*muz*mur(i)); da(i)=(1.0-haf)/(1.0+haf); db(i)=dt/muz/mur(i)/dx/(1.0+haf); end %*********************************************************************** % Geometry specification (main grid) %*********************************************************************** % Initialize entire main grid to free space caex(1:ie,1:jb)=ca(1); cbex(1:ie,1:jb)=cb(1); caey(1:ib,1:je)=ca(1); cbey(1:ib,1:je)=cb(1); dahz(1:ie,1:je)=da(1); dbhz(1:ie,1:je)=db(1); % Add metal cylinder diam=20; % diameter of cylinder: 6 cm rad=diam/2.0; % radius of cylinder: 3 cm icenter=4*ie/5; % i-coordinate of cylinder's center jcenter=je/2; % j-coordinate of cylinder's center for i=1:ie for j=1:je dist2=(i+0.5-icenter)^2 + (j-jcenter)^2; if dist2 <= rad^2 caex(i,j)=ca(2); cbex(i,j)=cb(2); end dist2=(i-icenter)^2 + (j+0.5-jcenter)^2; if dist2 <= rad^2 caey(i,j)=ca(2); cbey(i,j)=cb(2); end end end %*********************************************************************** % Fill the PML regions %*********************************************************************** delbc=iebc*dx; sigmam=-log(rmax/100.0)*epsz*cc*(orderbc+1)/(2*delbc); bcfactor=eps(1)*sigmam/(dx*(delbc^orderbc)*(orderbc+1)); % FRONT region caexbcf(1:iefbc,1)=1.0; cbexbcf(1:iefbc,1)=0.0; for j=2:jebc y1=(jebc-j+1.5)*dx; y2=(jebc-j+0.5)*dx; sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1)); ca1=exp(-sigmay*dt/(epsz*eps(1))); cb1=(1.0-ca1)/(sigmay*dx); caexbcf(1:iefbc,j)=ca1; cbexbcf(1:iefbc,j)=cb1; end sigmay = bcfactor*(0.5*dx)^(orderbc+1); ca1=exp(-sigmay*dt/(epsz*eps(1))); cb1=(1-ca1)/(sigmay*dx); caex(1:ie,1)=ca1; cbex(1:ie,1)=cb1; caexbcl(1:iebc,1)=ca1; cbexbcl(1:iebc,1)=cb1; caexbcr(1:iebc,1)=ca1; cbexbcr(1:iebc,1)=cb1; for j=1:jebc y1=(jebc-j+1)*dx; y2=(jebc-j)*dx; sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1)); sigmays=sigmay*(muz/(epsz*eps(1))); da1=exp(-sigmays*dt/muz); db1=(1-da1)/(sigmays*dx); dahzybcf(1:iefbc