FDTD(时域有限差分法)算法 下载本文

内容发布更新时间 : 2024/11/8 15:13:16星期一 下面是文章的全部内容请认真阅读。

% 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 three-dimensional % Cartesian space lattice comprised of uniform cubic grid cells. %

% To illustrate the algorithm, an air-filled rectangular cavity % resonator is modeled. The length, width, and height of the % cavity are 10.0 cm (x-direction), 4.8 cm (y-direction), and % 2.0 cm (z-direction), respectively. % conditions:

% ex(i,j,k)=0 on the j=1, j=jb, k=1, and k=kb planes % ey(i,j,k)=0 on the i=1, i=ib, k=1, and k=kb planes % ez(i,j,k)=0 on the i=1, i=ib, j=1, and j=jb planes

% These PEC boundaries form the outer lossless walls of the cavity. %

% The cavity is excited by an additive current source oriented % along the z-direction. The source waveform is a differentiated % Gaussian pulse given by

% J(t)=-J0*(t-t0)*exp(-(t-t0)^2/tau^2),

% where tau=50 ps. The FWHM spectral bandwidth of this zero-dc- % content pulse is approximately 7 GHz. The grid resolution % (dx = 2 mm) was chosen to provide at least 10 samples per % wavelength up through 15 GHz. %

% To execute this M-file, type \% This M-file displays the FDTD-computed Ez fields at every other % time step, and records those frames in a movie matrix, M, which % is played at the end of the simulation using the \%

%*********************************************************************** 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

%*********************************************************************** % Grid parameters

%***********************************************************************

ie=50; %number of grid cells in x-direction je=24; %number of grid cells in y-direction ke=10; %number of grid cells in z-direction ib=ie+1; jb=je+1; kb=ke+1;

is=26; %location of z-directed current source js=13; %location of z-directed current source kobs=5;

dx=0.002; %space increment of cubic lattice dt=dx/(2.0*cc); %time step

nmax=500; %total number of time steps

%*********************************************************************** % Differentiated Gaussian pulse excitation

%***********************************************************************

rtau=50.0e-12; tau=rtau/dt; ndelay=3*tau; srcconst=-dt*3.0e+11;

%*********************************************************************** % Material parameters

%*********************************************************************** eps=1.0; sig=0.0;

%*********************************************************************** % Updating coefficients

%***********************************************************************

ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps)); cb=(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps)); da=1.0; db=dt/muz/dx;

%*********************************************************************** % Field arrays

%***********************************************************************

ex=zeros(ie,jb,kb); ey=zeros(ib,je,kb); ez=zeros(ib,jb,ke); hx=zeros(ib,je,ke); hy=zeros(ie,jb,ke); hz=zeros(ie,je,kb);

%*********************************************************************** % Movie initialization

%***********************************************************************

tview(:,:)=ez(:,:,kobs); sview(:,:)=ez(:,js,:);

subplot('position',[0.15 0.45 0.7 0.45]),pcolor(tview'); shading flat; caxis([-1.0 1.0]); colorbar; axis image;

title(['Ez(i,j,k=5), time step = 0']); xlabel('i coordinate'); ylabel('j coordinate');

subplot('position',[0.15 0.10 0.7 0.25]),pcolor(sview'); shading flat; caxis([-1.0 1.0]); colorbar; axis image;

title(['Ez(i,j=13,k), time step = 0']); xlabel('i coordinate');