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

%*********************************************************************** SmXoNiM"y  
%    3-D FDTD code with PEC boundaries h:Pfiw]  
%*********************************************************************** T3w%y`K  
% *>m[ZJd%=  
%    Program author: Susan C. Hagness +t!]nE#  
%                    Department of Electrical and Computer Engineering %ZVYgtk;*  
%                    University of Wisconsin-Madison =*UK!y?n  
%                    1415 Engineering Drive 9q;O`&  
%                    Madison, WI 53706-1691 [LVXXjkFI  
%                    608-265-5739 axQ>~vWN/  
%                    hagness@engr.wisc.edu ;lMvxt:  
% VG5+u,U6>  
%    Date of this version:  February 2000 m6Dm1'+  
% -P]O t>%S  
%    This MATLAB M-file implements the finite-difference time-domain S_lGrk\j  
%    solution of Maxwell's curl equations over a three-dimensional e!u]l  
%    Cartesian space lattice comprised of uniform cubic grid cells. i"Jy>'  
% P\"kr?jZP  
%    To illustrate the algorithm, an air-filled rectangular cavity SioeIXU  
%    resonator is modeled.  The length, width, and height of the h.<f%&)F  
%    cavity are 10.0 cm (x-direction), 4.8 cm (y-direction), and !HTOE@  
%    2.0 cm (z-direction), respectively. [YvS#M3T  
% `d <`>  
%    The computational domain is truncated using PEC boundary i?T-6
{3I  
%    conditions: NaR} 0  
%          ex(i,j,k)=0 on the j=1, j=jb, k=1, and k=kb planes J(1Tl  
%          ey(i,j,k)=0 on the i=1, i=ib, k=1, and k=kb planes J*@pM  
%          ez(i,j,k)=0 on the i=1, i=ib, j=1, and j=jb planes sI,cX#h&Y  
%    These PEC boundaries form the outer lossless walls of the cavity. ^t0!Dbx3SE  
% !LKxZ"  
%    The cavity is excited by an additive current source oriented 196aYLE  
%    along the z-direction.  The source waveform is a differentiated P1F-Wy1  
%    Gaussian pulse given by C[FHqo9M?H  
%          J(t)=-J0*(t-t0)*exp(-(t-t0)^2/tau^2), {"o9pIh{~  
%    where tau=50 ps.  The FWHM spectral bandwidth of this zero-dc- f&-`+V}U  
%    content pulse is approximately 7 GHz. The grid resolution }*x1e_m}H  
%    (dx = 2 mm) was chosen to provide at least 10 samples per ((M,6Q}  
%    wavelength up through 15 GHz. ^ :VH?I=  
% yP"2.9\erH  
%    To execute this M-file, type "fdtd3D" at the MATLAB prompt. p6JTNx
D  
%    This M-file displays the FDTD-computed Ez fields at every other 'FW?   
%    time step, and records those frames in a movie matrix, M, which &h:4TaD  
%    is played at the end of the simulation using the "movie" command. &xGdKH  
% /-M:6  
%*********************************************************************** [-
(^>Y  
OVL
VsNg  
clear 6,t6~Uo/  
WwmYJl0  
%*********************************************************************** ~WA@YjQ]  
%    Fundamental constants ay8]"sa  
%*********************************************************************** 7?dWAUF  
'sEnh<  
cc=2.99792458e8;            %speed of light in free space }Rh\JDiQ  
muz=4.0*pi*1.0e-7;          %permeability of free space \M`qaFan5^  
epsz=1.0/(cc*cc*muz);      %permittivity of free space b9H(w%7ucU  
C'#KTp4!1  
%*********************************************************************** rWht},-|1  
%    Grid parameters {B+}LL!  
%*********************************************************************** 7Ae`>5B#  
2ZK]}&yC  
ie=50;      %number of grid cells in x-direction _S9)<RVI+  
je=24;      %number of grid cells in y-direction "uH>S+%|b  
ke=10;      %number of grid cells in z-direction s&73g0$$  
l b(  
ib=ie+1; q1ysT.{p,  
jb=je+1; CCTU-Xz/  
kb=ke+1;  +n1!xv]  
'_v~+  
is=26;      %location of z-directed current source t<rIg1  
js=13;      %location of z-directed current source _H$Lu4b)N  
;Gr {  
kobs=5; fD!c t;UK  
VMah3T!  
dx=0.002;          %space increment of cubic lattice w=GMQ8  
dt=dx/(2.0*cc);    %time step f:-)S8OJ  
&d6@SQ  
nmax=500;          %total number of time steps !Fl'?Kz  
66val"^W  
O:?3B!wF  
%*********************************************************************** O84]J:b  
%    Differentiated Gaussian pulse excitation ;L-=z]IR,  
%*********************************************************************** j\C6k  
P@N+jS`Vf  
rtau=50.0e-12; 2@S}x@^  
tau=rtau/dt; TPp]UG  
ndelay=3*tau; "0eX/rY%  
srcconst=-dt*3.0e+11; JM7mQ'`Ud  
rXB;#ypO  
%*********************************************************************** TpU\IQ  
%    Material parameters }nrjA0WN  
%*********************************************************************** tpWGmjfo>  
ny-:%A  
eps=1.0; JxIJxhA>  
sig=0.0; |k.'w<6mb9  
Cm[}DB  
%*********************************************************************** /+*"*Br/  
%    Updating coefficients LS~at.3zX  
%*********************************************************************** ytttF5-  
b 3x|Dq.  
ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps)); 3d}v?q78  
cb=(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps)); 0B.Gt&Oal  
da=1.0; 2^r~->  
db=dt/muz/dx; 2yA)SGri  
v3r3$(Hr  
%*********************************************************************** =
MNp;  
%    Field arrays x2!R&q8U>  
%*********************************************************************** -6[DQB  
"2mPWRItO  
ex=zeros(ie,jb,kb); R4p Pt  
ey=zeros(ib,je,kb); ImN'o4vo  
ez=zeros(ib,jb,ke); <-I69`  
hx=zeros(ib,je,ke); r!R-3LO0s  
hy=zeros(ie,jb,ke); 0pE>O7  
hz=zeros(ie,je,kb); 4aO/^Hl  
`Gio 2gl9  
%*********************************************************************** !q\8`ss  
%    Movie initialization &ijz'Sg3  
%*********************************************************************** 5`^@k<  
O`Tz^Q/D  
tview(:,:)=ez(:,:,kobs); P&0eu  
sview(:,:)=ez(:,js,:); bYsX?0T!p  
@!`__>K  
subplot('position',[0.15 0.45 0.7 0.45]),pcolor(tview'); =|oi0  
shading flat; %]+R>+  
caxis([-1.0 1.0]); atyvo0fNd  
colorbar; $a_y-lY  
axis image; B=0^Rysg  
title(['Ez(i,j,k=5), time step = 0']); c}(H*VY2n  
xlabel('i coordinate'); 5/"$_7"{a  
ylabel('j coordinate'); Fi``l)Tt  
y[m,t}gi  
subplot('position',[0.15 0.10 0.7 0.25]),pcolor(sview'); <YvXyIs  
shading flat; I}PI  
caxis([-1.0 1.0]); *(9Tl]w  
colorbar; <r}wQ\F#  
axis image; cx[^D,usf~  
title(['Ez(i,j=13,k), time step = 0']); !|[rh,e]  
xlabel('i coordinate'); J*;RL`  
ylabel('k coordinate'); </[.1&S+\  
%,\JTN|g|A  
rect=get(gcf,'Position'); {lf{0c$X.  
rect(1:2)=[0 0]; o=Vs)8W  
9n
R\7!_  
M=moviein(nmax/2,gcf,rect); 2+9VDf2  
"D'e  
%*********************************************************************** 6?a`'&  
%    BEGIN TIME-STEPPING LOOP |nNcV~%~  
%*********************************************************************** A-!e$yz>  
PYCN3s#Gi  
for n=1:nmax /2p*uv}IP  
a{`hAI${  
%*********************************************************************** v~mVf.j1  
%    Update electric fields ~nA k-toJ  
%*********************************************************************** c_>Gl8J  
|N/
d}  
ex(1:ie,2:je,2:ke)=ca*ex(1:ie,2:je,2:ke)+... G4(R/<J,BQ  
                  cb*(hz(1:ie,2:je,2:ke)-hz(1:ie,1:je-1,2:ke)+... re*}a)iL  
                      hy(1:ie,2:je,1:ke-1)-hy(1:ie,2:je,2:ke)); }Ulxt:}   
 \0)jWCK  
ey(2:ie,1:je,2:ke)=ca*ey(2:ie,1:je,2:ke)+... -Q20af-  
                  cb*(hx(2:ie,1:je,2:ke)-hx(2:ie,1:je,1:ke-1)+... AR[m+E  
                      hz(1:ie-1,1:je,2:ke)-hz(2:ie,1:je,2:ke)); yCxYFi  
IR-n:z  
ez(2:ie,2:je,1:ke)=ca*ez(2:ie,2:je,1:ke)+... hU$aZ  
                  cb*(hx(2:ie,1:je-1,1:ke)-hx(2:ie,2:je,1:ke)+... $cSUB  
                      hy(2:ie,2:je,1:ke)-hy(1:ie-1,2:je,1:ke)); O,r;-t4vYU  
I[?\Or  
ez(is,js,1:ke)=ez(is,js,1:ke)+... nXT`7  
              srcconst*(n-ndelay)*exp(-((n-ndelay)^2/tau^2)); Mq\=pxC@  
^MpMqm1?8;  
%*********************************************************************** {;{U@Z  
%    Update magnetic fields t\,Y<9{w  
%*********************************************************************** n{gEIUo
#  
F
Ct<h/  
hx(2:ie,1:je,1:ke)=hx(2:ie,1:je,1:ke)+... @1*
^ttC  
                  db*(ey(2:ie,1:je,2:kb)-ey(2:ie,1:je,1:ke)+... ;FqmZjm  
                      ez(2:ie,1:je,1:ke)-ez(2:ie,2:jb,1:ke)); *D|a`R!Y  
:;_
khno  
hy(1:ie,2:je,1:ke)=hy(1:ie,2:je,1:ke)+... M`Q$-#E:  
                  db*(ex(1:ie,2:je,1:ke)-ex(1:ie,2:je,2:kb)+... # 55>?  
                      ez(2:ib,2:je,1:ke)-ez(1:ie,2:je,1:ke)); sT?{  
h ;*x1BVE  
hz(1:ie,1:je,2:ke)=hz(1:ie,1:je,2:ke)+... ei5YxV6I  
                  db*(ex(1:ie,2:jb,2:ke)-ex(1:ie,1:je,2:ke)+... }5+^  
                      ey(1:ie,1:je,2:ke)-ey(2:ib,1:je,2:ke)); H~FI@Cf$L  
=+iY<~8  
%*********************************************************************** WZ-~F/:c%  
%    Visualize fields @} Ig*@  
%*********************************************************************** }#^Cj;  
 {b!{~q  

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

hao tie zi ,wo ding

顶~~~~~~~~~~~~

不错，我顶一下！

很不错 7^@ 1cA=S  
定了

确实是好材料啊！！

haodongxi,xiexie

真的是不错的东东呢！

先研究一下先

基本的迭代式，没有ABC