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

图片:cuwu.jpg

0gn@h/F2%  
利用3D-upml 仿真了自由空间，加的激励为高斯脉冲ez(is,js,ks)=ez(is,js,ks)+exp(-((n-45)^2/15^2)); )){xlFA}  
然后我就检测了一下ez(is,js,ks)随时间的变化图，结果是图中所示，不再是高斯脉冲形式 sIl33kmv  
这样是不是不对啊。 程序师用的下面的那个，基本没改，就是改变了一个激励源。 6UE(f@  
%*********************************************************************** -[Qvg49jy  
%     3-D FDTD code with UPML absorbing boundary conditions PjIeZ&p  
%*********************************************************************** ZCg`z
  
% sgr=w+",Q  
%     Program author: Keely J. Willis, Graduate Student 3
},Zlu  
%                     UW Computational Electromagnetics Laboratory 2 /y}a#s  
%                           Director: Susan C. Hagness I=O
y-  
%                     Department of Electrical and Computer Engineering  pAu72O?  
%                     University of Wisconsin-Madison ^MmC$U^n  
%                     1415 Engineering Drive )=
Q)BN[  
%                     Madison, WI 53706-1691 0~&"  
%                     kjwillis@wisc.edu mga6[E<  
% e0;  
%     Copyright 2005 %o}(sShS  
% =_[Z W  
%     This MATLAB M-file implements the finite-difference time-domain {FI\~q  
%     solution of Maxwell's curl equations over a three-dimensional w18RA#Zo/  
%     Cartesian space lattice comprised of uniform cubic grid cells. |7^^*UzSK:  
%     z\/53Sy<  
%     The dimensions of the computational domain are 8.2 cm R1-k3;v^  
%     (x-direction), 3.4 cm (y-direction), and 3.2 cm (z-direction).  X[W]=yJJ  
%     The grid is terminated with UPML absorbing boundary conditions. Lz\UZeq  
% K"2|[5  
%     An electric current source comprised of two collinear Jz components -K[782Q  
%     (realizing a Hertzian dipole) excites a radially propagating wave.  ?_`0G/xl  
%     The current source is located in the center of the grid.  The !?[oIQ)h  
%     source waveform is a differentiated Gaussian pulse given by ~B$b)`*  
%          J(t)=J0*(t-t0)*exp(-(t-t0)^2/tau^2), fB+b}aoV  
%     where tau=50 ps.  The FWHM spectral bandwidth of this zero-dc- !eJCM`cp  
%     content pulse is approximately 7 GHz. The grid resolution 7);:ZpDv%L  
%     (dx = 2 mm) was chosen to provide at least 10 samples per D^Ys)- d  
%     wavelength up through 15 GHz. lr2rQo>  
% f'6|OsVQ  
%     To execute this M-file, type "fdtd3D_UPML" at the MATLAB prompt.  s^T+5E&}  
% y)F!c29  
%     This code has been tested in the following Matlab environments: t?^9HP1b_  
%     Matlab version 6.1.0.450 Release 12.1 (May 18, 2001) Fpt-V  
%     Matlab version 6.5.1.199709 Release 13 Service Pack 1 (August 4, 2003) 9B+wYJp  
%     Matlab version 7.0.0.19920 R14 (May 6, 2004) uvA(Rn  
%     Matlab version 7.0.1.24704 R14 Service Pack 1 (September 13, 2004) =raA?Bp3;(  
%     Matlab version 7.0.4.365 R14 Service Pack 2 (January 29, 2005) :ZxLJK9x1  
% oY`qInM_  
%     Note: if you are using Matlab version 6.x, you may wish to make A.Bk/N1G  
%     one or more of the following modifications to this code: Bfd-:`Jk  
%       --uncomment line numbers 485 and 486 }xlKonk  
%       --comment out line numbers 552 and 561 hFvi5I-b  
% xiyxrR;  
%*********************************************************************** r?
l;I3~  
clear H0*5_OJ!i  
%*********************************************************************** ~@8+hnE]  
%     Fundamental constants o7+>G~i  
%*********************************************************************** K<v:-TjQZ:  
cc=2.99792458e8; 1'G8o=~  
muz=4.0*pi*1.0e-7; P.sg
RsL  
epsz=1.0/(cc*cc*muz); J`#`fX  
etaz=sqrt(muz/epsz); x:t<ZG&Xwg  
%*********************************************************************** s T3p
>8n  
%     Material parameters 0W>9'Rw  
%*********************************************************************** NfE.N&vI_c  
mur=1.0; &2EBk=X  
epsr=1.0; %M
cO6.M@  
eta=etaz*sqrt(mur/epsr); Pj7gGf6v  
%*********************************************************************** B(l-}|m_  
%     Grid parameters 5eX59:vtl  
% cbKL$|  
%     Each grid size variable name describes the number of sampled points !5x Ly6=}  
%     for a particular field component in the direction of that component. ["3d
f>!f  
%     Additionally, the variable names indicate the region of the grid Poa?Ej  
%     for which the dimension is relevant.  For example, ie_tot is the Qrz4}0  
%     number of sample points of Ex along the x-axis in the total J -Qh/d%]  
%     computational grid, and jh_bc is the number of sample points of Hy )'q%2%Ak  
%     along the y-axis in the y-normal UPML regions. xfSG~csoz  
% !k~z5z'=py  
%*********************************************************************** gY`Nr!O  
ie=60;          % Size of main grid zL s^,x  
je=60; !;>(ie\  
ke=60; a%hGZCI  
ih=ie+1; @XOi62(  
jh=je+1;   ]Qi,j#X  
kh=ke+1;   >kdM:MK  
upml=10;        % Thickness of PML boundaries wVi%
oSfM  
ih_bc=upml+1; s0{ NsK>  
jh_bc=upml+1;
z^sST  
kh_bc=upml+1; 9u wL{P&  
ie_tot=ie+2*upml;          % Size of total computational domain v5T9Y-{`  
je_tot=je+2*upml;        oVZ4bRl   
ke_tot=ke+2*upml;        .#^0pv!  
ih_tot=ie_tot+1; b"8FlZ$  
jh_tot=je_tot+1;          1a9w(
X  
kh_tot=ke_tot+1;          gZ(O)uzv  
is=round(ih_tot/2);         % Location of z-directed current source -Gsl[Rc0H;  
js=round(jh_tot/2); Nm8w/Q5D`  
ks=round(ke_tot/2); Q:S\0cI0  
%*********************************************************************** h'^FrWaU/  
%     Fundamental grid parameters 015Owi  
%*********************************************************************** @~jxG%y86  
delta=0.002;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ^i<}]c_|f  
dt=delta*sqrt(epsr*mur)/(2.0*cc); =LZ>su  
nmax=100; 7tH]*T9e>  
%*********************************************************************** J(d2:V{h  
%     Differentiated Gaussian pulse excitation I{dl%z73  
%*********************************************************************** j$eCe<.3  
rtau=50.0e-12; 2 X<nn  
tau=rtau/dt; L2Ux9_S  
ndelay=3*tau; W._G0b4}  
J0=1.0*epsz; ~HP LV  
%*********************************************************************** K="I<bK  
%     Initialize field arrays ]Sta]}VQ  
%*********************************************************************** w)S;J,Hv  
ex=zeros(ie_tot,jh_tot,kh_tot); _1JmjIH)M  
ey=zeros(ih_tot,je_tot,kh_tot); I
^fPk  
ez=zeros(ih_tot,jh_tot,ke_tot); ~4s-S3YzaM  
dx=zeros(ie_tot,jh_tot,kh_tot); oA3d^%(c  
dy=zeros(ih_tot,je_tot,kh_tot); '\*A"8;h  
dz=zeros(ih_tot,jh_tot,ke_tot); GhnE>d;i  
hx=zeros(ih_tot,je_tot,ke_tot); =|y|P80w  
hy=zeros(ie_tot,jh_tot,ke_tot); K[?R[  
hz=zeros(ie_tot,je_tot,kh_tot); fkW(Dt,  
bx=zeros(ih_tot,je_tot,ke_tot); xQz
#i-v  
by=zeros(ie_tot,jh_tot,ke_tot); 0&`}EXe<f  
bz=zeros(ie_tot,je_tot,kh_tot); f]}}yBte`  
%*********************************************************************** }(=ml7)v  
%     Initialize update coefficient arrays X}apxSd"  
%*********************************************************************** umDtp\  
C1ex=zeros(size(ex)); 7b.U!Ju  
C2ex=zeros(size(ex)); 0|2%# E  
C3ex=zeros(size(ex)); ez0\bym  
C4ex=zeros(size(ex)); |t\KsW  
C5ex=zeros(size(ex)); v7@
H\x*  
C6ex=zeros(size(ex)); Pi%tsKk%  
C1ey=zeros(size(ey)); \ j]~>9  
C2ey=zeros(size(ey)); ida*]+ ~  
C3ey=zeros(size(ey)); WzG07 2w  
C4ey=zeros(size(ey)); 8Nvr93T,  
C5ey=zeros(size(ey)); E:Y:X~vy  
C6ey=zeros(size(ey)); Tcs3>lJ}  
C1ez=zeros(size(ez)); *IlQ5+3I  
C2ez=zeros(size(ez)); -fM1$/]  
C3ez=zeros(size(ez)); \r]('x3S  
C4ez=zeros(size(ez)); q<}PM  
C5ez=zeros(size(ez)); hZ#\t  
C6ez=zeros(size(ez)); -]&<Sr-  
D1hx=zeros(size(hx)); 0=m&^Jpp  
D2hx=zeros(size(hx)); fI[dhd6  
D3hx=zeros(size(hx)); Z90Fcp:R  
D4hx=zeros(size(hx)); eH=c|m]!P  
D5hx=zeros(size(hx)); -q(:%;  
D6hx=zeros(size(hx)); L;C|ow^c  
D1hy=zeros(size(hy)); tG 7+7Z=  
D2hy=zeros(size(hy)); zZYH
c?Z  
D3hy=zeros(size(hy)); 5TET<f6R  
D4hy=zeros(size(hy)); >H?uuzi  
D5hy=zeros(size(hy)); 2$ VTu+  
D6hy=zeros(size(hy)); Wy)('EM  
D1hz=zeros(size(hz)); ~}%&p& p  
D2hz=zeros(size(hz)); Fhi5LhWe+.  
D3hz=zeros(size(hz)); `Y\QUj  
D4hz=zeros(size(hz)); 1OPfRDn.bk  
D5hz=zeros(size(hz)); +GPd  
D6hz=zeros(size(hz)); #f9qlM32  
%*********************************************************************** t|".=3%G  
%     Update coefficients, as described in Section 7.8.2. !>|`ly$6  
% cX"G7Bh  
%     In order to simplify the update equations used in the time-stepping vUfO4yfdg  
%     loop, we implement UPML update equations throughout the entire F=5kF/}x-z  
%     grid.  In the main grid, the electric-field update coefficients of 'WnpwY  
%     Equations 7.91a-f and the correponding magnetic field update %jgg59  
%     coefficients extracted from Equations 7.89 and 7.90 are simplified fSC.+,qk  
%     for the main grid (free space) and calculated below. gDc]^K4>  
% %9YA^ri  
%*********************************************************************** (lWKy9eTy`  
C1=1.0; &`sR){R  
C2=dt; {9:
hg9;E*  
C3=1.0; /^QFqM;  
C4=1.0/2.0/epsr/epsr/epsz/epsz; HtS#_y%(  
C5=2.0*epsr*epsz; /~"AG l.  
C6=2.0*epsr*epsz; '7=<#Blc  
D1=1.0; j~$)c)h"  
D2=dt; 2E([#Pzb  
D3=1.0; T@>63  
D4=1.0/2.0/epsr/epsz/mur/muz; 5YLho2h38!  
D5=2.0*epsr*epsz; DUPmq!A  
D6=2.0*epsr*epsz; 0, "ZV}  
%*********************************************************************** _6/Qp`s  
%     Initialize main grid update coefficients ^X(_zinN"  
%***********************************************************************
~s0P FS7  
C1ex(:,jh_bc:jh_tot-upml,:)=C1;     mlmnkgl ]  
C2ex(:,jh_bc:jh_tot-upml,:)=C2; 0Sz/c+ 6  
C3ex(:,:,kh_bc:kh_tot-upml)=C3; 1[# =,  
C4ex(:,:,kh_bc:kh_tot-upml)=C4; DMRs}Yz6  
C5ex(ih_bc:ie_tot-upml,:,:)=C5; ?xf~!D  
C6ex(ih_bc:ie_tot-upml,:,:)=C6; 9Iy[E,j  
C1ey(:,:,kh_bc:kh_tot-upml)=C1; wtUG^hV #_  
C2ey(:,:,kh_bc:kh_tot-upml)=C2; 3V!W@[ }:  
C3ey(ih_bc:ih_tot-upml,:,:)=C3; ^zkd{ov  
C4ey(ih_bc:ih_tot-upml,:,:)=C4; B4 <_"0  
C5ey(:,jh_bc:je_tot-upml,:)=C5; *ig5Q(b*N  
C6ey(:,jh_bc:je_tot-upml,:)=C6; I$F\(]"@  
C1ez(ih_bc:ih_tot-upml,:,:)=C1; n_ OUWvs  
C2ez(ih_bc:ih_tot-upml,:,:)=C2; s!=!A  
C3ez(:,jh_bc:jh_tot-upml,:)=C3; 0YHYxn  
C4ez(:,jh_bc:jh_tot-upml,:)=C4; I>8Bc  
C5ez(:,:,kh_bc:ke_tot-upml)=C5; {vCU^BN,k  
C6ez(:,:,kh_bc:ke_tot-upml)=C6; H|'$dO)W  
D1hx(:,jh_bc:je_tot-upml,:)=D1; "9.6\Y\*  
D2hx(:,jh_bc:je_tot-upml,:)=D2; $&NbLjeS  
D3hx(:,:,kh_bc:ke_tot-upml)=D3; <|wmjW/D  
D4hx(:,:,kh_bc:ke_tot-upml)=D4;  aeQ{_SK  
D5hx(ih_bc:ih_tot-upml,:,:)=D5; ?~]>H A:  
D6hx(ih_bc:ih_tot-upml,:,:)=D6; VN!^m]0  
D1hy(:,:,kh_bc:ke_tot-upml)=D1; H.f9d.<W%  
D2hy(:,:,kh_bc:ke_tot-upml)=D2; r$Kh3EEF`E  
D3hy(ih_bc:ie_tot-upml,:,:)=D3; s>V*=#L  
D4hy(ih_bc:ie_tot-upml,:,:)=D4; G+;g:_E=  
D5hy(:,jh_bc:jh_tot-upml,:)=D5; Ubp
g92  
D6hy(:,jh_bc:jh_tot-upml,:)=D6; OHQ3+WJ  
D1hz(ih_bc:ie_tot-upml,:,:)=D1; MQhYJ01i  
D2hz(ih_bc:ie_tot-upml,:,:)=D2; F}\[eFf[  
D3hz(:,jh_bc:je_tot-upml,:)=D3; X^9t  
D4hz(:,jh_bc:je_tot-upml,:)=D4; EywZIw?mjX  
D5hz(:,:,kh_bc:kh_tot-upml)=D5; MEDskvBG  
D6hz(:,:,kh_bc:kh_tot-upml)=D6; Psg +\14  
%*********************************************************************** !fif8kf  
%     Fill in PML regions X{4xm,B/  
% PwRNBb}6  
%     PML theory describes a continuous grading of the material properties c '/2F0y  
%     over the PML region.  In the FDTD grid it is necessary to discretize 7?*~oV
ZW  
%     the grading by averaging the material properties over a grid cell F'fM?!(  
%     width centered on each field component.  As an example of the pq@$&G  
%     implementation of this averaging, we take the integral of the '$lw[1  
%     continuous sigma(x) in the PML region HL|0d }  
%   "F04c|oR<X  
%         sigma_i = integral(sigma(x))/delta &<m WA]cAL  
%   C[^a/P`i  
%     where the integral is over a single grid cell width in x, and is WodF -bE  
%     bounded by x1 and x2.  Applying this to the polynomial grading of sCuQBZ h  
%     Equation 7.60a produces ]eORw$f  
%  ; \Y-  
%         sigma_i = (x2^(m+1)-x1^(m+1))*sigmam/(delta*(m+1)*d^m) s;~J2h[  
% aj% `x4eA  
%     where sigmam is the maximum value of sigma as described by Equation if\k[O 1T6  
%     7.62. 1Q3%!~<\s  
%         (d1V1t2r6  
%     The definitions of x1 and x2 depend on the position of the field F&-5&'6G+  
%     component within the grid cell.  We have either ok8JnQC
 
% ;o&_:]S  
%         x1 = (i-0.5)*delta,  x2 = (i+0.5)*delta 2I9{
+>k  
%  Ebmqq#SHjX  
%     or %pQdq[J={  
%  REeD?u j  
%         x1 = (i)*delta,      x2 = (i+1)*delta 7ZrJ#n8?ih  
% a^(S!I  
%     where i varies over the PML region. 8'>.#vyMGv  
%  0R.Gjz*Q  
%*********************************************************************** G>9'5Lt  
rmax=exp(-16);  %desired reflection error, designated as R(0) in Equation 7.62 ,|=iv  
orderbc=4;      %order of the polynomial grading, designated as m in Equation 7.60a,b "Vs Nyy  
%   x-varying material properties tVK?VNW  
delbc=upml*delta; `=WzG"  
sigmam=-log(rmax)*(orderbc+1.0)/(2.0*eta*delbc); <PH3gyC  
sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1.0)); qM+!f2t  
kmax=1; %@)U/G6s}  
kfactor=(kmax-1.0)/delta/(orderbc+1.0)/delbc^orderbc; ~0:$G?fz  
for i=1:upml
5LnB]dW  
    ?d,acm  
    % Coefficients for field components in the center of the grid cell a .B\=3xn  
    x1=(upml-i+1)*delta; JT}dor  
    x2=(upml-i)*delta; jCdZ}M($  
    sigma=sigfactor*(x1^(orderbc+1)-x2^(orderbc+1)); !:`Ra  
    ki=1+kfactor*(x1^(orderbc+1)-x2^(orderbc+1)); -uS7~Ww.a  
    facm=(2*epsr*epsz*ki-sigma*dt); j>0SE  
    facp=(2*epsr*epsz*ki+sigma*dt); B]CS2LEqh  
    C5ex(i,:,:)=facp; \L&qfMjW"Z  
    C5ex(ie_tot-i+1,:,:)=facp; 7~QwlU3n<F  
    C6ex(i,:,:)=facm; $Ykp8u,(  
    C6ex(ie_tot-i+1,:,:)=facm; V1AEjh  
    D1hz(i,:,:)=facm/facp; #<^/yoH7C6  
    D1hz(ie_tot-i+1,:,:)=facm/facp; .T[!!z#^  
    D2hz(i,:,:)=2.0*epsr*epsz*dt/facp; u&Ie%@:h9R  
    D2hz(ie_tot-i+1,:,:)=2.0*epsr*epsz*dt/facp; E}E7VQjM  
    D3hy(i,:,:)=facm/facp; \:18Uoe7  
    D3hy(ie_tot-i+1,:,:)=facm/facp; t2`X!`  
    D4hy(i,:,:)=1.0/facp/mur/muz; EC:x,i  
    D4hy(ie_tot-i+1,:,:)=1.0/facp/mur/muz; .'+|>6eU  
    % Coefficients for field components on the grid cell boundary [7CH(o1a&  
    x1=(upml-i+1.5)*delta; :zS>^RE  
    x2=(upml-i+0.5)*delta; gb^UFD L  
    sigma=sigfactor*(x1^(orderbc+1)-x2^(orderbc+1)); Qt)7m
f  
    ki=1.0+kfactor*(x1^(orderbc+1)-x2^(orderbc+1)); CgaB)`.  
    facm=(2.0*epsr*epsz*ki-sigma*dt); oc,U4+T  
    facp=(2.0*epsr*epsz*ki+sigma*dt); c>%z)uY>/  
    C1ez(i,:,:)=facm/facp; $/-wgyP3m+  
    C1ez(ih_tot-i+1,:,:)=facm/facp; +$Ddd`J'  
    C2ez(i,:,:)=2.0*epsr*epsz*dt/facp; &St~!y6M?  
    C2ez(ih_tot-i+1,:,:)=2.0*epsr*epsz*dt/facp; gPn%`_d5  
    C3ey(i,:,:)=facm/facp; 'ocwXyP,  
    C3ey(ih_tot-i+1,:,:)=facm/facp; gw*yIZ@3)  
    C4ey(i,:,:)=1.0/facp/epsr/epsz; 'R-JQE-]  
    C4ey(ih_tot-i+1,:,:)=1.0/facp/epsr/epsz; cPa 0n4  
    D5hx(i,:,:)=facp; !~VR|n-  
    D5hx(ih_tot-i+1,:,:)=facp; "fX8xZdS  
    D6hx(i,:,:)=facm; QB oZCLv  
    D6hx(ih_tot-i+1,:,:)=facm; +;4AG::GN  
    j/; @P  
end ~gcst;  
%   PEC walls co(fGp#!  
C1ez(1,:,:)=-1.0; P,+0   
C1ez(ih_tot,:,:)=-1.0; 8M5!5Jzv  
C2ez(1,:,:)=0.0; UeC%Wa<[  
C2ez(ih_tot,:,:)=0.0; {jCu9 ]c!  
C3ey(1,:,:)=-1.0; rA&#>R`  
C3ey(ih_tot,:,:)=-1.0; WL*W=(  
C4ey(1,:,:)=0.0; *U5>j#,  
C4ey(ih_tot,:,:)=0.0; G> 5=`  
%   y-varying material properties tleK(^  
delbc=upml*delta; P"[l86:  
sigmam=-log(rmax)*epsr*epsz*cc*(orderbc+1.0)/(2.0*delbc); Z{|.xgsY  
sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1.0)); 2Q;Y@%G  
kmax=1.0; 6H. L!tUI  
kfactor=(kmax-1.0)/delta/(orderbc+1.0)/delbc^orderbc; h)aWerzL  
for j=1:upml lFzQG:k@  
    aL$c).hq0  
    % Coefficients for field components in the center of the grid cell cnjj) c  
    y1=(upml-j+1)*delta; EIVQu~,H  
    y2=(upml-j)*delta; g~WNL^GGS  
    sigma=sigfactor*(y1^(orderbc+1)-y2^(orderbc+1)); 9x 6ca  
    ki=1+kfactor*(y1^(orderbc+1)-y2^(orderbc+1)); hBs>2u|z9  
    facm=(2*epsr*epsz*ki-sigma*dt); tkUW)ScJ  
    facp=(2*epsr*epsz*ki+sigma*dt); UO7a}Tz<  
    2TevdyI  
    C5ey(:,j,:)=facp; GurE7J^=  
    C5ey(:,je_tot-j+1,:)=facp; {?kKpMNNn  
    C6ey(:,j,:)=facm; `)xU;-  
    C6ey(:,je_tot-j+1,:)=facm; ^FyvaO  
    D1hx(:,j,:)=facm/facp; +{ ,w#@  
    D1hx(:,je_tot-j+1,:)=facm/facp; <o"D/<XnB3  
    D2hx(:,j,:)=2*epsr*epsz*dt/facp; _tR.RAaa"  
    D2hx(:,je_tot-j+1,:)=2*epsr*epsz*dt/facp; vYTPZ@RL  
    D3hz(:,j,:)=facm/facp; bx%hizb  
    D3hz(:,je_tot-j+1,:)=facm/facp; ct|'I]nB.h  
    D4hz(:,j,:)=1/facp/mur/muz; { <ao4w6B  
    D4hz(:,je_tot-j+1,:)=1/facp/mur/muz; 9_CA5?y$:  
    %V" +}Dr  
    % Coefficients for field components on the grid cell boundary 9lazo  
    y1=(upml-j+1.5)*delta; Wy0a2Ve  
    y2=(upml-j+0.5)*delta; N!4xP.Ps  
    sigma=sigfactor*(y1^(orderbc+1)-y2^(orderbc+1)); uk>/Il  
    ki=1+kfactor*(y1^(orderbc+1)-y2^(orderbc+1)); emDvy2uA#  
    facm=(2*epsr*epsz*ki-sigma*dt); XiZ Zo  
    facp=(2*epsr*epsz*ki+sigma*dt);    }Rf:DmPE  
     \c@qtIc  
    C1ex(:,j,:)=facm/facp; Qe=Q8c
T  
    C1ex(:,jh_tot-j+1,:)=facm/facp; $dWl A<u  
    C2ex(:,j,:)=2*epsr*epsz*dt/facp; th>yi)m  
    C2ex(:,jh_tot-j+1,:)=2*epsr*epsz*dt/facp; [laL6  
    C3ez(:,j,:)=facm/facp; N_WA4?rB  
    C3ez(:,jh_tot-j+1,:)=facm/facp; 7;#dX~>@{  
    C4ez(:,j,:)=1/facp/epsr/epsz; xF:poi  
    C4ez(:,jh_tot-j+1,:)=1/facp/epsr/epsz;   1rzq$,O  
    D5hy(:,j,:)=facp; C">=2OO  
    D5hy(:,jh_tot-j+1,:)=facp; h-v&I>  
    D6hy(:,j,:)=facm; h0EGhJs  
    D6hy(:,jh_tot-j+1,:)=facm; )Q pP1
[  
end ]mGsNQ ].H  
%   PEC walls =B 4gEWR  
C1ex(:,1,:)=-1; ,uD*FSp>  
C1ex(:,jh_tot,:)=-1; a[{QlD^D  
C2ex(:,1,:)=0; L)//- k9  
C2ex(:,jh_tot,:)=0; }C2i#;b  
C3ez(:,1,:)=-1; QW"6]  
C3ez(:,jh_tot,:)=-1; @_weMz8}  
C4ez(:,1,:)=0; :Pp;{=J  
C4ez(:,jh_tot,:)=0;   \~1zAiSd>#  
%   z-varying material properties E'kQ  
delbc=upml*delta; m3v*,~  
sigmam=-log(rmax)*epsr*epsz*cc*(orderbc+1)/(2*delbc); 3YNkT"~T  
sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1)); VRN9yn2  
kmax=1; *R~(:z>>  
kfactor=(kmax-1)/delta/(orderbc+1)/delbc^orderbc; F|ML$  
for k=1:upml x /Ky: Ky  
    % Coefficients for field components in the center of the grid cell p2_Zsq  
    z1=(upml-k+1)*delta; uB)6\fkTB  
    z2=(upml-k)*delta; e
z<wEtS  
    sigma=sigfactor*(z1^(orderbc+1)-z2^(orderbc+1)); MBTt'6M  
    ki=1+kfactor*(z1^(orderbc+1)-z2^(orderbc+1)); <t|9`l_XW  
    facm=(2*epsr*epsz*ki-sigma*dt); 2g07wJ6x  
    facp=(2*epsr*epsz*ki+sigma*dt); !`S%l1[Z  
    Q)|LiCR,  
    C5ez(:,:,k)=facp; keq[6Lv  
    C5ez(:,:,ke_tot-k+1)=facp; 0 VgnN
  
    C6ez(:,:,k)=facm; pQ2)M8 gf  
    C6ez(:,:,ke_tot-k+1)=facm; QCY{D@7T  
    D1hy(:,:,k)=facm/facp; Pc-8L]2oaF  
    D1hy(:,:,ke_tot-k+1)=facm/facp; TH}ycue  
    D2hy(:,:,k)=2*epsr*epsz*dt/facp; Q24:G  
    D2hy(:,:,ke_tot-k+1)=2*epsr*epsz*dt/facp; (Vv[  
    D3hx(:,:,k)=facm/facp; 01d26`G$i~  
    D3hx(:,:,ke_tot-k+1)=facm/facp; {&j{V-}f  
    D4hx(:,:,k)=1/facp/mur/muz; $T%<'=u|E  
    D4hx(:,:,ke_tot-k+1)=1/facp/mur/muz; o>,z %+  
    .j4ziRa-  
    % Coefficients for field components on the grid cell boundary r.G/f{=<@  
    z1=(upml-k+1.5)*delta; ;s/b_RN  
    z2=(upml-k+0.5)*delta; })u}PQ  
    sigma=sigfactor*(z1^(orderbc+1)-z2^(orderbc+1)); R*5;J`TW  
    ki=1+kfactor*(z1^(orderbc+1)-z2^(orderbc+1)); %D6Wlf+^n  
    facm=(2*epsr*epsz*ki-sigma*dt); JFk|Uqs(  
    facp=(2*epsr*epsz*ki+sigma*dt); Vv"wf;#  
    ).`a-Pv  
    C1ey(:,:,k)=facm/facp; )p_LkX(  
    C1ey(:,:,kh_tot-k+1)=facm/facp; s&_O2
(l  
    C2ey(:,:,k)=2*epsr*epsz*dt/facp; X}QmeY[0I  
    C2ey(:,:,kh_tot-k+1)=2*epsr*epsz*dt/facp; m_U6"\n 5  
    C3ex(:,:,k)=facm/facp; tMk>Bx9[  
    C3ex(:,:,kh_tot-k+1)=facm/facp; <TROs!x$a  
    C4ex(:,:,k)=1/facp/epsr/epsz; u9^;~i,  
    C4ex(:,:,kh_tot-k+1)=1/facp/epsr/epsz; HG6{`
i  
    D5hz(:,:,k)=facp; (uxQBy  
    D5hz(:,:,kh_tot-k+1)=facp; 6V?&hq&t  
    D6hz(:,:,k)=facm; GB0b|9(6D"  
    D6hz(:,:,kh_tot-k+1)=facm; Ch!Q?4  
end 21EUP6}8j  
%   PEC walls TF3q?0  
C1ey(:,:,1)=-1; x$b[m20  
C1ey(:,:,kh_tot)=-1; <T% hfW
 
C2ey(:,:,1)=0; J?ZVzKTb>}  
C2ey(:,:,kh_tot)=0; GSb)|mj  
C3ex(:,:,1)=-1; ,* vnt6C*  
C3ex(:,:,kh_tot)=-1; Tb6x@MorP  
C4ex(:,:,1)=0; I;4Cv
oT  
C4ex(:,:,kh_tot)=0; (tz]!Aa{s  
%figure !?M_%fNE  
%set(gcf,'DoubleBuffer','on') {3uSg)  
%*********************************************************************** umrI4.1c  
%     Begin time stepping loop HNyDWD)_  
%*********************************************************************** A;nmua-Fv  
d0=zeros(1,nmax);d1=zeros(1,nmax);d2=zeros(1,nmax); -
&$%m)wN  
for n=1:nmax m%(JRh  
     bQQ/7KM  
    % Update magnetic field !eoec2h#5  
    bstore=bx; B148wh#r  
    bx(2:ie_tot,:,:)=D1hx(2:ie_tot,:,:).*  bx(2:ie_tot,:,:)-... `.-k%2?/  
                     D2hx(2:ie_tot,:,:).*((ez(2:ie_tot,2:jh_tot,:)-ez(2:ie_tot,1:je_tot,:))-... ;= @-j@?  
                                          (ey(2:ie_tot,:,2:kh_tot)-ey(2:ie_tot,:,1:ke_tot)))./delta; snNg:rTL  
    hx(2:ie_tot,:,:)= D3hx(2:ie_tot,:,:).*hx(2:ie_tot,:,:)+... 5BO!K$6  
                      D4hx(2:ie_tot,:,:).*(D5hx(2:ie_tot,:,:).*bx(2:ie_tot,:,:)-... 10a*7 L  
                                           D6hx(2:ie_tot,:,:).*bstore(2:ie_tot,:,:)); y#0Z[[I0  
    bstore=by; F"P:9`/  
    by(:,2:je_tot,:)=D1hy(:,2:je_tot,:).*  by(:,2:je_tot,:)-... +VCo=oA  
                     D2hy(:,2:je_tot,:).*((ex(:,2:je_tot,2:kh_tot)-ex(:,2:je_tot,1:ke_tot))-... CN>};>WlG  
                                          (ez(2:ih_tot,2:je_tot,:)-ez(1:ie_tot,2:je_tot,:)))./delta; pXlBKJmW  
    hy(:,2:je_tot,:)= D3hy(:,2:je_tot,:).*hy(:,2:je_tot,:)+... !!#ale&  
                      D4hy(:,2:je_tot,:).*(D5hy(:,2:je_tot,:).*by(:,2:je_tot,:)-... 1mPS)X_  
                                           D6hy(:,2:je_tot,:).*bstore(:,2:je_tot,:)); 3Qr!?=nf  
    bstore=bz; )Nd:PnA  
    bz(:,:,2:ke_tot)=D1hz(:,:,2:ke_tot).*  bz(:,:,2:ke_tot)-... M!DoR6  
                     D2hz(:,:,2:ke_tot).*((ey(2:ih_tot,:,2:ke_tot)-ey(1:ie_tot,:,2:ke_tot))-... eQIi}\`  
                                          (ex(:,2:jh_tot,2:ke_tot)-ex(:,1:je_tot,2:ke_tot)))./delta; *B}R4Y|g  
    hz(:,:,2:ke_tot)= D3hz(:,:,2:ke_tot).*hz(:,:,2:ke_tot)+... <RsKV$Je I  
                      D4hz(:,:,2:ke_tot).*(D5hz(:,:,2:ke_tot).*bz(:,:,2:ke_tot)-... ojqX#>0K  
                                           D6hz(:,:,2:ke_tot).*bstore(:,:,2:ke_tot)); >-+X;0&  
    `B0*/ml  
    % Update electric field =w?cp}HW  
    dstore=dx; IL"#TKKv  
    dx(:,2:je_tot,2:ke_tot)=C1ex(:,2:je_tot,2:ke_tot).*  dx(:,2:je_tot,2:ke_tot)+... Elk$9 <<  
                            C2ex(:,2:je_tot,2:ke_tot).*((hz(:,2:je_tot,2:ke_tot)-hz(:,1:je_tot-1,2:ke_tot))-... Sf:lN4  
                                                        (hy(:,2:je_tot,2:ke_tot)-hy(:,2:je_tot,1:ke_tot-1)))./delta; zate%y  
    ex(:,2:je_tot,2:ke_tot)=C3ex(:,2:je_tot,2:ke_tot).*ex(:,2:je_tot,2:ke_tot)+... htSk2N/  
                            C4ex(:,2:je_tot,2:ke_tot).*(C5ex(:,2:je_tot,2:ke_tot).*dx(:,2:je_tot,2:ke_tot)-... rAdcMFW  
                                                        C6ex(:,2:je_tot,2:ke_tot).*dstore(:,2:je_tot,2:ke_tot)); Ym2
m1  
    dstore=dy; ;mxT>|z  
    dy(2:ie_tot,:,2:ke_tot)=C1ey(2:ie_tot,:,2:ke_tot).*  dy(2:ie_tot,:,2:ke_tot)+... SLB iQd.  
                            C2ey(2:ie_tot,:,2:ke_tot).*((hx(2:ie_tot,:,2:ke_tot)-hx(2:ie_tot,:,1:ke_tot-1))-... 6oBt<r?CJ  
                                                        (hz(2:ie_tot,:,2:ke_tot)-hz(1:ie_tot-1,:,2:ke_tot)))./delta; VQ3&  
    ey(2:ie_tot,:,2:ke_tot)=C3ey(2:ie_tot,:,2:ke_tot).*ey(2:ie_tot,:,2:ke_tot)+... W>IKy#  
                            C4ey(2:ie_tot,:,2:ke_tot).*(C5ey(2:ie_tot,:,2:ke_tot).*dy(2:ie_tot,:,2:ke_tot)-... qr;" K?NX  
                                                        C6ey(2:ie_tot,:,2:ke_tot).*dstore(2:ie_tot,:,2:ke_tot)); XdVC>6  
    dstore=dz; `8bp6}OD,  
    dz(2:ie_tot,2:je_tot,:)=C1ez(2:ie_tot,2:je_tot,:).*  dz(2:ie_tot,2:je_tot,:)+... $iJ #%&D  
                            C2ez(2:ie_tot,2:je_to .. 1T?%i  
LMzYsXG*[  

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

哎，连个回答俺的也没有 NFPkK?+  

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

高手们都说说 B Xp3u|t  

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

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

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

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

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