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

图片:cuwu.jpg

o1-_BlZ  
利用3D-upml 仿真了自由空间，加的激励为高斯脉冲ez(is,js,ks)=ez(is,js,ks)+exp(-((n-45)^2/15^2)); 10bv%ZX7  
然后我就检测了一下ez(is,js,ks)随时间的变化图，结果是图中所示，不再是高斯脉冲形式 _c}# f\ +_  
这样是不是不对啊。 程序师用的下面的那个，基本没改，就是改变了一个激励源。 +AFBTJ  
%*********************************************************************** +/" \.wYv  
%     3-D FDTD code with UPML absorbing boundary conditions |.-Muv  
%*********************************************************************** *55unc  
% NYzBfL x  
%     Program author: Keely J. Willis, Graduate Student /a6i`  
%                     UW Computational Electromagnetics Laboratory I<+:Ho=6  
%                           Director: Susan C. Hagness iqN?'8  
%                     Department of Electrical and Computer Engineering BKgCuz:y  
%                     University of Wisconsin-Madison QFgKEUNgl  
%                     1415 Engineering Drive vTIRydg2b  
%                     Madison, WI 53706-1691 \m:('^\6o  
%                     kjwillis@wisc.edu 1jaK N*  
% OGG9f??  
%     Copyright 2005 3.KNAObO  
% 4Tb"+Y
}  
%     This MATLAB M-file implements the finite-difference time-domain wt
i  
%     solution of Maxwell's curl equations over a three-dimensional $\M];S=CY  
%     Cartesian space lattice comprised of uniform cubic grid cells. wrsr
U  
%     i)$<j!L  
%     The dimensions of the computational domain are 8.2 cm ${gO=Z  
%     (x-direction), 3.4 cm (y-direction), and 3.2 cm (z-direction).  n9-WZsc1  
%     The grid is terminated with UPML absorbing boundary conditions. 8NTE`l=>/  
% 7<Y aw,G  
%     An electric current source comprised of two collinear Jz components /w2-Pgm-[\  
%     (realizing a Hertzian dipole) excites a radially propagating wave.  4U u`1gtz  
%     The current source is located in the center of the grid.  The Zq5~M bldh  
%     source waveform is a differentiated Gaussian pulse given by S6fbwZZMG  
%          J(t)=J0*(t-t0)*exp(-(t-t0)^2/tau^2), =1<v1s|)q  
%     where tau=50 ps.  The FWHM spectral bandwidth of this zero-dc- QbY@{"" `  
%     content pulse is approximately 7 GHz. The grid resolution <Pi#-r.,  
%     (dx = 2 mm) was chosen to provide at least 10 samples per fVdu9 l  
%     wavelength up through 15 GHz. o^r\7g6\  
% 0sB[]E|7[s  
%     To execute this M-file, type "fdtd3D_UPML" at the MATLAB prompt.  wyXQP+9G  
% Dv&K3^~Rfb  
%     This code has been tested in the following Matlab environments: p%K(dA  
%     Matlab version 6.1.0.450 Release 12.1 (May 18, 2001) ,,BWWFg~  
%     Matlab version 6.5.1.199709 Release 13 Service Pack 1 (August 4, 2003) [khXAf1{Q  
%     Matlab version 7.0.0.19920 R14 (May 6, 2004) h 9}x6t,  
%     Matlab version 7.0.1.24704 R14 Service Pack 1 (September 13, 2004) i?7?I  
%     Matlab version 7.0.4.365 R14 Service Pack 2 (January 29, 2005) ~qK/w0=j  
% : LT'#Q8  
%     Note: if you are using Matlab version 6.x, you may wish to make Aq\K N.  
%     one or more of the following modifications to this code: Eh$1piJG  
%       --uncomment line numbers 485 and 486 Ds#BfP7a  
%       --comment out line numbers 552 and 561 3Vak C  
% KKWvV4u  
%*********************************************************************** b[:{\!I  
clear k|U2Mp  
%*********************************************************************** aK7}}  
%     Fundamental constants 2.MY8}&WBu  
%*********************************************************************** Kx?8HA[5  
cc=2.99792458e8; '}"&JO~vPj  
muz=4.0*pi*1.0e-7; {<?8Y  
epsz=1.0/(cc*cc*muz); e^$JGh2  
etaz=sqrt(muz/epsz); 8pZOgh  
%*********************************************************************** ;k,@^f8  
%     Material parameters 'K#ndCGJ$  
%*********************************************************************** y*p02\)  
mur=1.0; JV_VM{w{K  
epsr=1.0; c5:X$k\  
eta=etaz*sqrt(mur/epsr); 0sTR`Xk  
%*********************************************************************** (L(n%  
%     Grid parameters Xg*](>/\,  
% 8 VhU)fY  
%     Each grid size variable name describes the number of sampled points N,3iSH=cN[  
%     for a particular field component in the direction of that component. ?0?3yD-!9  
%     Additionally, the variable names indicate the region of the grid \N$)Q.M  
%     for which the dimension is relevant.  For example, ie_tot is the E>`|?DE@  
%     number of sample points of Ex along the x-axis in the total <1 ;pyw y  
%     computational grid, and jh_bc is the number of sample points of Hy NB+/S;`  
%     along the y-axis in the y-normal UPML regions. B&6lG!K'?  
% LWhPd\  
%*********************************************************************** !e*T. 1Kz  
ie=60;          % Size of main grid <XN=v!2;  
je=60; |=MhI5gsx  
ke=60; G\B+bBz  
ih=ie+1; (t@:dW  
jh=je+1;   Q|e-)FS)  
kh=ke+1;   FZLx.3k4  
upml=10;        % Thickness of PML boundaries a,r B7aD  
ih_bc=upml+1; II!~"-WH  
jh_bc=upml+1; Qkhor-f0  
kh_bc=upml+1; 8_"NF%%(n  
ie_tot=ie+2*upml;          % Size of total computational domain .t/@d(R  
je_tot=je+2*upml;        qI${7  
ke_tot=ke+2*upml;        ww #kc!'  
ih_tot=ie_tot+1; mrr~#Bb>  
jh_tot=je_tot+1;          BRM!g9  
kh_tot=ke_tot+1;          P+wpX  
is=round(ih_tot/2);         % Location of z-directed current source `!AI:c*3p1  
js=round(jh_tot/2); +T8MQ[(4  
ks=round(ke_tot/2); gga}mqMv=  
%*********************************************************************** VqxK5  
%     Fundamental grid parameters yc`*zLWh  
%*********************************************************************** > >KCd  
delta=0.002;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \ Ce
*5
h  
dt=delta*sqrt(epsr*mur)/(2.0*cc); }}D32TVN  
nmax=100;
b&dv("e 4  
%*********************************************************************** {8oGWQgrj  
%     Differentiated Gaussian pulse excitation HrfS^B  
%*********************************************************************** OA(.&5]  
rtau=50.0e-12; 't5`Ni  
tau=rtau/dt; , 2xv  
ndelay=3*tau; _l"nwE
s  
J0=1.0*epsz; =O-irGms*  
%*********************************************************************** Z!7
xRy  
%     Initialize field arrays ?~!9\dek,  
%*********************************************************************** R<(xWH  
ex=zeros(ie_tot,jh_tot,kh_tot); e<[ ] W4"A  
ey=zeros(ih_tot,je_tot,kh_tot); To5hVL<Ex"  
ez=zeros(ih_tot,jh_tot,ke_tot); v+8Ybq  
dx=zeros(ie_tot,jh_tot,kh_tot); tC5-^5[y  
dy=zeros(ih_tot,je_tot,kh_tot); u05Yy&(f  
dz=zeros(ih_tot,jh_tot,ke_tot); t,IOq[Vtk  
hx=zeros(ih_tot,je_tot,ke_tot); ra>2<  
hy=zeros(ie_tot,jh_tot,ke_tot); PB?2{Cj
 
hz=zeros(ie_tot,je_tot,kh_tot); -V;BkE
76  
bx=zeros(ih_tot,je_tot,ke_tot); 
=I
@I  
by=zeros(ie_tot,jh_tot,ke_tot); v^vi *c  
bz=zeros(ie_tot,je_tot,kh_tot); =0!j"z=  
%*********************************************************************** aT(_c/t.  
%     Initialize update coefficient arrays wy0?*)~  
%*********************************************************************** D W^Zuu/)  
C1ex=zeros(size(ex)); C/'w  
C2ex=zeros(size(ex)); S(?A3 H  
C3ex=zeros(size(ex)); VpSpj/\m)'  
C4ex=zeros(size(ex)); B?- poB&  
C5ex=zeros(size(ex)); G$%F`R[  
C6ex=zeros(size(ex)); zn7)>cQ905  
C1ey=zeros(size(ey)); A.dbb'^  
C2ey=zeros(size(ey)); ,?k1if(0[  
C3ey=zeros(size(ey)); k~ByICE  
C4ey=zeros(size(ey)); C4P<GtR9  
C5ey=zeros(size(ey)); T~(Sc'8  
C6ey=zeros(size(ey)); :O$bsw:3w<  
C1ez=zeros(size(ez)); RHMXPsj  
C2ez=zeros(size(ez)); RjVmHhX
 
C3ez=zeros(size(ez)); X2rKH$<g  
C4ez=zeros(size(ez)); o:fe`#t  
C5ez=zeros(size(ez)); [.1MElM  
C6ez=zeros(size(ez)); ;i'[c`  
D1hx=zeros(size(hx)); [/%N2mj  
D2hx=zeros(size(hx)); G\TO]c  
D3hx=zeros(size(hx)); [~2imS  
D4hx=zeros(size(hx)); j49Uj}:j  
D5hx=zeros(size(hx)); d7 H*F  
D6hx=zeros(size(hx)); /XEW]/4  
D1hy=zeros(size(hy)); x3Y)l1gh  
D2hy=zeros(size(hy)); -Ou.C7ol  
D3hy=zeros(size(hy)); kS:#|yY8%  
D4hy=zeros(size(hy)); Dfa3&##{  
D5hy=zeros(size(hy)); c38XM]Jeq  
D6hy=zeros(size(hy)); ]z/R?SM  
D1hz=zeros(size(hz)); stBe ^C  
D2hz=zeros(size(hz)); lg~7[=%k#  
D3hz=zeros(size(hz)); G{E`5KIvm  
D4hz=zeros(size(hz)); v{fcQb  
D5hz=zeros(size(hz)); 2wHbhW[  
D6hz=zeros(size(hz)); X<P <-e9  
%*********************************************************************** UL{J%Ze=~  
%     Update coefficients, as described in Section 7.8.2. |E.BGdS  
% \r[u>7I  
%     In order to simplify the update equations used in the time-stepping F_jHi0A  
%     loop, we implement UPML update equations throughout the entire AyOibnoZ2E  
%     grid.  In the main grid, the electric-field update coefficients of $2L6:&.P,  
%     Equations 7.91a-f and the correponding magnetic field update 1{ %y(?`  
%     coefficients extracted from Equations 7.89 and 7.90 are simplified 3m`>D e  
%     for the main grid (free space) and calculated below. ,<r&] eC  
% )AQ
^PBwp  
%*********************************************************************** DQm%=ON7  
C1=1.0;
c$%*p (zY  
C2=dt; 9S*"={}%  
C3=1.0; $[n:IDa*@1  
C4=1.0/2.0/epsr/epsr/epsz/epsz; *{!Y_FrL  
C5=2.0*epsr*epsz; OmO#} k<  
C6=2.0*epsr*epsz; (rkg0  
D1=1.0; I4{xQI  
D2=dt; ~~Ezt*lH  
D3=1.0; HOF$(86zqA  
D4=1.0/2.0/epsr/epsz/mur/muz; y{>f^S<  
D5=2.0*epsr*epsz; Jt@lH  
D6=2.0*epsr*epsz; #c>GjUJ.w  
%*********************************************************************** W%-XN  
%     Initialize main grid update coefficients Xa?O)Bq.  
%*********************************************************************** |f#hGk6  
C1ex(:,jh_bc:jh_tot-upml,:)=C1;     PD-&(ka.  
C2ex(:,jh_bc:jh_tot-upml,:)=C2; hN &?x5aC>  
C3ex(:,:,kh_bc:kh_tot-upml)=C3; a[(OeVQ5  
C4ex(:,:,kh_bc:kh_tot-upml)=C4; *_
o(~5w-K  
C5ex(ih_bc:ie_tot-upml,:,:)=C5; .'gm2  
C6ex(ih_bc:ie_tot-upml,:,:)=C6; I}3F'}JV<  
C1ey(:,:,kh_bc:kh_tot-upml)=C1; )J}v.8  
C2ey(:,:,kh_bc:kh_tot-upml)=C2; e12QYoh  
C3ey(ih_bc:ih_tot-upml,:,:)=C3; UI+6\ 3  
C4ey(ih_bc:ih_tot-upml,:,:)=C4; '#Au~5  
C5ey(:,jh_bc:je_tot-upml,:)=C5; g-~ _gt7  
C6ey(:,jh_bc:je_tot-upml,:)=C6; bYnq,JRA  
C1ez(ih_bc:ih_tot-upml,:,:)=C1; D5D *$IC  
C2ez(ih_bc:ih_tot-upml,:,:)=C2; .
Dr!\.hL  
C3ez(:,jh_bc:jh_tot-upml,:)=C3; c{BAQZVc  
C4ez(:,jh_bc:jh_tot-upml,:)=C4; 7MLLx#U  
C5ez(:,:,kh_bc:ke_tot-upml)=C5; "J1A9|  
C6ez(:,:,kh_bc:ke_tot-upml)=C6; a3tcLd|7J  
D1hx(:,jh_bc:je_tot-upml,:)=D1; 7RL J  
D2hx(:,jh_bc:je_tot-upml,:)=D2; 2!Dz
9m3  
D3hx(:,:,kh_bc:ke_tot-upml)=D3; -HG.GA  
D4hx(:,:,kh_bc:ke_tot-upml)=D4; VTM* 1uXS>  
D5hx(ih_bc:ih_tot-upml,:,:)=D5; _9 ]:0bDUo  
D6hx(ih_bc:ih_tot-upml,:,:)=D6; 7GYf#} N  
D1hy(:,:,kh_bc:ke_tot-upml)=D1; j)?M  
D2hy(:,:,kh_bc:ke_tot-upml)=D2; \}s/<Q  
D3hy(ih_bc:ie_tot-upml,:,:)=D3; % D  
D4hy(ih_bc:ie_tot-upml,:,:)=D4; Wye* ~t  
D5hy(:,jh_bc:jh_tot-upml,:)=D5; Pc`d]*BYi  
D6hy(:,jh_bc:jh_tot-upml,:)=D6; pOc2V  
D1hz(ih_bc:ie_tot-upml,:,:)=D1; SG&,o=I$  
D2hz(ih_bc:ie_tot-upml,:,:)=D2; QLWnP-  
D3hz(:,jh_bc:je_tot-upml,:)=D3; |Ev|A9J!  
D4hz(:,jh_bc:je_tot-upml,:)=D4; zVq!M-e  
D5hz(:,:,kh_bc:kh_tot-upml)=D5; Pwl*5/l  
D6hz(:,:,kh_bc:kh_tot-upml)=D6; yu6{6[  
%*********************************************************************** qTr P@F4`g  
%     Fill in PML regions 30!DraW8  
% FklR!*oL,)  
%     PML theory describes a continuous grading of the material properties yx :^*/  
%     over the PML region.  In the FDTD grid it is necessary to discretize vSH,fS-n  
%     the grading by averaging the material properties over a grid cell 8(L$a1#5W  
%     width centered on each field component.  As an example of the Z~~6y6p  
%     implementation of this averaging, we take the integral of the }kqh[`:  
%     continuous sigma(x) in the PML region oLT#'42+H  
%   ,X4+i8Yc  
%         sigma_i = integral(sigma(x))/delta .d]/:T -0  
%   W2 -%/  
%     where the integral is over a single grid cell width in x, and is 6]?mjG6  
%     bounded by x1 and x2.  Applying this to the polynomial grading of oS fr5 i  
%     Equation 7.60a produces %N*[{j= ^  
% (Xh<F  
%         sigma_i = (x2^(m+1)-x1^(m+1))*sigmam/(delta*(m+1)*d^m) c$Kc,`2m
7  
% sFTAE1|  
%     where sigmam is the maximum value of sigma as described by Equation S\g9@g.  
%     7.62. WiS3W;  
%         lFjz*g2'  
%     The definitions of x1 and x2 depend on the position of the field 7~e,"^>T
 
%     component within the grid cell.  We have either OlOOg  
% 4,I,f>V  
%         x1 = (i-0.5)*delta,  x2 = (i+0.5)*delta \9)5b8  
%  XB7Aa)  
%     or ^v'kEsE^*  
%  k&ci5MpN  
%         x1 = (i)*delta,      x2 = (i+1)*delta _x,X0ncv]@  
% ny5P*yWEh  
%     where i varies over the PML region. .h-mFcjy  
%  ]#)(D-i  
%*********************************************************************** ga5Q  
rmax=exp(-16);  %desired reflection error, designated as R(0) in Equation 7.62
@:C)^f"  
orderbc=4;      %order of the polynomial grading, designated as m in Equation 7.60a,b }qn>#ETi  
%   x-varying material properties 4>*=q*<V5E  
delbc=upml*delta; Zv;nY7B  
sigmam=-log(rmax)*(orderbc+1.0)/(2.0*eta*delbc); M:/NW-:  
sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1.0)); 79v+ze  
kmax=1; &<VU}c^!  
kfactor=(kmax-1.0)/delta/(orderbc+1.0)/delbc^orderbc; s6,~JF^  
for i=1:upml {dpC;jsW1  
    NPT-d  
    % Coefficients for field components in the center of the grid cell .\R9tt}  
    x1=(upml-i+1)*delta; HAxLYun(3w  
    x2=(upml-i)*delta; '~D4%WKT  
    sigma=sigfactor*(x1^(orderbc+1)-x2^(orderbc+1)); `Nx@MPo  
    ki=1+kfactor*(x1^(orderbc+1)-x2^(orderbc+1)); P%aqY~yF3  
    facm=(2*epsr*epsz*ki-sigma*dt); Kjd3!%4mB  
    facp=(2*epsr*epsz*ki+sigma*dt); i%K6<1R;y{  
    C5ex(i,:,:)=facp; _QL|pLf-  
    C5ex(ie_tot-i+1,:,:)=facp; :?6HG_9X  
    C6ex(i,:,:)=facm; fEHFlgN3Ap  
    C6ex(ie_tot-i+1,:,:)=facm; >G6kF!V  
    D1hz(i,:,:)=facm/facp; K%v:giN$l`  
    D1hz(ie_tot-i+1,:,:)=facm/facp; D&%8JL  
    D2hz(i,:,:)=2.0*epsr*epsz*dt/facp; GY%9V5GB  
    D2hz(ie_tot-i+1,:,:)=2.0*epsr*epsz*dt/facp; _z@/~M(  
    D3hy(i,:,:)=facm/facp; NfV|c~
?d  
    D3hy(ie_tot-i+1,:,:)=facm/facp; TEz;:*,CG  
    D4hy(i,:,:)=1.0/facp/mur/muz; {N4 'g_  
    D4hy(ie_tot-i+1,:,:)=1.0/facp/mur/muz; 23gN;eD+m6  
    % Coefficients for field components on the grid cell boundary P0l fK}  
    x1=(upml-i+1.5)*delta; qVC+q8  
    x2=(upml-i+0.5)*delta; ~T_|?lU`R  
    sigma=sigfactor*(x1^(orderbc+1)-x2^(orderbc+1)); [ohLG_9  
    ki=1.0+kfactor*(x1^(orderbc+1)-x2^(orderbc+1)); l=CAr  
    facm=(2.0*epsr*epsz*ki-sigma*dt); r1L@p[>  
    facp=(2.0*epsr*epsz*ki+sigma*dt); r%U6,7d=)  
    C1ez(i,:,:)=facm/facp; T+Z[&|  
    C1ez(ih_tot-i+1,:,:)=facm/facp; %R0 Wq4}  
    C2ez(i,:,:)=2.0*epsr*epsz*dt/facp; @]l|-xGCWn  
    C2ez(ih_tot-i+1,:,:)=2.0*epsr*epsz*dt/facp; Hd~g\  
    C3ey(i,:,:)=facm/facp; u#76w74  
    C3ey(ih_tot-i+1,:,:)=facm/facp; t*IePz]/  
    C4ey(i,:,:)=1.0/facp/epsr/epsz; q<&1,^A  
    C4ey(ih_tot-i+1,:,:)=1.0/facp/epsr/epsz; wQ+pVu?6_  
    D5hx(i,:,:)=facp; {A0jkU  
    D5hx(ih_tot-i+1,:,:)=facp; fDy*dp4z  
    D6hx(i,:,:)=facm; |Ea%nghl  
    D6hx(ih_tot-i+1,:,:)=facm; d&+]@ Ii
 
    qLEYBv-3  
end QLY;@-jF$  
%   PEC walls ^Arv6kD,  
C1ez(1,:,:)=-1.0; Nny*C`uDF  
C1ez(ih_tot,:,:)=-1.0; FK^xZ?G  
C2ez(1,:,:)=0.0; *9
\j1Nd  
C2ez(ih_tot,:,:)=0.0; Cn~VJ,l g  
C3ey(1,:,:)=-1.0; @xWWN  
C3ey(ih_tot,:,:)=-1.0; @_ %RQO_X  
C4ey(1,:,:)=0.0; ?'> .>  
C4ey(ih_tot,:,:)=0.0; [c,V=:Cq  
%   y-varying material properties A9xeOy8e  
delbc=upml*delta; MQoA\  
sigmam=-log(rmax)*epsr*epsz*cc*(orderbc+1.0)/(2.0*delbc); m_)-  
sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1.0)); YV.' L  
kmax=1.0;  d$$5&a  
kfactor=(kmax-1.0)/delta/(orderbc+1.0)/delbc^orderbc; `UsJaoR#f  
for j=1:upml jIs>>  
    7{m>W!  
    % Coefficients for field components in the center of the grid cell 2;v:Z^&  
    y1=(upml-j+1)*delta; D6bYg `  
    y2=(upml-j)*delta; 'nTlCYT  
    sigma=sigfactor*(y1^(orderbc+1)-y2^(orderbc+1)); z!g$#hmL>  
    ki=1+kfactor*(y1^(orderbc+1)-y2^(orderbc+1)); w6vbYPCN  
    facm=(2*epsr*epsz*ki-sigma*dt); /e2zH  
    facp=(2*epsr*epsz*ki+sigma*dt); w-K A~  
    ]?y~;-^  
    C5ey(:,j,:)=facp; 1-y8Hy_a2  
    C5ey(:,je_tot-j+1,:)=facp; rCPIz<  
    C6ey(:,j,:)=facm; dA)T>  
    C6ey(:,je_tot-j+1,:)=facm; wH~A> 4*(  
    D1hx(:,j,:)=facm/facp; ?X|)0o  
    D1hx(:,je_tot-j+1,:)=facm/facp; ]Y[N=
G  
    D2hx(:,j,:)=2*epsr*epsz*dt/facp; ##jJaSxG  
    D2hx(:,je_tot-j+1,:)=2*epsr*epsz*dt/facp; 2I
B{FO/  
    D3hz(:,j,:)=facm/facp; xuXPVJdi  
    D3hz(:,je_tot-j+1,:)=facm/facp; (<Cq_Kw  
    D4hz(:,j,:)=1/facp/mur/muz; !n-Sh<8  
    D4hz(:,je_tot-j+1,:)=1/facp/mur/muz; j\ y!  
    ]o] VS  
    % Coefficients for field components on the grid cell boundary vb>F)X?b_  
    y1=(upml-j+1.5)*delta; v9f+ {Y%-  
    y2=(upml-j+0.5)*delta; +=($mcw#[  
    sigma=sigfactor*(y1^(orderbc+1)-y2^(orderbc+1)); s5*4<VxQN.  
    ki=1+kfactor*(y1^(orderbc+1)-y2^(orderbc+1)); '[$KG  
    facm=(2*epsr*epsz*ki-sigma*dt); )g@+ MR  
    facp=(2*epsr*epsz*ki+sigma*dt);    @[r[l#4yUi  
     BN9e S   
    C1ex(:,j,:)=facm/facp;
7KIekL  
    C1ex(:,jh_tot-j+1,:)=facm/facp; #*iUZo  
    C2ex(:,j,:)=2*epsr*epsz*dt/facp; c(Dp`f,  
    C2ex(:,jh_tot-j+1,:)=2*epsr*epsz*dt/facp; Bp^LLH  
    C3ez(:,j,:)=facm/facp; v'hc-Q9+>  
    C3ez(:,jh_tot-j+1,:)=facm/facp; VIF43/>(  
    C4ez(:,j,:)=1/facp/epsr/epsz; ;%n'k  
    C4ez(:,jh_tot-j+1,:)=1/facp/epsr/epsz;   FyEKqYl  
    D5hy(:,j,:)=facp; h@ lz  
    D5hy(:,jh_tot-j+1,:)=facp; yj:@Fg-3g  
    D6hy(:,j,:)=facm; YS|Dw'%g /  
    D6hy(:,jh_tot-j+1,:)=facm; &~_F2]oM  
end Mq0MtC6-  
%   PEC walls u> {aF{  
C1ex(:,1,:)=-1; _[6sr7H!  
C1ex(:,jh_tot,:)=-1; 0|AgmW_7 .  
C2ex(:,1,:)=0; E)*ht;u  
C2ex(:,jh_tot,:)=0; l[E^nh>  
C3ez(:,1,:)=-1; 01mu6)  
C3ez(:,jh_tot,:)=-1; fu!T4{2  
C4ez(:,1,:)=0; }b2YX+/e$f  
C4ez(:,jh_tot,:)=0;   Jqxd92 bI  
%   z-varying material properties Biv)s@"f-Q  
delbc=upml*delta; \TP$2i%W  
sigmam=-log(rmax)*epsr*epsz*cc*(orderbc+1)/(2*delbc); 9`f@"%h  
sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1)); pT,8E(*l2  
kmax=1; `3\aX|4@  
kfactor=(kmax-1)/delta/(orderbc+1)/delbc^orderbc; _#{*I(l  
for k=1:upml kK75(x  
    % Coefficients for field components in the center of the grid cell ys`-QlkB  
    z1=(upml-k+1)*delta; n-9xfn0U~#  
    z2=(upml-k)*delta; [<XYU,{R  
    sigma=sigfactor*(z1^(orderbc+1)-z2^(orderbc+1)); HT.,BF  
    ki=1+kfactor*(z1^(orderbc+1)-z2^(orderbc+1)); sa.H,<;  
    facm=(2*epsr*epsz*ki-sigma*dt); do8[wej<:  
    facp=(2*epsr*epsz*ki+sigma*dt); }TTghE!  
    O\<zQ2m  
    C5ez(:,:,k)=facp; cSPQ NYU:  
    C5ez(:,:,ke_tot-k+1)=facp; y.Z_\@  
    C6ez(:,:,k)=facm; 3q%z
 
    C6ez(:,:,ke_tot-k+1)=facm; Q/|.=:~FO  
    D1hy(:,:,k)=facm/facp; 9QU\J0c/  
    D1hy(:,:,ke_tot-k+1)=facm/facp; z6`0Uv~  
    D2hy(:,:,k)=2*epsr*epsz*dt/facp; F3bTFFt  
    D2hy(:,:,ke_tot-k+1)=2*epsr*epsz*dt/facp; m7k }k)  
    D3hx(:,:,k)=facm/facp; 9^/Y7Wp/@  
    D3hx(:,:,ke_tot-k+1)=facm/facp; fw&*;az  
    D4hx(:,:,k)=1/facp/mur/muz; Z1fY' f  
    D4hx(:,:,ke_tot-k+1)=1/facp/mur/muz; 9dNB_  
    Wc@ ,#v  
    % Coefficients for field components on the grid cell boundary 7T/BzXr,B  
    z1=(upml-k+1.5)*delta; b+~_/;Y9  
    z2=(upml-k+0.5)*delta; ~xqiasE#K  
    sigma=sigfactor*(z1^(orderbc+1)-z2^(orderbc+1)); qm=U<'b^  
    ki=1+kfactor*(z1^(orderbc+1)-z2^(orderbc+1)); O
i\ s  
    facm=(2*epsr*epsz*ki-sigma*dt); KS*,'hvY  
    facp=(2*epsr*epsz*ki+sigma*dt); vEI{AmogRx  
    W`x.qumN
 
    C1ey(:,:,k)=facm/facp; |^1g*fy?  
    C1ey(:,:,kh_tot-k+1)=facm/facp; ]Za[]E8MD  
    C2ey(:,:,k)=2*epsr*epsz*dt/facp; dfFw6R  
    C2ey(:,:,kh_tot-k+1)=2*epsr*epsz*dt/facp; zQ+Mu^|u+  
    C3ex(:,:,k)=facm/facp; jRg/N_2'2  
    C3ex(:,:,kh_tot-k+1)=facm/facp; cX9o'e:C  
    C4ex(:,:,k)=1/facp/epsr/epsz; y]Nk^ga:U6  
    C4ex(:,:,kh_tot-k+1)=1/facp/epsr/epsz; qt L]x -O  
    D5hz(:,:,k)=facp; N}#Rw2Vl  
    D5hz(:,:,kh_tot-k+1)=facp; y,:WLk~  
    D6hz(:,:,k)=facm; 1'f_C<.0  
    D6hz(:,:,kh_tot-k+1)=facm; uX-^9t  
end )%^l+w+&  
%   PEC walls =w3A{h"^  
C1ey(:,:,1)=-1; RI*n]HNgy+  
C1ey(:,:,kh_tot)=-1; =AO (  
C2ey(:,:,1)=0; 6Amt75RY  
C2ey(:,:,kh_tot)=0; f&CQn.K"  
C3ex(:,:,1)=-1; \ ITd\)F%N  
C3ex(:,:,kh_tot)=-1; 1%_RXQVG  
C4ex(:,:,1)=0; >5t! Xt  
C4ex(:,:,kh_tot)=0; !yv>e7g^  
%figure AFi_P\X  
%set(gcf,'DoubleBuffer','on') 4$iS@o|  
%*********************************************************************** =DdPwr 0Op  
%     Begin time stepping loop @XJ7ff&  
%*********************************************************************** d92Z;FWb  
d0=zeros(1,nmax);d1=zeros(1,nmax);d2=zeros(1,nmax); W.^zN'a  
for n=1:nmax VJ\qp%  
    K+)3 LR^  
    % Update magnetic field ~u%$ 9IhM  
    bstore=bx; l73% y  
    bx(2:ie_tot,:,:)=D1hx(2:ie_tot,:,:).*  bx(2:ie_tot,:,:)-... B.y}S  
                     D2hx(2:ie_tot,:,:).*((ez(2:ie_tot,2:jh_tot,:)-ez(2:ie_tot,1:je_tot,:))-... @aQ:3/  
                                          (ey(2:ie_tot,:,2:kh_tot)-ey(2:ie_tot,:,1:ke_tot)))./delta; 4w+AOWjd  
    hx(2:ie_tot,:,:)= D3hx(2:ie_tot,:,:).*hx(2:ie_tot,:,:)+... [7}3k?42X  
                      D4hx(2:ie_tot,:,:).*(D5hx(2:ie_tot,:,:).*bx(2:ie_tot,:,:)-... e:HORc~U  
                                           D6hx(2:ie_tot,:,:).*bstore(2:ie_tot,:,:)); }.o.*N  
    bstore=by; zr!7*, p  
    by(:,2:je_tot,:)=D1hy(:,2:je_tot,:).*  by(:,2:je_tot,:)-... rN9qH  
                     D2hy(:,2:je_tot,:).*((ex(:,2:je_tot,2:kh_tot)-ex(:,2:je_tot,1:ke_tot))-... 9D14/9*(dU  
                                          (ez(2:ih_tot,2:je_tot,:)-ez(1:ie_tot,2:je_tot,:)))./delta; YC1Bgz  
    hy(:,2:je_tot,:)= D3hy(:,2:je_tot,:).*hy(:,2:je_tot,:)+... _.9 5>`  
                      D4hy(:,2:je_tot,:).*(D5hy(:,2:je_tot,:).*by(:,2:je_tot,:)-... U,!qNi}  
                                           D6hy(:,2:je_tot,:).*bstore(:,2:je_tot,:)); j4;^5 Dy^  
    bstore=bz; g(pr.Dw6  
    bz(:,:,2:ke_tot)=D1hz(:,:,2:ke_tot).*  bz(:,:,2:ke_tot)-... dU9;sx  
                     D2hz(:,:,2:ke_tot).*((ey(2:ih_tot,:,2:ke_tot)-ey(1:ie_tot,:,2:ke_tot))-... 9`Qa/Y!  
                                          (ex(:,2:jh_tot,2:ke_tot)-ex(:,1:je_tot,2:ke_tot)))./delta; CcUF)$kz  
    hz(:,:,2:ke_tot)= D3hz(:,:,2:ke_tot).*hz(:,:,2:ke_tot)+... s?HK2b^;D  
                      D4hz(:,:,2:ke_tot).*(D5hz(:,:,2:ke_tot).*bz(:,:,2:ke_tot)-... E n7~wKF  
                                           D6hz(:,:,2:ke_tot).*bstore(:,:,2:ke_tot)); GTLS0l)  
    .F,l>wUNe  
    % Update electric field Tw';;euw  
    dstore=dx; XncX2E4E  
    dx(:,2:je_tot,2:ke_tot)=C1ex(:,2:je_tot,2:ke_tot).*  dx(:,2:je_tot,2:ke_tot)+... deAV:c  
                            C2ex(:,2:je_tot,2:ke_tot).*((hz(:,2:je_tot,2:ke_tot)-hz(:,1:je_tot-1,2:ke_tot))-... X|\`\[  
                                                        (hy(:,2:je_tot,2:ke_tot)-hy(:,2:je_tot,1:ke_tot-1)))./delta; c>$d!IKCL  
    ex(:,2:je_tot,2:ke_tot)=C3ex(:,2:je_tot,2:ke_tot).*ex(:,2:je_tot,2:ke_tot)+... (m'-1wX.  
                            C4ex(:,2:je_tot,2:ke_tot).*(C5ex(:,2:je_tot,2:ke_tot).*dx(:,2:je_tot,2:ke_tot)-... x\'3UKQP+^  
                                                        C6ex(:,2:je_tot,2:ke_tot).*dstore(:,2:je_tot,2:ke_tot)); _".h(  
    dstore=dy; AZ(zM.y!#_  
    dy(2:ie_tot,:,2:ke_tot)=C1ey(2:ie_tot,:,2:ke_tot).*  dy(2:ie_tot,:,2:ke_tot)+... wn@~80)$  
                            C2ey(2:ie_tot,:,2:ke_tot).*((hx(2:ie_tot,:,2:ke_tot)-hx(2:ie_tot,:,1:ke_tot-1))-... 6?u`u t  
                                                        (hz(2:ie_tot,:,2:ke_tot)-hz(1:ie_tot-1,:,2:ke_tot)))./delta; ou-#+Sdd  
    ey(2:ie_tot,:,2:ke_tot)=C3ey(2:ie_tot,:,2:ke_tot).*ey(2:ie_tot,:,2:ke_tot)+... Z7bJ<TpZ  
                            C4ey(2:ie_tot,:,2:ke_tot).*(C5ey(2:ie_tot,:,2:ke_tot).*dy(2:ie_tot,:,2:ke_tot)-... Fj`k3~tUw  
                                                        C6ey(2:ie_tot,:,2:ke_tot).*dstore(2:ie_tot,:,2:ke_tot)); s'yR2JYv  
    dstore=dz; `<g]p-=":  
    dz(2:ie_tot,2:je_tot,:)=C1ez(2:ie_tot,2:je_tot,:).*  dz(2:ie_tot,2:je_tot,:)+... X
MS:F]HN  
                            C2ez(2:ie_tot,2:je_to .. r
Ql9SUs  
)(,O~w  

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

哎，连个回答俺的也没有 ^oav-R&  

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

高手们都说说 +kd1q  

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

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

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

其实这个结果是对了 4c< s"2F  
原因解释如下： uem-fTG  
你采用的一般高斯源，它主要能量集中在低频，即直流分量很大，观察曲线后面部分是直线，是直流源的影响，相当于在观察处的静场值。

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