[求助]请大家帮我分析一下仿真三维自由空间的结果，让这UPML愁死了

图片:cuwu.jpg

h7J4 p  
利用3D-upml 仿真了自由空间，加的激励为高斯脉冲ez(is,js,ks)=ez(is,js,ks)+exp(-((n-45)^2/15^2)); Huf;A1.  
然后我就检测了一下ez(is,js,ks)随时间的变化图，结果是图中所示，不再是高斯脉冲形式 :ioD*k  
这样是不是不对啊。 程序师用的下面的那个，基本没改，就是改变了一个激励源。 mO&zE;/[  
%*********************************************************************** Ah_0o_Di  
%     3-D FDTD code with UPML absorbing boundary conditions ^ wb9n  
%*********************************************************************** dyVfDF  
% Vw+RRi(  
%     Program author: Keely J. Willis, Graduate Student pReSvF}}C  
%                     UW Computational Electromagnetics Laboratory Ti&v9re%wO  
%                           Director: Susan C. Hagness Y InPmR  
%                     Department of Electrical and Computer Engineering A 1B_EX.  
%                     University of Wisconsin-Madison b]
~  
%                     1415 Engineering Drive >DoP2]  
%                     Madison, WI 53706-1691 |[X-i["y  
%                     kjwillis@wisc.edu t(Gg 1  
% mOntc6&]  
%     Copyright 2005 S=qx,<J 39  
% 7.bPPr&  
%     This MATLAB M-file implements the finite-difference time-domain Mi047-% (  
%     solution of Maxwell's curl equations over a three-dimensional `gz/?q  
%     Cartesian space lattice comprised of uniform cubic grid cells. q+W*?a)  
%     ~'0W(~Q8  
%     The dimensions of the computational domain are 8.2 cm /A93mY[  
%     (x-direction), 3.4 cm (y-direction), and 3.2 cm (z-direction).  FQM9>l@6)>  
%     The grid is terminated with UPML absorbing boundary conditions. 2q ~y\fe  
% @6%o0p9zz  
%     An electric current source comprised of two collinear Jz components zj2l&)N  
%     (realizing a Hertzian dipole) excites a radially propagating wave.  Ir6g"kwCKq  
%     The current source is located in the center of the grid.  The A>%mJ3M  
%     source waveform is a differentiated Gaussian pulse given by 8y'.H21:;  
%          J(t)=J0*(t-t0)*exp(-(t-t0)^2/tau^2), 2<tU  
%     where tau=50 ps.  The FWHM spectral bandwidth of this zero-dc- Yz? 8n  
%     content pulse is approximately 7 GHz. The grid resolution Vr hd\  
%     (dx = 2 mm) was chosen to provide at least 10 samples per "=!sZO?3  
%     wavelength up through 15 GHz. XPT@ LM  
% La r9}nx0  
%     To execute this M-file, type "fdtd3D_UPML" at the MATLAB prompt.  ]}L tf,9  
% v)<|@TD)  
%     This code has been tested in the following Matlab environments: oa6&?4K?F  
%     Matlab version 6.1.0.450 Release 12.1 (May 18, 2001) 6-QcHJ>m6U  
%     Matlab version 6.5.1.199709 Release 13 Service Pack 1 (August 4, 2003) g[(Eh?]Sc  
%     Matlab version 7.0.0.19920 R14 (May 6, 2004) PG@6*E  
%     Matlab version 7.0.1.24704 R14 Service Pack 1 (September 13, 2004) Wd?=RO`a  
%     Matlab version 7.0.4.365 R14 Service Pack 2 (January 29, 2005) j.+}Z |  
% l
!j
,9wz7  
%     Note: if you are using Matlab version 6.x, you may wish to make e@{Rlz  
%     one or more of the following modifications to this code:
.~fov8  
%       --uncomment line numbers 485 and 486 OPC8fX5.  
%       --comment out line numbers 552 and 561 tgC)vZ&a  
% {9-n3j}  
%*********************************************************************** 3M5wF6nY[[  
clear 4DsHUc6  
%*********************************************************************** o_n 3.O=  
%     Fundamental constants BVu{To:g  
%*********************************************************************** #7=- zda5  
cc=2.99792458e8; U &
RZx&W  
muz=4.0*pi*1.0e-7; R})b%y`]  
epsz=1.0/(cc*cc*muz); gbziEjRe  
etaz=sqrt(muz/epsz); eDY)i9"W  
%*********************************************************************** 2#R8}\  
%     Material parameters R
<;;Ph  
%*********************************************************************** fT.MglJcb  
mur=1.0; G(gZL%M6  
epsr=1.0; vfdTGM`3  
eta=etaz*sqrt(mur/epsr); 23/;W|   
%*********************************************************************** Z+2 j(  
%     Grid parameters 75eZhs[b  
% ?<1~KLPMhY  
%     Each grid size variable name describes the number of sampled points Kon|TeC>d  
%     for a particular field component in the direction of that component. @0-vf>e3-  
%     Additionally, the variable names indicate the region of the grid - *v)sP"@  
%     for which the dimension is relevant.  For example, ie_tot is the
,M=s3D8C  
%     number of sample points of Ex along the x-axis in the total \{;3'<  
%     computational grid, and jh_bc is the number of sample points of Hy sf8F h  
%     along the y-axis in the y-normal UPML regions. G{gc]7\=Cd  
% @@H?w7y?&  
%*********************************************************************** hkx(r5o  
ie=60;          % Size of main grid Lmw
4  
je=60; i0rh{Ko  
ke=60; Y=Om0=v  
ih=ie+1; 7'Gkip  
jh=je+1;   5=/j  
kh=ke+1;    bU$M)  
upml=10;        % Thickness of PML boundaries >1hhz  
ih_bc=upml+1; O3tw@ &k
 
jh_bc=upml+1; Kiq
[PK  
kh_bc=upml+1; fPq)Lx1'  
ie_tot=ie+2*upml;          % Size of total computational domain zoZ10?ojC  
je_tot=je+2*upml;        ~nb1c:F  
ke_tot=ke+2*upml;        pyp0SGCM:  
ih_tot=ie_tot+1; ]nq/yAF%  
jh_tot=je_tot+1;          v(=E R%  
kh_tot=ke_tot+1;          xc,Wm/[  
is=round(ih_tot/2);         % Location of z-directed current source @4=Az1W*  
js=round(jh_tot/2); _ O;R  
ks=round(ke_tot/2); \`R8s_S  
%*********************************************************************** dP=,<H#]m  
%     Fundamental grid parameters #~[{*[B+  
%*********************************************************************** VM<$!Aaz  
delta=0.002;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3,1HD_  
dt=delta*sqrt(epsr*mur)/(2.0*cc); d fSj= 4  
nmax=100; t"P:}ps{?  
%*********************************************************************** s&1
}^'|  
%     Differentiated Gaussian pulse excitation iH8V]%  
%*********************************************************************** j7~Rw"(XQc  
rtau=50.0e-12; a(lmm@;V<  
tau=rtau/dt; BH0s` K"  
ndelay=3*tau; 8=OpX,t(  
J0=1.0*epsz; C4)m4r%  
%*********************************************************************** \i3)/sZ?l  
%     Initialize field arrays P DwBSj  
%*********************************************************************** !Y10UmMu  
ex=zeros(ie_tot,jh_tot,kh_tot); '<xV]k|v  
ey=zeros(ih_tot,je_tot,kh_tot); #wp~lW9!s9  
ez=zeros(ih_tot,jh_tot,ke_tot); yiMqe^zy  
dx=zeros(ie_tot,jh_tot,kh_tot); ^w_\D?  
dy=zeros(ih_tot,je_tot,kh_tot); C(kL
=WD   
dz=zeros(ih_tot,jh_tot,ke_tot); +mC?.B2D  
hx=zeros(ih_tot,je_tot,ke_tot);  rp=Y }  
hy=zeros(ie_tot,jh_tot,ke_tot); ?{\h`+A  
hz=zeros(ie_tot,je_tot,kh_tot); f/t`B^}@  
bx=zeros(ih_tot,je_tot,ke_tot); g0#w 4rGF)  
by=zeros(ie_tot,jh_tot,ke_tot); ^.c<b_(=h  
bz=zeros(ie_tot,je_tot,kh_tot); -Eoq#ULvR  
%*********************************************************************** ATjE8!gO!  
%     Initialize update coefficient arrays ydMSL25<+  
%*********************************************************************** N1hj[G[H"  
C1ex=zeros(size(ex)); .0p'G}1  
C2ex=zeros(size(ex)); tgY/8&$M  
C3ex=zeros(size(ex));  vURgR  
C4ex=zeros(size(ex)); [DvQk?,t  
C5ex=zeros(size(ex)); i5SDy(?r  
C6ex=zeros(size(ex)); u{S"NEc  
C1ey=zeros(size(ey)); kAq#cLprG  
C2ey=zeros(size(ey)); \LM.>vJ  
C3ey=zeros(size(ey)); r?/!VO-*N  
C4ey=zeros(size(ey)); }^2'@y!(  
C5ey=zeros(size(ey)); ,n$NF0^l  
C6ey=zeros(size(ey)); N
v6=[_D  
C1ez=zeros(size(ez)); CW -[c  
C2ez=zeros(size(ez)); ~l]g4iEp  
C3ez=zeros(size(ez)); M2I*_pI  
C4ez=zeros(size(ez)); "(Nt9K%P)  
C5ez=zeros(size(ez)); b-RuUfUn0  
C6ez=zeros(size(ez)); d<[L^s9  
D1hx=zeros(size(hx)); vbfQy2q  
D2hx=zeros(size(hx)); \-
8aTF  
D3hx=zeros(size(hx)); k:URP`w[X=  
D4hx=zeros(size(hx)); t\-|J SZ  
D5hx=zeros(size(hx)); ._q<~_~R  
D6hx=zeros(size(hx)); a|z@5r%  
D1hy=zeros(size(hy)); c$x>6&&L  
D2hy=zeros(size(hy)); J3$@: S'  
D3hy=zeros(size(hy)); U:bnX51D4  
D4hy=zeros(size(hy)); "O~kIT?/v  
D5hy=zeros(size(hy)); xtut S  
D6hy=zeros(size(hy)); 49YN@PXC  
D1hz=zeros(size(hz));  F |aLF{  
D2hz=zeros(size(hz)); klnk{R.>|  
D3hz=zeros(size(hz)); <GLn!~Px@5  
D4hz=zeros(size(hz)); -_= m j  
D5hz=zeros(size(hz)); gt\kTn."  
D6hz=zeros(size(hz)); Ey%KbvNv  
%*********************************************************************** vsJDVJ +=  
%     Update coefficients, as described in Section 7.8.2. e-#Vs{?|r  
% w-3Lw<  
%     In order to simplify the update equations used in the time-stepping 'V\V=yc1  
%     loop, we implement UPML update equations throughout the entire e xkPu-[W  
%     grid.  In the main grid, the electric-field update coefficients of a%5/Oc[[  
%     Equations 7.91a-f and the correponding magnetic field update vRH2[{KQ9  
%     coefficients extracted from Equations 7.89 and 7.90 are simplified 1u"#rC>7.4  
%     for the main grid (free space) and calculated below. *,28@_EwY  
% EI496bsRHm  
%*********************************************************************** `pF7B6[B  
C1=1.0; ] !n3j=*  
C2=dt; Nh\o39=  
C3=1.0; $laUkD#vz  
C4=1.0/2.0/epsr/epsr/epsz/epsz; C7K]c4T  
C5=2.0*epsr*epsz; =MT'e,T  
C6=2.0*epsr*epsz; qG +PqK;  
D1=1.0; M0$E_*  
D2=dt; ^I)+u>fJ  
D3=1.0; -
d^'-s  
D4=1.0/2.0/epsr/epsz/mur/muz; U+zntB  
D5=2.0*epsr*epsz; 5v^tPGg4  
D6=2.0*epsr*epsz; tG~[E,/`  
%*********************************************************************** =_CH$F!U  
%     Initialize main grid update coefficients ;SX~u*`R  
%*********************************************************************** w~yC^`  
C1ex(:,jh_bc:jh_tot-upml,:)=C1;     :6]qr86  
C2ex(:,jh_bc:jh_tot-upml,:)=C2; [&kz4_  
C3ex(:,:,kh_bc:kh_tot-upml)=C3; KDb`g}1Q  
C4ex(:,:,kh_bc:kh_tot-upml)=C4; *K BaKS  
C5ex(ih_bc:ie_tot-upml,:,:)=C5; !t}yoN n|  
C6ex(ih_bc:ie_tot-upml,:,:)=C6; _
fSBb<  
C1ey(:,:,kh_bc:kh_tot-upml)=C1; ]^,!;do  
C2ey(:,:,kh_bc:kh_tot-upml)=C2; Ss_}@p ^  
C3ey(ih_bc:ih_tot-upml,:,:)=C3; | ^G38  
C4ey(ih_bc:ih_tot-upml,:,:)=C4; O{0it6  
C5ey(:,jh_bc:je_tot-upml,:)=C5; $9@AwS@Uu  
C6ey(:,jh_bc:je_tot-upml,:)=C6; 'V&Tlw|  
C1ez(ih_bc:ih_tot-upml,:,:)=C1; 1 J}ML}h)  
C2ez(ih_bc:ih_tot-upml,:,:)=C2; ?!O4ia3nFk  
C3ez(:,jh_bc:jh_tot-upml,:)=C3; <W+9h0c  
C4ez(:,jh_bc:jh_tot-upml,:)=C4; Jt0U`
_  
C5ez(:,:,kh_bc:ke_tot-upml)=C5; F x^X(!)~]  
C6ez(:,:,kh_bc:ke_tot-upml)=C6; F@mxd  
D1hx(:,jh_bc:je_tot-upml,:)=D1; Tgh?=]H  
D2hx(:,jh_bc:je_tot-upml,:)=D2; Hg$7[um  
D3hx(:,:,kh_bc:ke_tot-upml)=D3; a ,"   
D4hx(:,:,kh_bc:ke_tot-upml)=D4; J 5xMA-  
D5hx(ih_bc:ih_tot-upml,:,:)=D5; VfUHqdg-  
D6hx(ih_bc:ih_tot-upml,:,:)=D6; BWNI|pq)v  
D1hy(:,:,kh_bc:ke_tot-upml)=D1; RC?vU  
D2hy(:,:,kh_bc:ke_tot-upml)=D2; >P]gjYN  
D3hy(ih_bc:ie_tot-upml,:,:)=D3; pw`'q(ad  
D4hy(ih_bc:ie_tot-upml,:,:)=D4; )j\_*SoH  
D5hy(:,jh_bc:jh_tot-upml,:)=D5; *F!1xyg  
D6hy(:,jh_bc:jh_tot-upml,:)=D6; Vf<q-3q  
D1hz(ih_bc:ie_tot-upml,:,:)=D1; %>5>wP   
D2hz(ih_bc:ie_tot-upml,:,:)=D2; =-,'LOE  
D3hz(:,jh_bc:je_tot-upml,:)=D3; p$uPj*  
D4hz(:,jh_bc:je_tot-upml,:)=D4; 7O:g;UI#  
D5hz(:,:,kh_bc:kh_tot-upml)=D5; (][-()YV  
D6hz(:,:,kh_bc:kh_tot-upml)=D6; ~@(
C+3,  
%*********************************************************************** XUnw*3tPJ  
%     Fill in PML regions 4c[/%e:\-  
% N9*:]a  
%     PML theory describes a continuous grading of the material properties v\_\bT1  
%     over the PML region.  In the FDTD grid it is necessary to discretize =3`|D0E  
%     the grading by averaging the material properties over a grid cell E^Q J50  
%     width centered on each field component.  As an example of the 9M5W4&  
%     implementation of this averaging, we take the integral of the |* ^LsuFb  
%     continuous sigma(x) in the PML region #sxv?r  
%   qDU4W7|T`  
%         sigma_i = integral(sigma(x))/delta vn!3Z!dm(  
%   E(vO^)#  
%     where the integral is over a single grid cell width in x, and is 0a bQY  
%     bounded by x1 and x2.  Applying this to the polynomial grading of k1ja ([Q  
%     Equation 7.60a produces K/j u=>  
% tY: Nq*@  
%         sigma_i = (x2^(m+1)-x1^(m+1))*sigmam/(delta*(m+1)*d^m) crN*eFeW  
% 57=d;Yg e  
%     where sigmam is the maximum value of sigma as described by Equation noM=8C&U  
%     7.62. "eqzn KT%u  
%         lQBEq"7$  
%     The definitions of x1 and x2 depend on the position of the field LNyrIk/1  
%     component within the grid cell.  We have either ]^T-X/v9  
% %V+"i_{m  
%         x1 = (i-0.5)*delta,  x2 = (i+0.5)*delta v1Q78P  
%  QZIzddwp  
%     or >239SyC-,  
%  Sc/$2gSG  
%         x1 = (i)*delta,      x2 = (i+1)*delta }2iR=$2  
% +vxOCN4}v  
%     where i varies over the PML region. W6_rSVm  
%  `( w"{8laB  
%*********************************************************************** 2pU'&8  
rmax=exp(-16);  %desired reflection error, designated as R(0) in Equation 7.62 ?5/7 @V  
orderbc=4;      %order of the polynomial grading, designated as m in Equation 7.60a,b iJZNSRQJ}r  
%   x-varying material properties ,aa 4Kh  
delbc=upml*delta; g7a446QR\K  
sigmam=-log(rmax)*(orderbc+1.0)/(2.0*eta*delbc); cpALs1j:  
sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1.0)); O6v
xp?:^  
kmax=1; ^9]iUx  
kfactor=(kmax-1.0)/delta/(orderbc+1.0)/delbc^orderbc; 3W]gn8  
for i=1:upml U|VL+9#hd  
    .g1x$cQ1<  
    % Coefficients for field components in the center of the grid cell j--byk6PB  
    x1=(upml-i+1)*delta; VBw5[  
    x2=(upml-i)*delta; 'nB
P%  
    sigma=sigfactor*(x1^(orderbc+1)-x2^(orderbc+1)); 0{vH.b @  
    ki=1+kfactor*(x1^(orderbc+1)-x2^(orderbc+1)); XH@(V4J(.  
    facm=(2*epsr*epsz*ki-sigma*dt); j/fniyJ)  
    facp=(2*epsr*epsz*ki+sigma*dt); B\e*-:pq>  
    C5ex(i,:,:)=facp; x
)GheM^  
    C5ex(ie_tot-i+1,:,:)=facp; Pq8oK'z-  
    C6ex(i,:,:)=facm; g6;smtu_T  
    C6ex(ie_tot-i+1,:,:)=facm; p$qk\efv*4  
    D1hz(i,:,:)=facm/facp; gPb.%^p  
    D1hz(ie_tot-i+1,:,:)=facm/facp; YB<nz<;JR  
    D2hz(i,:,:)=2.0*epsr*epsz*dt/facp; k!}(a0h  
    D2hz(ie_tot-i+1,:,:)=2.0*epsr*epsz*dt/facp; {L~dER  
    D3hy(i,:,:)=facm/facp; qw?(^uZNW  
    D3hy(ie_tot-i+1,:,:)=facm/facp; g(;OUkj$Zp  
    D4hy(i,:,:)=1.0/facp/mur/muz; yj#*H  
    D4hy(ie_tot-i+1,:,:)=1.0/facp/mur/muz; w;j<$<4=7  
    % Coefficients for field components on the grid cell boundary gPT_}#_GxM  
    x1=(upml-i+1.5)*delta; =B5{7g\  
    x2=(upml-i+0.5)*delta; MIn_?r  
    sigma=sigfactor*(x1^(orderbc+1)-x2^(orderbc+1)); 4dcm)Xr  
    ki=1.0+kfactor*(x1^(orderbc+1)-x2^(orderbc+1)); @!u{>!~0  
    facm=(2.0*epsr*epsz*ki-sigma*dt); 6@t&  
    facp=(2.0*epsr*epsz*ki+sigma*dt); p)ig~kk`  
    C1ez(i,:,:)=facm/facp; *YL86R+U  
    C1ez(ih_tot-i+1,:,:)=facm/facp; /t`\b [  
    C2ez(i,:,:)=2.0*epsr*epsz*dt/facp; yNowhh
 
    C2ez(ih_tot-i+1,:,:)=2.0*epsr*epsz*dt/facp; 8Znr1=1   
    C3ey(i,:,:)=facm/facp; J-<_e??  
    C3ey(ih_tot-i+1,:,:)=facm/facp; 0c`nk\vUy  
    C4ey(i,:,:)=1.0/facp/epsr/epsz; Z\x
nPhV  
    C4ey(ih_tot-i+1,:,:)=1.0/facp/epsr/epsz; GK!@|Kk8q7  
    D5hx(i,:,:)=facp; Bv!{V)$  
    D5hx(ih_tot-i+1,:,:)=facp; T!ZjgCY}  
    D6hx(i,:,:)=facm; Dmr*Lh~  
    D6hx(ih_tot-i+1,:,:)=facm; 8V 4e\q  
    p:5NMo  
end 2l+L96  
%   PEC walls
AvB=/p@]  
C1ez(1,:,:)=-1.0; iC
/*d  
C1ez(ih_tot,:,:)=-1.0; cv-rEHT  
C2ez(1,:,:)=0.0; lwrh4<~\,*  
C2ez(ih_tot,:,:)=0.0; 2#sFY/@  
C3ey(1,:,:)=-1.0; ..X_nF  
C3ey(ih_tot,:,:)=-1.0; At?|[%<`  
C4ey(1,:,:)=0.0; iw3\`,5   
C4ey(ih_tot,:,:)=0.0; X:$vP'B>  
%   y-varying material properties NsP=l]   
delbc=upml*delta; CuvY^["
 
sigmam=-log(rmax)*epsr*epsz*cc*(orderbc+1.0)/(2.0*delbc); PfRA\  
sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1.0)); P.C?/7$7Z+  
kmax=1.0; *DC/O( 0  
kfactor=(kmax-1.0)/delta/(orderbc+1.0)/delbc^orderbc; q,O_y<uw  
for j=1:upml A|
S)cr8z  
    mL
2
J  
    % Coefficients for field components in the center of the grid cell _uLpU4# ?  
    y1=(upml-j+1)*delta; _:=\h5}8  
    y2=(upml-j)*delta; ?]$<Ufr  
    sigma=sigfactor*(y1^(orderbc+1)-y2^(orderbc+1)); s_XCKhN:  
    ki=1+kfactor*(y1^(orderbc+1)-y2^(orderbc+1)); &1yJrj9y  
    facm=(2*epsr*epsz*ki-sigma*dt); 0NGth(2  
    facp=(2*epsr*epsz*ki+sigma*dt); ( _6j@?u  
    D 5n\h5  
    C5ey(:,j,:)=facp; iD G&Muc  
    C5ey(:,je_tot-j+1,:)=facp; Z4Dx:m
-  
    C6ey(:,j,:)=facm; tWTHyL  
    C6ey(:,je_tot-j+1,:)=facm; 5WJ ~%"O  
    D1hx(:,j,:)=facm/facp; boGdZ2$h4  
    D1hx(:,je_tot-j+1,:)=facm/facp; &:&l+  
    D2hx(:,j,:)=2*epsr*epsz*dt/facp; ]vvA]e  
    D2hx(:,je_tot-j+1,:)=2*epsr*epsz*dt/facp; sTl^j gV7j  
    D3hz(:,j,:)=facm/facp; 7@\.()  
    D3hz(:,je_tot-j+1,:)=facm/facp; /MY's&D(  
    D4hz(:,j,:)=1/facp/mur/muz; [`b,SX x  
    D4hz(:,je_tot-j+1,:)=1/facp/mur/muz; BVNJas  
    Q=Mv"~2>B  
    % Coefficients for field components on the grid cell boundary dGm%If9P  
    y1=(upml-j+1.5)*delta; zGR,}v%%  
    y2=(upml-j+0.5)*delta; vJ}WNvncVF  
    sigma=sigfactor*(y1^(orderbc+1)-y2^(orderbc+1)); {)lZfj}l  
    ki=1+kfactor*(y1^(orderbc+1)-y2^(orderbc+1)); L6PgWc;m  
    facm=(2*epsr*epsz*ki-sigma*dt); &q
G/\  
    facp=(2*epsr*epsz*ki+sigma*dt);    W&)f#/M8  
     yuBRYy#E|%  
    C1ex(:,j,:)=facm/facp; <O0tg[ub  
    C1ex(:,jh_tot-j+1,:)=facm/facp; ;_c&J&I  
    C2ex(:,j,:)=2*epsr*epsz*dt/facp; BqX"La,  
    C2ex(:,jh_tot-j+1,:)=2*epsr*epsz*dt/facp; -Pp{aFe  
    C3ez(:,j,:)=facm/facp; Mf#@8"l  
    C3ez(:,jh_tot-j+1,:)=facm/facp; Zf3(! a[  
    C4ez(:,j,:)=1/facp/epsr/epsz; r9WR1&T)  
    C4ez(:,jh_tot-j+1,:)=1/facp/epsr/epsz;   hmo4H3g!N  
    D5hy(:,j,:)=facp; N,-C+r5}<4  
    D5hy(:,jh_tot-j+1,:)=facp; AXHY
$f|  
    D6hy(:,j,:)=facm; Os'E7;:1h  
    D6hy(:,jh_tot-j+1,:)=facm; k
C%H E  
end py:L-5  
%   PEC walls 5>.ATfAsV  
C1ex(:,1,:)=-1; RU7+$Z0K  
C1ex(:,jh_tot,:)=-1; Si]?4:E7=  
C2ex(:,1,:)=0; XYuX+&XW/  
C2ex(:,jh_tot,:)=0; Lw,}wM5X  
C3ez(:,1,:)=-1; ) hs&?:)  
C3ez(:,jh_tot,:)=-1; &uRT/+18W3  
C4ez(:,1,:)=0; v#&;z_I+  
C4ez(:,jh_tot,:)=0;   O}zHkcL  
%   z-varying material properties +;lDU}$  
delbc=upml*delta; P
|@[D=y  
sigmam=-log(rmax)*epsr*epsz*cc*(orderbc+1)/(2*delbc); 5?SE?VC=t  
sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1)); mMwV5\(  
kmax=1; '1LN)Yw  
kfactor=(kmax-1)/delta/(orderbc+1)/delbc^orderbc; )/@KdEA:  
for k=1:upml 4"kc(J`c  
    % Coefficients for field components in the center of the grid cell #%k_V+o3  
    z1=(upml-k+1)*delta; H`)eT6:|/  
    z2=(upml-k)*delta; Z9$pY=8^?  
    sigma=sigfactor*(z1^(orderbc+1)-z2^(orderbc+1));  94PI  
    ki=1+kfactor*(z1^(orderbc+1)-z2^(orderbc+1)); `WOoC  
    facm=(2*epsr*epsz*ki-sigma*dt); (0q`eO2  
    facp=(2*epsr*epsz*ki+sigma*dt); y EfAa6  
    9e K~g0m  
    C5ez(:,:,k)=facp; :XFQ}Cl  
    C5ez(:,:,ke_tot-k+1)=facp; 5@%Gq)z5  
    C6ez(:,:,k)=facm; 3WF]%P%  
    C6ez(:,:,ke_tot-k+1)=facm; ~6kF`}5  
    D1hy(:,:,k)=facm/facp; dY&v(~&;]  
    D1hy(:,:,ke_tot-k+1)=facm/facp; T$"~Vu  
    D2hy(:,:,k)=2*epsr*epsz*dt/facp;
X}(X\rp  
    D2hy(:,:,ke_tot-k+1)=2*epsr*epsz*dt/facp; @[9  
    D3hx(:,:,k)=facm/facp; Nuot[1kS  
    D3hx(:,:,ke_tot-k+1)=facm/facp; YGOkq
I  
    D4hx(:,:,k)=1/facp/mur/muz; yZ,pH1  
    D4hx(:,:,ke_tot-k+1)=1/facp/mur/muz; oQ r.cKD ?  
    M?sax+'  
    % Coefficients for field components on the grid cell boundary / ^d9At614  
    z1=(upml-k+1.5)*delta; kE'p=dXx  
    z2=(upml-k+0.5)*delta; G<-KwGy,D  
    sigma=sigfactor*(z1^(orderbc+1)-z2^(orderbc+1)); Yc]k<tQ  
    ki=1+kfactor*(z1^(orderbc+1)-z2^(orderbc+1)); 8%m\J:eR  
    facm=(2*epsr*epsz*ki-sigma*dt); 0SYJ*7lPX  
    facp=(2*epsr*epsz*ki+sigma*dt); x}K|\KXy  
    a5/, O4Q  
    C1ey(:,:,k)=facm/facp; S1^/W-yoc~  
    C1ey(:,:,kh_tot-k+1)=facm/facp; wi7Br&bGi  
    C2ey(:,:,k)=2*epsr*epsz*dt/facp; -c?x5
/@3  
    C2ey(:,:,kh_tot-k+1)=2*epsr*epsz*dt/facp; w/o^OjwQ  
    C3ex(:,:,k)=facm/facp; *W\3cS  
    C3ex(:,:,kh_tot-k+1)=facm/facp; [k\VUg:P  
    C4ex(:,:,k)=1/facp/epsr/epsz; {AD-p!6G  
    C4ex(:,:,kh_tot-k+1)=1/facp/epsr/epsz; z3i`O La  
    D5hz(:,:,k)=facp; ,\"x#Cc f  
    D5hz(:,:,kh_tot-k+1)=facp; 'uL$j
=vB  
    D6hz(:,:,k)=facm; ?%~p@  
    D6hz(:,:,kh_tot-k+1)=facm; -@>]iBl  
end 51|s2+GG  
%   PEC walls #oTVfY#  
C1ey(:,:,1)=-1; gPB=Z!  
C1ey(:,:,kh_tot)=-1; }EwE#sZ#  
C2ey(:,:,1)=0; \-]Jm[]^  
C2ey(:,:,kh_tot)=0; ##Pzc~xSn  
C3ex(:,:,1)=-1; 1
$03:ve1  
C3ex(:,:,kh_tot)=-1; 13/,^?  
C4ex(:,:,1)=0; YQC.jnb2  
C4ex(:,:,kh_tot)=0; eR5q3E/;G  
%figure &%j`WF4p  
%set(gcf,'DoubleBuffer','on') TxK v!-1  
%*********************************************************************** {l$)X  
%     Begin time stepping loop iu=Mq|t0  
%***********************************************************************
6jc5B#  
d0=zeros(1,nmax);d1=zeros(1,nmax);d2=zeros(1,nmax); k8}fKVU;  
for n=1:nmax yJC: bD1xi  
    h SS9mQ  
    % Update magnetic field 9D\E0YG X/  
    bstore=bx; cQ-#]  
    bx(2:ie_tot,:,:)=D1hx(2:ie_tot,:,:).*  bx(2:ie_tot,:,:)-... @JhkUGG]p  
                     D2hx(2:ie_tot,:,:).*((ez(2:ie_tot,2:jh_tot,:)-ez(2:ie_tot,1:je_tot,:))-... ]I,&Bme  
                                          (ey(2:ie_tot,:,2:kh_tot)-ey(2:ie_tot,:,1:ke_tot)))./delta; 4h!yh2c..  
    hx(2:ie_tot,:,:)= D3hx(2:ie_tot,:,:).*hx(2:ie_tot,:,:)+... seK;TQ3/7  
                      D4hx(2:ie_tot,:,:).*(D5hx(2:ie_tot,:,:).*bx(2:ie_tot,:,:)-... ;~:Z~8+{c  
                                           D6hx(2:ie_tot,:,:).*bstore(2:ie_tot,:,:)); \Qah*1  
    bstore=by; vzI>:Bf  
    by(:,2:je_tot,:)=D1hy(:,2:je_tot,:).*  by(:,2:je_tot,:)-... D.{vuftu  
                     D2hy(:,2:je_tot,:).*((ex(:,2:je_tot,2:kh_tot)-ex(:,2:je_tot,1:ke_tot))-... ) 6QJZ$  
                                          (ez(2:ih_tot,2:je_tot,:)-ez(1:ie_tot,2:je_tot,:)))./delta; 
;hq_}.  
    hy(:,2:je_tot,:)= D3hy(:,2:je_tot,:).*hy(:,2:je_tot,:)+... &2 Yo  
                      D4hy(:,2:je_tot,:).*(D5hy(:,2:je_tot,:).*by(:,2:je_tot,:)-... |+Rx)  
                                           D6hy(:,2:je_tot,:).*bstore(:,2:je_tot,:)); 3lWGa7<4Z  
    bstore=bz; 5,MM`:{{  
    bz(:,:,2:ke_tot)=D1hz(:,:,2:ke_tot).*  bz(:,:,2:ke_tot)-... npkT>dB+  
                     D2hz(:,:,2:ke_tot).*((ey(2:ih_tot,:,2:ke_tot)-ey(1:ie_tot,:,2:ke_tot))-... O
XM=@B<"  
                                          (ex(:,2:jh_tot,2:ke_tot)-ex(:,1:je_tot,2:ke_tot)))./delta; nw/g[/<;  
    hz(:,:,2:ke_tot)= D3hz(:,:,2:ke_tot).*hz(:,:,2:ke_tot)+... paV1o>_Rd  
                      D4hz(:,:,2:ke_tot).*(D5hz(:,:,2:ke_tot).*bz(:,:,2:ke_tot)-... VO:  
                                           D6hz(:,:,2:ke_tot).*bstore(:,:,2:ke_tot)); )fXw~  
    %1k"K~eu  
    % Update electric field <`SA>P  
    dstore=dx; i`l;k~rP  
    dx(:,2:je_tot,2:ke_tot)=C1ex(:,2:je_tot,2:ke_tot).*  dx(:,2:je_tot,2:ke_tot)+... h!(# /  
                            C2ex(:,2:je_tot,2:ke_tot).*((hz(:,2:je_tot,2:ke_tot)-hz(:,1:je_tot-1,2:ke_tot))-... 4uO88[=  
                                                        (hy(:,2:je_tot,2:ke_tot)-hy(:,2:je_tot,1:ke_tot-1)))./delta; F-:AT$Ok  
    ex(:,2:je_tot,2:ke_tot)=C3ex(:,2:je_tot,2:ke_tot).*ex(:,2:je_tot,2:ke_tot)+... &&<l}E  
                            C4ex(:,2:je_tot,2:ke_tot).*(C5ex(:,2:je_tot,2:ke_tot).*dx(:,2:je_tot,2:ke_tot)-... <(Ar[Rp  
                                                        C6ex(:,2:je_tot,2:ke_tot).*dstore(:,2:je_tot,2:ke_tot)); G2wSd'n*y  
    dstore=dy; W#g!Usf:/  
    dy(2:ie_tot,:,2:ke_tot)=C1ey(2:ie_tot,:,2:ke_tot).*  dy(2:ie_tot,:,2:ke_tot)+... X1B)(|7$  
                            C2ey(2:ie_tot,:,2:ke_tot).*((hx(2:ie_tot,:,2:ke_tot)-hx(2:ie_tot,:,1:ke_tot-1))-... K{I"2c  
                                                        (hz(2:ie_tot,:,2:ke_tot)-hz(1:ie_tot-1,:,2:ke_tot)))./delta; 3t9+YdNKU  
    ey(2:ie_tot,:,2:ke_tot)=C3ey(2:ie_tot,:,2:ke_tot).*ey(2:ie_tot,:,2:ke_tot)+... vqT)=ZC1  
                            C4ey(2:ie_tot,:,2:ke_tot).*(C5ey(2:ie_tot,:,2:ke_tot).*dy(2:ie_tot,:,2:ke_tot)-... j7g>r/1eE  
                                                        C6ey(2:ie_tot,:,2:ke_tot).*dstore(2:ie_tot,:,2:ke_tot)); ]#shuZ##>0  
    dstore=dz; $> QJ%v9+  
    dz(2:ie_tot,2:je_tot,:)=C1ez(2:ie_tot,2:je_tot,:).*  dz(2:ie_tot,2:je_tot,:)+... uaP5(hUI  
                            C2ez(2:ie_tot,2:je_to .. ?H_@/?  
\Kui`X  

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

哎，连个回答俺的也没有 EY)?hJS,  

这个程序我调过，没问题的。这个加的是软源，Ez分量是上一步的值和高斯脉冲叠加的结果

高手们都说说 I\[z(CHg@  

好的，谢谢啊，可让俺放心了

这个图应该是高斯脉冲的形式  但你的不是 qIh9? |`U  
检查一下你的空间步长设定是不是太短了 也可以适当的增加时间步 让程序多迭代一些时间步 看看有没有什么变化

同意5楼的看法，可能是时间步长太小了，不过如果增大时间步长的间隔有可能会发散，造成不稳定，所以妥帖的办法是增大时间步数。

其实这个结果是对了 +WguWLO"  
原因解释如下： Z!4
B=?(  
你采用的一般高斯源，它主要能量集中在低频，即直流分量很大，观察曲线后面部分是直线，是直流源的影响，相当于在观察处的静场值。

如果采用一阶以上的Gauss源或调制高斯源，就没有直流分量，仿真出来没有这个现象