3 FDTD学习代码.zip_FDTD file_fdtd_it_waveguide_waveguide FDTD

%*********************************************************************** % 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 % hagness@engr.wisc.edu % % Copyright 2005 % % 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. % Exponential time stepping is implemented in the PML regions. % 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 code has been tested in the following Matlab environments: % Matlab version 6.1.0.450 Release 12.1 (May 18, 2001) % Matlab version 6.5.1.199709 Release 13 Service Pack 1 (August 4, 2003) % Matlab version 7.0.0.19920 R14 (May 6, 2004) % Matlab version 7.0.1.24704 R14 Service Pack 1 (September 13, 2004) % Matlab version 7.0.4.365 R14 Service Pack 2 (January 29, 2005) % % Note: if you are using Matlab version 6.x, you may wish to make % one or more of the following modifications to this code: % --uncomment line numbers 418 and 419 % --comment out line numbers 610, 619, and 628 % %*********************************************************************** 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 etaz=sqrt(muz/epsz); 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=1.0e-7; 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 % % PML theory describes a continuous grading of the material properties % over the PML region. In the FDTD grid it is necessary to discretize % the grading by averaging the material properties over a grid cell % width centered on each field component. As an example of the % implementation of this averaging, we take the integral of the % continuous sigma(x) in the PML region % % sigma_i = integral(sigma(x))/dx % % where the integral is over a single grid cell width in x,