[资料共享]Susan C. Hagness的1D/2D/3D FDTD（matlab）源码

学习下，初学者啊

wo xu yao de .duoxie! "vF7b|I  

好东东 (;2]`D [x  
好东东 +`+r\*C5  
您正在看的文章

谢谢lz分享

这个代码有了，谢谢楼主

谢谢谢谢谢谢谢谢谢谢谢谢

很好的资料，，看看啊 谢谢

%*********************************************************************** 2(-J9y|  
%     3-D FDTD code with PEC boundaries 5/v,|  
%*********************************************************************** y^rcUPLT  
% d},IQ,Az:Z  
%     Program author: Susan C. Hagness B?Vr9H7n  
%                     Department of Electrical and Computer Engineering @R5jUPUVV  
%                     University of Wisconsin-Madison B8jSdlvz  
%                     1415 Engineering Drive ){AtV&{$  
%                     Madison, WI 53706-1691 0{ZYYB&"~J  
%                     608-265-5739 6m.ChlO/  
%                     hagness@engr.wisc.edu 'IorjR@40  
% /Hox]
r]'e  
%     Date of this version:  February 2000 9`"o,wGX3  
% [79 eq=  
%     This MATLAB M-file implements the finite-difference time-domain jWn!96NhlL  
%     solution of Maxwell's curl equations over a three-dimensional .6,+q2tyk,  
%     Cartesian space lattice comprised of uniform cubic grid cells. *K9I+t"g  
%     LvWl*:z  
%     To illustrate the algorithm, an air-filled rectangular cavity Zz"b&`K  
%     resonator is modeled.  The length, width, and height of the bIFKP  
%     cavity are 10.0 cm (x-direction), 4.8 cm (y-direction), and 7x/S4Gs'4  
%     2.0 cm (z-direction), respectively. E+Jh4$x{  
% vv2N;/;I  
%     The computational domain is truncated using PEC boundary Qz T>h  
%     conditions: z (N3oBW  
%          ex(i,j,k)=0 on the j=1, j=jb, k=1, and k=kb planes 7eQE[C  
%          ey(i,j,k)=0 on the i=1, i=ib, k=1, and k=kb planes ^2gDhoO_  
%          ez(i,j,k)=0 on the i=1, i=ib, j=1, and j=jb planes vybQ}dscn  
%     These PEC boundaries form the outer lossless walls of the cavity. 1g_(xwUp+  
% !]%M  
%     The cavity is excited by an additive current source oriented &1O!guq%  
%     along the z-direction.  The source waveform is a differentiated ]=s!cfu  
%     Gaussian pulse given by AN!MFsk  
%          J(t)=-J0*(t-t0)*exp(-(t-t0)^2/tau^2), \2?p  
%     where tau=50 ps.  The FWHM spectral bandwidth of this zero-dc- L<kIzB !  
%     content pulse is approximately 7 GHz. The grid resolution U({20  
%     (dx = 2 mm) was chosen to provide at least 10 samples per s6#@S4^=\  
%     wavelength up through 15 GHz. \9U4V>p  
% ~R;9a"nr  
%     To execute this M-file, type "fdtd3D" at the MATLAB prompt. g/)$-Z)Nu  
%     This M-file displays the FDTD-computed Ez fields at every other +3HukoR(  
%     time step, and records those frames in a movie matrix, M, which :$bp4+
3>  
%     is played at the end of the simulation using the "movie" command. 5yvaY "B  
% L53qQej<  
%*********************************************************************** /jR]sC)xs  
;2o+|U@  
clear a,o>E4#c  
;l^4/BR  
%*********************************************************************** A)5-w`1  
%     Fundamental constants %+WIv+<  
%*********************************************************************** g$&uD  
45Z"U<I,9  
cc=2.99792458e8;            %speed of light in free space  R^%uEP  
muz=4.0*pi*1.0e-7;          %permeability of free space ]gN]Cw\L  
epsz=1.0/(cc*cc*muz);       %permittivity of free space V#j|_N1hm  
U=?hT&w\S  
%*********************************************************************** {K{&__Nk  
%     Grid parameters '%W'HqVcG1  
%*********************************************************************** 8q
)wT0A~  
h)%}O.ueB  
ie=50;       %number of grid cells in x-direction JvA6kw,  
je=24;       %number of grid cells in y-direction IrJPP2Q  
ke=10;       %number of grid cells in z-direction Mo+mO&B  
=)56]ki}  
ib=ie+1;     OL,3Jh% x  
jb=je+1;   <O`yM2/pS  
kb=ke+1;   ?
e? mg  
@ljvTgZ(X  
is=26;       %location of z-directed current source P;&rh U^[  
js=13;       %location of z-directed current source -Jo8jE~>V  
Km~\^(a '  
kobs=5; }3tbqFiH  
!J34yro+s  
dx=0.002;          %space increment of cubic lattice 2gn*B$a  
dt=dx/(2.0*cc);    %time step ? x1"uH  
xS~OAcxg  
nmax=500;          %total number of time steps "W_C%elg  
5K1WfdBX7)  
%*********************************************************************** bYr;~ ^  
%     Differentiated Gaussian pulse excitation L*ZC` .h  
%*********************************************************************** ]!>ThBMa  
dg%Orvuz  
rtau=50.0e-12; j@1rVOmK  
tau=rtau/dt; 6E9y[ %+  
ndelay=3*tau; KFCL|9P  
srcconst=-dt*3.0e+11; 8\-Q(9q(  
o<`)cb }  
%*********************************************************************** =AFTB<7-^  
%     Material parameters jL$&]sQ`O)  
%*********************************************************************** T[w]w  
*>Z|!{bI  
eps=1.0;
DRldRm/  
sig=0.0;         r]x;JBy  
O[p;IG
`  
%*********************************************************************** hcQvL>  
%     Updating coefficients KRS_6G],{  
%*********************************************************************** 4<S*gu*W  
)I80Nq  
ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps)); >*xa\ve  
cb=(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps)); ,&S^Ryc  
da=1.0; *7"R[!9  
db=dt/muz/dx; ,"U8Fgf[r  
^ <Z^3c>/  
%*********************************************************************** :Zy7h7P,lT  
%     Field arrays %21i#R`E  
%*********************************************************************** UcCkn7}  
R
P]hW{:U  
ex=zeros(ie,jb,kb); njk.$]M|nf  
ey=zeros(ib,je,kb); 32_{nLV$[  
ez=zeros(ib,jb,ke); MCamc  
hx=zeros(ib,je,ke); xl%!7?G|$>  
hy=zeros(ie,jb,ke); SnK j:|bV  
hz=zeros(ie,je,kb); Q AJX7  

qp(F}@  
%*********************************************************************** 1a#oJU  
%     Movie initialization r7t
N(2;5  
%*********************************************************************** yVThbL_YJ  
=d/\8\4  
tview(:,:)=ez(:,:,kobs); ETP}mo  
sview(:,:)=ez(:,js,:); 6qA48:/F=  
|>(@n{  
subplot('position',[0.15 0.45 0.7 0.45]),pcolor(tview'); Az.k6)~  
shading flat; im&E\`L7  
caxis([-1.0 1.0]); -b>"2B?  
colorbar; 8uyUvSB  
axis image; R5LzqT,/N:  
title(['Ez(i,j,k=5), time step = 0']); 0\tk/<w2  
xlabel('i coordinate'); 9wP_dJvb  
ylabel('j coordinate'); 7s%DM6li 6  
%K^l]tWa@  
subplot('position',[0.15 0.10 0.7 0.25]),pcolor(sview'); W11Wv&  
shading flat; ?^i$} .%W  
caxis([-1.0 1.0]); .p%p_  
colorbar; Q]_3 #_'  
axis image; ^'S0A=1  
title(['Ez(i,j=13,k), time step = 0']); ~Mv@Bl  
xlabel('i coordinate'); +w Oa  
ylabel('k coordinate'); '3Ir(]Wfd  
@Taj++ua  
rect=get(gcf,'Position'); T ,, Ao36  
rect(1:2)=[0 0]; L&1VPli  
pv2_A   
M=moviein(nmax/2,gcf,rect); _A*5BAB:h(  
o56_t{<  
%*********************************************************************** _S:6;_bz  
%     BEGIN TIME-STEPPING LOOP EG5'kYw2  
%*********************************************************************** `BKo`@  
_W41
;OY  
for n=1:nmax W G2 E3y  
   u,0N[.&N  
%*********************************************************************** `d#l o  
%     Update electric fields \Q"o\:IoIT  
%*********************************************************************** \H$Ps9Xh  
\14"Bgj1  
ex(1:ie,2:je,2:ke)=ca*ex(1:ie,2:je,2:ke)+... F_ ~L&jHP  
                   cb*(hz(1:ie,2:je,2:ke)-hz(1:ie,1:je-1,2:ke)+... k%\_UYa  
                       hy(1:ie,2:je,1:ke-1)-hy(1:ie,2:je,2:ke)); ?2VY^7N[  
r}
OK3J  
ey(2:ie,1:je,2:ke)=ca*ey(2:ie,1:je,2:ke)+... &x{CC@g/  
                   cb*(hx(2:ie,1:je,2:ke)-hx(2:ie,1:je,1:ke-1)+... !j8h$+:K  
                       hz(1:ie-1,1:je,2:ke)-hz(2:ie,1:je,2:ke)); @;}bBHQz{p  
                     @R
=gJ:&a  
ez(2:ie,2:je,1:ke)=ca*ez(2:ie,2:je,1:ke)+... .s31D%N  
                   cb*(hx(2:ie,1:je-1,1:ke)-hx(2:ie,2:je,1:ke)+... jsSxjf;O  
                       hy(2:ie,2:je,1:ke)-hy(1:ie-1,2:je,1:ke)); &
6A'}9Ch  
                     nAn/Vu  
ez(is,js,1:ke)=ez(is,js,1:ke)+... -|mWi  
               srcconst*(n-ndelay)*exp(-((n-ndelay)^2/tau^2)); ~LPxVYhK  
X>mY`$!/  
%*********************************************************************** [B9'/:  
%     Update magnetic fields &~P5[[Q  
%*********************************************************************** f[3DK
A  
>9c$2d|>  
hx(2:ie,1:je,1:ke)=hx(2:ie,1:je,1:ke)+... G l+[|?N  
                   db*(ey(2:ie,1:je,2:kb)-ey(2:ie,1:je,1:ke)+... bkkhx,Oi[G  
                       ez(2:ie,1:je,1:ke)-ez(2:ie,2:jb,1:ke)); @NGK2J  
                 PF
@+~FI  
hy(1:ie,2:je,1:ke)=hy(1:ie,2:je,1:ke)+... Y)Y`9u<?  
                   db*(ex(1:ie,2:je,1:ke)-ex(1:ie,2:je,2:kb)+... H4LZNko  
                       ez(2:ib,2:je,1:ke)-ez(1:ie,2:je,1:ke)); Efi@hdEV  
                 W=M`Bkw{  
hz(1:ie,1:je,2:ke)=hz(1:ie,1:je,2:ke)+... :mLcb.E  
                   db*(ex(1:ie,2:jb,2:ke)-ex(1:ie,1:je,2:ke)+...
oxE'u<  
                       ey(1:ie,1:je,2:ke)-ey(2:ib,1:je,2:ke)); ^yUel.N5"  
                     .lfKS!m2  
%*********************************************************************** HlqvXt\  
%     Visualize fields @cPb*  
%*********************************************************************** 2wE?O^J  
c[]_gUp8  
if mod(n,2)==0; ; >3q@9\D  
U'S}7g
ya  
timestep=int2str(n); X3#/|>
 
tview(:,:)=ez(:,:,kobs); FL!W oTB  
sview(:,:)=ez(:,js,:); ~<Z;)e  
TEl:;4  
subplot('position',[0.15 0.45 0.7 0.45]),pcolor(tview'); E4RvVfA0F  
shading flat; LRBcW;.Su  
caxis([-1.0 1.0]); jC1mui|Y^  
colorbar; pN/)$6=  
axis image; R6HMi#eF  
title(['Ez(i,j,k=5), time step = ',timestep]); LZm6\x  
xlabel('i coordinate'); ?Y2ZqI  
ylabel('j coordinate'); !=@Lyt)_b  
9Vz1*4Ln  
subplot('position',[0.15 0.10 0.7 0.25]),pcolor(sview'); vz:P2TkM  
shading flat; t4pc2b  
caxis([-1.0 1.0]); *@fR36  
colorbar; b:/;  
axis image; iph>"b$D
 
title(['Ez(i,j=13,k), time step = ',timestep]); Zh`lC1l'  
xlabel('i coordinate'); \|nF55W [  
ylabel('k coordinate');
R*IO%9O  
-frmvNJ F  
nn=n/2; Mk"+*G  
M(:,nn)=getframe(gcf,rect); 7PO3{I  
LziEF-_  
end; %)T>Wn%b]v  
V0 x[sEW  
%*********************************************************************** <qN0Q7  
%     END TIME-STEPPING LOOP \y^Od7F  
%*********************************************************************** pkpD1c^  
9K~2!<  
end $u_0"sUV  
wEENN_w  
movie(gcf,M,0,10,rect);

有没有一维的啊，我急需啊！

3D有没有PML边界的？ X{4jyi-<