[FDTD]FDTD（时域有限差分法）算法的Matlab源程序

离线alexis

UID ：1909

注册：2007-04-18
登录：2009-06-01
发帖：24
等级：仿真二级

0楼发表于: 2007-04-28 00:51:58

— 本帖被 casey 执行取消置顶操作(2008-09-03) —

[post]%***********************************************************************
%    3-D FDTD code with PEC boundaries
%***********************************************************************
%
%    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.
%
%    The computational domain is truncated using PEC boundary
%    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 "fdtd3D" at the MATLAB prompt.
%    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 "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

%***********************************************************************
%    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');
ylabel('k coordinate');

rect=get(gcf,'Position');
rect(1:2)=[0 0];

M=moviein(nmax/2,gcf,rect);

%***********************************************************************
%    BEGIN TIME-STEPPING LOOP
%***********************************************************************

for n=1:nmax

%***********************************************************************
%    Update electric fields
%***********************************************************************

ex(1:ie,2:je,2:ke)=ca*ex(1:ie,2:je,2:ke)+...
   cb*(hz(1:ie,2:je,2:ke)-hz(1:ie,1:je-1,2:ke)+...
   hy(1:ie,2:je,1:ke-1)-hy(1:ie,2:je,2:ke));

ey(2:ie,1:je,2:ke)=ca*ey(2:ie,1:je,2:ke)+...
   cb*(hx(2:ie,1:je,2:ke)-hx(2:ie,1:je,1:ke-1)+...
   hz(1:ie-1,1:je,2:ke)-hz(2:ie,1:je,2:ke));

ez(2:ie,2:je,1:ke)=ca*ez(2:ie,2:je,1:ke)+...
   cb*(hx(2:ie,1:je-1,1:ke)-hx(2:ie,2:je,1:ke)+...
   hy(2:ie,2:je,1:ke)-hy(1:ie-1,2:je,1:ke));

ez(is,js,1:ke)=ez(is,js,1:ke)+...
   srcconst*(n-ndelay)*exp(-((n-ndelay)^2/tau^2));

%***********************************************************************
%    Update magnetic fields
%***********************************************************************

hx(2:ie,1:je,1:ke)=hx(2:ie,1:je,1:ke)+...
   db*(ey(2:ie,1:je,2:kb)-ey(2:ie,1:je,1:ke)+...
   ez(2:ie,1:je,1:ke)-ez(2:ie,2:jb,1:ke));

hy(1:ie,2:je,1:ke)=hy(1:ie,2:je,1:ke)+...
   db*(ex(1:ie,2:je,1:ke)-ex(1:ie,2:je,2:kb)+...
   ez(2:ib,2:je,1:ke)-ez(1:ie,2:je,1:ke));

hz(1:ie,1:je,2:ke)=hz(1:ie,1:je,2:ke)+...
   db*(ex(1:ie,2:jb,2:ke)-ex(1:ie,1:je,2:ke)+...
   ey(1:ie,1:je,2:ke)-ey(2:ib,1:je,2:ke));

%***********************************************************************
%    Visualize fields
%************************** ..

未注册仅能浏览部分内容，查看全部内容及附件请先登录或注册

共2条评分