%2D FDTD TE mode with a plane wave source and a PML abc
%program edited by: K.B.A., NIP,U.P.,Diliman
%Date: January 2008
%this program borrowed the (2D FDTD)code by susan hagness
%***********************************************************************
%CHANGES: The metallic cylindrical scatterer is removed.
% The whole computational region becomes free space only.
% Program implements the TF/SF source formulation
%***********************************************************************
%Note:if you want to include the cylindrical metal, replace the sig
% expression with the commented sig in line 124
%
%***********************************************************************
% 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;
is1=7; %connecting boundaries
is2=ie-is1-1;
js1=7;
js2=je-js1-1;
%***********************************************************************
% Material parameters
%***********************************************************************
media=2;
eps=[1.0 1.0];
sig=[0.0 0.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
ey_inc=zeros(ib,1);
hz_inc=zeros(ie,1);
ey_low_m1=0;
ey_low_m2=0;
ey_high_m1=0;
ey_high_m2=0;
%***********************************************************************
% 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.
fdtd1.rar_fdtd
版权申诉
20 浏览量
2022-07-14
22:30:27
上传
评论
收藏 5KB RAR 举报
钱亚锋
- 粉丝: 87
- 资源: 1万+