[求助]请教个UPML程序中的问题

我在看该3D的UPML程序时，有个地方看不懂，不知道是从哪个公式写出来的，请高人指点，非常谢谢了~ v2B0q4*BS?  
问题是： 9~p[  
1.sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1.0));这句的公式是对应Allen书上的（7.60a）式吗，它是怎么对应X的M次方的？ L<ue$'  
2.C1ex(:,jh_bc:jh_tot-upml,:)=C1;等一系列的赋值的范围如jh_bc:jh_tot-upml是根据什么定的啊？ Dp)=0<$y  
!HnXXVW  
该程序如下： 7R6ry(6N  
%*********************************************************************** A-ZN F4  
%     3-D FDTD code with UPML absorbing boundary conditions Dpl A?  
%*********************************************************************** /ro=?QYb  
% U56G.  
%     Program author: Keely J. Willis, Graduate Student S/9DtXQ  
%                     UW Computational Electromagnetics Laboratory +VO-oFE|  
%                           Director: Susan C. Hagness yXHUJgjl/  
%                     Department of Electrical and Computer Engineering 2/"u
5  
%                     University of Wisconsin-Madison (NF~Ck$#q  
%                     1415 Engineering Drive G+X Sfr  
%                     Madison, WI 53706-1691 ";3zXk[#  
%                     kjwillis@wisc.edu  %-c*C$  
% hw= F
t4L  
%     Copyright 2005 I'uSp-Sfy  
% 2E}*v5b,  
%     This MATLAB M-file implements the finite-difference time-domain VXR>]HUF  
%     solution of Maxwell's curl equations over a three-dimensional ;[M}MFc/`  
%     Cartesian (笛卡尔)space lattice comprised of uniform cubic grid cells. BT}!W`  
%     hRUhX[  
%     The dimensions of the computational domain are 8.2 cm #,":vr  
%     (x-direction), 3.4 cm (y-direction), and 3.2 cm (z-direction).   C~o7X^[R\  
%     The grid is terminated with UPML absorbing boundary conditions. W g02 A\  
% 7f0lQ  
%     An electric current source comprised of two collinea（共线的）r Jz components ;#vKi0V7  
%     (realizing a Hertzian dipole) excites a radially propagating wave.   vb<oi&X  
%     The current source is located in the center of the grid.  The BYVY)<v/  
%     source waveform is a differentiated Gaussian pulse given by uVJDne,R  
%          J(t)=J0*(t-t0)*exp(-(t-t0)^2/tau^2), FaDjLo2'o  
%     where tau=50 ps.  The FWHM spectral bandwidth of this zero-dc- ^eo|P~w g  
%     content pulse is approximately 7 GHz. The grid resolution J&.{7YF  
%     (dx = 2 mm) was chosen to provide at least 10 samples per 3|'>`!hb  
%     wavelength up through 15 GHz. C|}iCB  
% Hn5|B 3vN  
%     To execute this M-file, type "fdtd3D_UPML" at the MATLAB prompt.   #a'Ex=%rM  
% Gn ~6X-l  
%     This code has been tested in the following Matlab environments: G 8g<>d{j  
%     Matlab version 6.1.0.450 Release 12.1 (May 18, 2001) 'VzP};  
%     Matlab version 6.5.1.199709 Release 13 Service Pack 1 (August 4, 2003) kXi6lh  
%     Matlab version 7.0.0.19920 R14 (May 6, 2004) kBD>-5Sn_T  
%     Matlab version 7.0.1.24704 R14 Service Pack 1 (September 13, 2004) j4|N-:  
%     Matlab version 7.0.4.365 R14 Service Pack 2 (January 29, 2005) inh=WUEW  
% MP_ ~<Q  
%     Note: if you are using Matlab version 6.x, you may wish to make 05b_)&4R  
%     one or more of the following modifications to this code: <vONmE a  
%       --uncomment line numbers 485 and 486 6W]9$n\"?  
%       --comment out line numbers 552 and 561 +@p% p  
% i}.&0Fp  
%*********************************************************************** ;"dV"W  
i36eBjT  
clear i{`FmrPO~  
I%sFqh>  
%*********************************************************************** &#!4XOyB  
%     Fundamental constants +l9!Fl{MK\  
%*********************************************************************** pe]A5\4c  
:h\Q;?  
cc=2.99792458e8; :,'wVS8"]  
muz=4.0*pi*1.0e-7; E4|jOz^j4\  
epsz=1.0/(cc*cc*muz); nxWY7hU  
etaz=sqrt(muz/epsz); a@@)6FM  
9~]~#Uj  
%*********************************************************************** NQ(1   
%     Material parameters ]
n_ k`  
%*********************************************************************** psg)*'r  
BJM.iXU)[  
mur=1.0; k1{K*O$e  
epsr=1.0; Ju96#v+:  
eta=etaz*sqrt(mur/epsr); g#Sl
%Y  
/aZ+T5O  
%*********************************************************************** 1(I6.BHW  
%     Grid parameters Y }$/e  
% RS)tO0  
%     Each grid size variable name describes the number of sampled points a yCY~=i  
%     for a particular field component in the direction of that component. 9Iwe2lu  
%     Additionally, the variable names indicate the region of the grid pTPi@SBaP{  
%     for which the dimension is relevant.  For example, ie_tot is the Jk7|{W\OA  
%     number of sample points of Ex along the x-axis in the total bdC8zDD  
%     computational grid, and jh_bc is the number of sample points of Hy [5G6VNh=  
%     along the y-axis in the y-normal UPML regions. DW5Y@;[  
% m[~V/N3  
%*********************************************************************** y9q8i(E0  
S- pV_Ff  
ie=41;          % Size of main grid oSyyd  
je=17; /& Jan:  
ke=16; *h!28Ya(~  
ih=ie+1; U+:m4a  
jh=je+1;   r'^Hg/Jzt  
kh=ke+1;   pEBM3r!X  
o4m\~as)Y  
upml=10;        % Thickness of PML boundaries k5:G-BQ:  
ih_bc=upml+1; %?}33yV  
jh_bc=upml+1; OSs&r$  
kh_bc=upml+1; LhOa{1SY  
B@&4i?yJ  
ie_tot=ie+2*upml;          % Size of total computational domain _jLL_GD  
je_tot=je+2*upml;         WY?[,_4U  
ke_tot=ke+2*upml;         t&H?\)!4  
ih_tot=ie_tot+1; $c0h.t  
jh_tot=je_tot+1;           7gj4j^a^]{  
kh_tot=ke_tot+1;           }+.}J  
wQ^EYKD  
is=round(ih_tot/2);         % Location of z-directed current source \fG#7_wt  
js=round(jh_tot/2); ';3{T:I  
ks=round(ke_tot/2); "e.jZcN*  
cAY:AtD  
%*********************************************************************** (7*%K&x  
%     Fundamental grid parameters F5P[dp-`1  
%*********************************************************************** >,F bX8Zz  
B:'J`M"N  
delta=0.002; A#.edVj.g4  
dt=delta*sqrt(epsr*mur)/(2.0*cc); @;*Ksy@1O  
nmax=100; =/s>Q l  
h"X;3b^ m  
%*********************************************************************** ,?>s>bHV  
%     Differentiated Gaussian pulse excitation ,]9P{k]O  
%*********************************************************************** WI-&x '  
?[@J8  
rtau=50.0e-12; ZrPbl"`7  
tau=rtau/dt; K#6P}tf  
ndelay=3*tau; o /j*d3  
J0=-1.0*epsz; w>pq+og&  
6:b!F  
%***********************************************************************
jrr EAp  
%     Initialize field arrays w65K[l;2  
%*********************************************************************** F*IzQ(#HW  
)
J2mM  
ex=zeros(ie_tot,jh_tot,kh_tot); B2
>H_dmQ  
ey=zeros(ih_tot,je_tot,kh_tot); &]`(v}`]  
ez=zeros(ih_tot,jh_tot,ke_tot); \.MR""@y`{  
dx=zeros(ie_tot,jh_tot,kh_tot); ,:%CB"J  
dy=zeros(ih_tot,je_tot,kh_tot); "I3@m%qv  
dz=zeros(ih_tot,jh_tot,ke_tot); U{2BVqM  
OLyf8&AU@  
hx=zeros(ih_tot,je_tot,ke_tot); [@VP?74  
hy=zeros(ie_tot,jh_tot,ke_tot); ;_c;0)  
hz=zeros(ie_tot,je_tot,kh_tot); \k.{-nh  
bx=zeros(ih_tot,je_tot,ke_tot); yA)/Q Yge  
by=zeros(ie_tot,jh_tot,ke_tot); '(U-(wTC'/  
bz=zeros(ie_tot,je_tot,kh_tot);  _zY#U9  
>0/i[k-dk  
%*********************************************************************** ]{3)^axW;  
%     Initialize update coefficient arrays EMY/~bQW  
%*********************************************************************** 24 [+pu  
4ezEW|S  
C1ex=zeros(size(ex)); 4{y)TZ  
C2ex=zeros(size(ex)); q]T1dz?  
C3ex=zeros(size(ex)); 8IlunJ  
C4ex=zeros(size(ex)); VCV"S>aVf  
C5ex=zeros(size(ex)); bo2H]PL*  
C6ex=zeros(size(ex)); J1( 9QN[w  
A,JmX  
C1ey=zeros(size(ey)); cqQ#p2<%  
C2ey=zeros(size(ey)); u(@$a4z  
C3ey=zeros(size(ey)); kZR8a(4D  
C4ey=zeros(size(ey)); wxK
X{Bs  
C5ey=zeros(size(ey)); k.uH~S_  
C6ey=zeros(size(ey)); ??^5;P{yx  
~ksi</s  
C1ez=zeros(size(ez)); 6a[}'/  
C2ez=zeros(size(ez)); dq(uVW^&ae  
C3ez=zeros(size(ez)); Ro\8ZXUQa  
C4ez=zeros(size(ez)); ;&9)I8Us  
C5ez=zeros(size(ez)); 5eLtCsHz  
C6ez=zeros(size(ez)); ]D?"aX'q>  
'~5LY!H(pT  
D1hx=zeros(size(hx)); YWe{juXSw  
D2hx=zeros(size(hx)); m8A#~i .  
D3hx=zeros(size(hx)); KI)M JG:t  
D4hx=zeros(size(hx)); |,S+@"0#  
D5hx=zeros(size(hx)); &v56#lG  
D6hx=zeros(size(hx)); Lt u'W22  
gwLf'  
D1hy=zeros(size(hy)); }tRm]w  
D2hy=zeros(size(hy)); Nv#t:J9f  
D3hy=zeros(size(hy)); Bo)3!wO8  
D4hy=zeros(size(hy)); T$>WE= Y  
D5hy=zeros(size(hy)); sd*p/Q|4  
D6hy=zeros(size(hy)); 6</xL9#/  
v[ .cd*b  
D1hz=zeros(size(hz)); e%svrJ2   
D2hz=zeros(size(hz)); CSJdvxb  
D3hz=zeros(size(hz)); e^8 O_VB  
D4hz=zeros(size(hz)); {j9{n  
D5hz=zeros(size(hz)); joFm]3$;  
D6hz=zeros(size(hz)); RSfQNc9Z  

zwhe  
%*********************************************************************** &Y!-%{e  
%     Update coefficients, as described in Section 7.8.2. gqZ'$7So  
% I3aNFa}  
%     In order to simplify the update equations used in the time-stepping
D9<!mH  
%     loop, we implement UPML update equations throughout the entire .CL[_;}  
%     grid.  In the main grid, the electric-field update coefficients of Lf16j*}-Q  
%     Equations 7.91a-f and the correponding magnetic field update Pr3qo4t.L  
%     coefficients extracted from Equations 7.89 and 7.90 are simplified BBaQ}{F8>2  
%     for the main grid (free space) and calculated below. = I:.X ;  
% .!Oo|m`V@  
%*********************************************************************** FtpK)9/4  
P@Hs`=  
C1=1.0; Z?6%;n^ 54  
C2=dt; DVG(Vw  
C3=1.0; 5&QJ7B,!  
C4=1.0/2.0/epsr/epsr/epsz/epsz; {:Orn%Q  
C5=2.0*epsr*epsz; |t^E~HLm,  
C6=2.0*epsr*epsz; p, h9D_  
FEW14U'O  
D1=1.0; $|L Sx  
D2=dt; @Pm>sY}d<I  
D3=1.0; VrHv)lUr  
D4=1.0/2.0/epsr/epsz/mur/muz; Bc(Y(X$PK  
D5=2.0*epsr*epsz; B;M?,<%FRU  
D6=2.0*epsr*epsz; x6BuF_.  
vUN22;Z\  
%*********************************************************************** 3-[q4R  
%     Initialize main grid update coefficients L5Ebc#  
%*********************************************************************** 8NxM4$nQX  
b+%f+zz*h  
C1ex(:,jh_bc:jh_tot-upml,:)=C1;     :y
+2*lV  
C2ex(:,jh_bc:jh_tot-upml,:)=C2; ;ic3).H  
C3ex(:,:,kh_bc:kh_tot-upml)=C3; 34`'M+3  
C4ex(:,:,kh_bc:kh_tot-upml)=C4; ?Y
$JWEPJ  
C5ex(ih_bc:ie_tot-upml,:,:)=C5; P2NQHX  
C6ex(ih_bc:ie_tot-upml,:,:)=C6; iJ-23_D  
!<j'Ea  
C1ey(:,:,kh_bc:kh_tot-upml)=C1; UvJ}b  
C2ey(:,:,kh_bc:kh_tot-upml)=C2; VQ(jpns5  
C3ey(ih_bc:ih_tot-upml,:,:)=C3; 4QH3fTv   
C4ey(ih_bc:ih_tot-upml,:,:)=C4; B\=L3eL<D  
C5ey(:,jh_bc:je_tot-upml,:)=C5; Vy6qbC-Kt  
C6ey(:,jh_bc:je_tot-upml,:)=C6; p#&h=,
W}  
E-4b[xNj*+  
C1ez(ih_bc:ih_tot-upml,:,:)=C1; lsJSYJG&  
C2ez(ih_bc:ih_tot-upml,:,:)=C2; sQ=]NF)\  
C3ez(:,jh_bc:jh_tot-upml,:)=C3; oQLq&zRH`f  
C4ez(:,jh_bc:jh_tot-upml,:)=C4; sGi"rg#  
C5ez(:,:,kh_bc:ke_tot-upml)=C5; AS!?q  
C6ez(:,:,kh_bc:ke_tot-upml)=C6; 9Z|jxy  
>W%Em
nLK  
D1hx(:,jh_bc:je_tot-upml,:)=D1; ?ME6+Z\  
D2hx(:,jh_bc:je_tot-upml,:)=D2; P9GN}GN%v  
D3hx(:,:,kh_bc:ke_tot-upml)=D3; _Us#\+]_:  
D4hx(:,:,kh_bc:ke_tot-upml)=D4; u7Y WnD  
D5hx(ih_bc:ih_tot-upml,:,:)=D5; |WQBDB`W  
D6hx(ih_bc:ih_tot-upml,:,:)=D6; $ZUdT  
KxI&G%z  
D1hy(:,:,kh_bc:ke_tot-upml)=D1; ot,jp|N>f~  
D2hy(:,:,kh_bc:ke_tot-upml)=D2; Tre]"2l  
D3hy(ih_bc:ie_tot-upml,:,:)=D3; v]'ztFA  
D4hy(ih_bc:ie_tot-upml,:,:)=D4; iNWw;_|1  
D5hy(:,jh_bc:jh_tot-upml,:)=D5; %D UH
@j  
D6hy(:,jh_bc:jh_tot-upml,:)=D6; ?5+.`L9H  
wu7Lk3  
D1hz(ih_bc:ie_tot-upml,:,:)=D1; $3W;=Id=+  
D2hz(ih_bc:ie_tot-upml,:,:)=D2; Pnk5mK$  
D3hz(:,jh_bc:je_tot-upml,:)=D3; AV:hBoO  
D4hz(:,jh_bc:je_tot-upml,:)=D4; xmNB29#  
D5hz(:,:,kh_bc:kh_tot-upml)=D5; /@wg>&L]  
D6hz(:,:,kh_bc:kh_tot-upml)=D6; f~t:L,\,  
JUsQ,ETn  
%*********************************************************************** i
/65v  
%     Fill in PML regions 9/{(%XwX  
% _ q(ko/T  
%     PML theory describes a continuous grading of the material properties ?[VM6- &  
%     over the PML region.  In the FDTD grid it is necessary to discretize -j+UMlkB  
%     the grading by averaging the material properties over a grid cell bbtGXfI+SB  
%     width centered on each field component.  As an example of the )YYf1o[+  
%     implementation of this averaging, we take the integral of the +dWDxguE{w  
%     continuous sigma(x) in the PML region XtXEB<4Z  
%   &VhroHO  
%         sigma_i = integral(sigma(x))/delta O%Scjm-^X  
%   +
@A  
%     where the integral is over a single grid cell width in x, and is HeK/7IAqp  
%     bounded by x1 and x2.  Applying this to the polynomial grading of 7yG#Z)VE  
%     Equation 7.60a produces ED2a}Tt>Z  
% Nu0C;B66  
%         sigma_i = (x2^(m+1)-x1^(m+1))*sigmam/(delta*(m+1)*d^m) AhCW'.  
% $'0u|Xy`  
%     where sigmam is the maximum value of sigma as described by Equation dWM'fg  
%     7.62. H~oail{EQ  
%         ySk R>y  
%     The definitions of x1 and x2 depend on the position of the field 3YeG$^y"  
%     component within the grid cell.  We have either vV'EZ?  
% .\\DKh%  
%         x1 = (i-0.5)*delta,  x2 = (i+0.5)*delta 
x5|I  
%   qPWP&k  
%     or #Z}Rfk(~  
%   Htay-PB }  
%         x1 = (i)*delta,      x2 = (i+1)*delta ,QOG!T4  
% _A M*@|p,  
%     where i varies over the PML region. OtY`@\hy  
%   XpIklL7  
%*********************************************************************** dl.N.P7}4  
6 +Sxr  
rmax=exp(-16);  %desired reflection error, designated as R(0) in Equation 7.62 u9:`4b  
GSY(  
orderbc=4;      %order of the polynomial grading, designated as m in Equation 7.60a,b @y e4q.m  
[=x[ w70  
%   x-varying material properties ~<P 0]ju  
delbc=upml*delta; &qLf@1AD  
sigmam=-log(rmax)*(orderbc+1.0)/(2.0*eta*delbc); qsF<!'m7`  
%sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1.0)); q?,).x nN  
:Gv1?M  
kmax=1; R]Vt Y7}i,  
kfactor=(kmax-1.0)/delta/(orderbc+1.0)/delbc^orderbc;
x[4`fM.m*  
%&$Tz1"  
for i=1:upml qBU-~"2t  
     LnFdhrB@x  
    % Coefficients for field components in the center of the grid cell 5(Cl1Yse=r  
    x1=(upml-i+1)*delta; OGBHos  
    x2=(upml-i)*delta; g6W)4cC8a  
    sigma=sigfactor*(x1^(orderbc+1)-x2^(orderbc+1)); 4]rnY~  
    ki=1+kfactor*(x1^(orderbc+1)-x2^(orderbc+1)); HvUxsdT  
    facm=(2*epsr*epsz*ki-sigma*dt); `^91%f  
    facp=(2*epsr*epsz*ki+sigma*dt); A]y`7jJ  
"MxnFeLM#  
    C5ex(i,:,:)=facp; Xk:OL,c  
    C5ex(ie_tot-i+1,:,:)=facp; =AsEZ)" _  
    C6ex(i,:,:)=facm; cO-7ke  
    C6ex(ie_tot-i+1,:,:)=facm; osciZ'~  
    D1hz(i,:,:)=facm/facp; y6@0O%TDN  
    D1hz(ie_tot-i+1,:,:)=facm/facp;
k=2Lo  
    D2hz(i,:,:)=2.0*epsr*epsz*dt/facp; 8{Q<N%Jnu  
    D2hz(ie_tot-i+1,:,:)=2.0*epsr*epsz*dt/facp; Om\o#{D  
    D3hy(i,:,:)=facm/facp; t^VwR=i  
    D3hy(ie_tot-i+1,:,:)=facm/facp; ,c$,!.r  
    D4hy(i,:,:)=1.0/facp/mur/muz; %`k6w3qI  
    D4hy(ie_tot-i+1,:,:)=1.0/facp/mur/muz; Q.bXM?V)  
@(
l^]9(V\  
    % Coefficients for field components on the grid cell boundary H12Fw'2  
    x1=(upml-i+1.5)*delta; r~[Ia!U?  
    x2=(upml-i+0.5)*delta; ~aw.(A?MI  
    sigma=sigfactor*(x1^(orderbc+1)-x2^(orderbc+1)); sD<a+Lw}x  
    ki=1.0+kfactor*(x1^(orderbc+1)-x2^(orderbc+1)); SEORSS  
    facm=(2.0*epsr*epsz*ki-sigma*dt); fTzvmC:g7  
    facp=(2.0*epsr*epsz*ki+sigma*dt); -1Jg?cPzk  
1ofKt=|=  
    C1ez(i,:,:)=facm/facp; l_3`G-`2  
    C1ez(ih_tot-i+1,:,:)=facm/facp; "B8Q:  
    C2ez(i,:,:)=2.0*epsr*epsz*dt/facp; tWo{7)Eb  
    C2ez(ih_tot-i+1,:,:)=2.0*epsr*epsz*dt/facp; cD@(/$wt  
    C3ey(i,:,:)=facm/facp; 909?_v  
    C3ey(ih_tot-i+1,:,:)=facm/facp; KTK <gV9:  
    C4ey(i,:,:)=1.0/facp/epsr/epsz; Gk967pC  
    C4ey(ih_tot-i+1,:,:)=1.0/facp/epsr/epsz; 6~OoFm5  
    D5hx(i,:,:)=facp; 4pe'06:  
    D5hx(ih_tot-i+1,:,:)=facp; < |e,05aM  
    D6hx(i,:,:)=facm; K7$x<5+)  
    D6hx(ih_tot-i+1,:,:)=facm; ~Xr=4V:a+  
     $v,dz_O*\  
end <DpevoF  
2`.cK 3  
%   PEC walls d[r#-h>dS  
C1ez(1,:,:)=-1.0; ?xK8#  
C1ez(ih_tot,:,:)=-1.0; XV!6dh!  
C2ez(1,:,:)=0.0; b>_o xK  
C2ez(ih_tot,:,:)=0.0; S(QpM.9*  
C3ey(1,:,:)=-1.0; z+x\(/  
C3ey(ih_tot,:,:)=-1.0; U!T~!C^  
C4ey(1,:,:)=0.0; TPVVck-T8  
C4ey(ih_tot,:,:)=0.0; KjV:|  
M]
<?k]_p  
%   y-varying material properties mrTlXXz  
delbc=upml*delta; UsgK  
sigmam=-log(rmax)*epsr*epsz*cc*(orderbc+1.0)/(2.0*delbc); ,/[6e\0~  
sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1.0)); ^*S ,xP  
kmax=1.0; 9s_vL9u  
kfactor=(kmax-1.0)/delta/(orderbc+1.0)/delbc^orderbc; ADZ};:]  
=5aDM\L$&  
for j=1:upml ~7Y+2FZ  
     g"Ljm7  
    % Coefficients for field components in the center of the grid cell PiYY6i0  
    y1=(upml-j+1)*delta; DvME1]7)  
    y2=(upml-j)*delta; Kfm5i Q  
    sigma=sigfactor*(y1^(orderbc+1)-y2^(orderbc+1)); P D4Tz!F  
    ki=1+kfactor*(y1^(orderbc+1)-y2^(orderbc+1)); avjpA?Vz  
    facm=(2*epsr*epsz*ki-sigma*dt); 4B=2>k  
    facp=(2*epsr*epsz*ki+sigma*dt); k $M]3}$U  
     -7m:91x  
    C5ey(:,j,:)=facp; 9Kr+\F  
    C5ey(:,je_tot-j+1,:)=facp; UYFwS/ RW}  
    C6ey(:,j,:)=facm; =b38(\  
    C6ey(:,je_tot-j+1,:)=facm; Y_}mYvJW  
    D1hx(:,j,:)=facm/facp; mt9.x  
    D1hx(:,je_tot-j+1,:)=facm/facp; nJbtS#`G4  
    D2hx(:,j,:)=2*epsr*epsz*dt/facp; ygOd69  
    D2hx(:,je_tot-j+1,:)=2*epsr*epsz*dt/facp; $`APHjijN  
    D3hz(:,j,:)=facm/facp; ktI/3Mb@  
    D3hz(:,je_tot-j+1,:)=facm/facp; W>!_|[a  
    D4hz(:,j,:)=1/facp/mur/muz; '"y|p+=j:  
    D4hz(:,je_tot-j+1,:)=1/facp/mur/muz; jQk*8  
     D@G\7KH@  
    % Coefficients for field components on the grid cell boundary f @8mS   
    y1=(upml-j+1.5)*delta;  R=.4  
    y2=(upml-j+0.5)*delta; ttXXy3G#  
    sigma=sigfactor*(y1^(orderbc+1)-y2^(orderbc+1)); ^ K|;~}P  
    ki=1+kfactor*(y1^(orderbc+1)-y2^(orderbc+1)); .id)VF-l  
    facm=(2*epsr*epsz*ki-sigma*dt); ]FD'5
p{  
    facp=(2*epsr*epsz*ki+sigma*dt);     U QE qX  
     1D16  
    C1ex(:,j,:)=facm/facp; ag$Vgl  
    C1ex(:,jh_tot-j+1,:)=facm/facp; Z:ni$7<.  
    C2ex(:,j,:)=2*epsr*epsz*dt/facp; 8iW;y2qF  
    C2ex(:,jh_tot-j+1,:)=2*epsr*epsz*dt/facp; ^Y<|F!0  
    C3ez(:,j,:)=facm/facp; Fd?"-  
    C3ez(:,jh_tot-j+1,:)=facm/facp; -<Hu!V`+  
    C4ez(:,j,:)=1/facp/epsr/epsz; YRv&1!VLE  
    C4ez(:,jh_tot-j+1,:)=1/facp/epsr/epsz;   [qdRUV'  
    D5hy(:,j,:)=facp; Jm|+-F@I  
    D5hy(:,jh_tot-j+1,:)=facp; CQZgMY1{  
    D6hy(:,j,:)=facm; ?!wgH9?8  
    D6hy(:,jh_tot-j+1,:)=facm; &P.4(1sC  
x??pBhJH  
end jlp:lX  
~UyV<  
%   PEC walls D*Ik7Pe  
C1ex(:,1,:)=-1; $f,n8]  
C1ex(:,jh_tot,:)=-1; iY`%SmB  
C2ex(:,1,:)=0; *J$=.fF1  
C2ex(:,jh_tot,:)=0; 9k9_mjLZ  
C3ez(:,1,:)=-1; dp++%:j  
C3ez(:,jh_tot,:)=-1; nM\eDNK  
C4ez(:,1,:)=0; VmCW6 G#M  
C4ez(:,jh_tot,:)=0;   1h
>yu3O  
UUF;p2{f  
%   z-varying material properties u583_k%  
delbc=upml*delta; RbCPmiZcH  
sigmam=-log(rmax)*epsr*epsz*cc*(orderbc+1)/(2*delbc); DKfE.p)  
sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1)); z?>D_NLX6  
%kmax=1; h\7fp.  
%kfactor=(kmax-1)/delta/(orderbc+1)/delbc^orderbc; J\J?yo 6  
kfactor=0; _tSAI  
vgD {qg@  
for k=1:upml ;GVV~.7/  
1J6,]M  
    % Coefficients for field components in the center of the grid cell .U"8mP=&  
    z1=(upml-k+1)*delta; G+F#n6Vx  
    z2=(upml-k)*delta; !E,A7s  
    sigma=sigfactor*(z1^(orderbc+1)-z2^(orderbc+1)); O_cbP59Y.  
    ki=1+kfactor*(z1^(orderbc+1)-z2^(orderbc+1)); {*[\'!d--.  
    facm=(2*epsr*epsz*ki-sigma*dt); Qhs/E`k4  
    facp=(2*epsr*epsz*ki+sigma*dt); $p#%G#T  
     Q9Uf.Lh2  
    C5ez(:,:,k)=facp; DjI3?NN  
    C5ez(:,:,ke_tot-k+1)=facp; HQ|MhM/"  
    C6ez(:,:,k)=facm; 76wc,+  
    C6ez(:,:,ke_tot-k+1)=facm; xBUya4w  
    D1hy(:,:,k)=facm/facp; L,SGT
8lL  
    D1hy(:,:,ke_tot-k+1)=facm/facp; )>b.;  
    D2hy(:,:,k)=2*epsr*epsz*dt/facp; V?Z.\~
 
    D2hy(:,:,ke_tot-k+1)=2*epsr*epsz*dt/facp; 0E?jW7yr  
    D3hx(:,:,k)=facm/facp; J
o$G,Q  
    D3hx(:,:,ke_tot-k+1)=facm/facp; C|d\3S\(  
    D4hx(:,:,k)=1/facp/mur/muz; ikSF)r;*t  
    D4hx(:,:,ke_tot-k+1)=1/facp/mur/muz; ET_W-  
     9Rn? :B~W:  
    % Coefficients for field components on the grid cell boundary ^M%uV  
    z1=(upml-k+1.5)*delta; DVah
 
    z2=(upml-k+0.5)*delta; M~WijDj  
    sigma=sigfactor*(z1^(orderbc+1)-z2^(orderbc+1)); w$H^q
!(  
    ki=1+kfactor*(z1^(orderbc+1)-z2^(orderbc+1)); H~GQ;PhRx  
    facm=(2*epsr*epsz*ki-sigma*dt); SF}<{x_  
    facp=(2*epsr*epsz*ki+sigma*dt); a
\IP12F?  
     _ ):d`O e  
    C1ey(:,:,k)=facm/facp; Q?8R[i  
    C1ey(:,:,kh_tot-k+1)=facm/facp; F.-R r  
    C2ey(:,:,k)=2*epsr*epsz*dt/facp; umF Z?a  
    C2ey(:,:,kh_tot-k+1)=2*epsr*epsz*dt/facp; nt;haeJ  
    C3ex(:,:,k)=facm/facp; ebS0qo[oLH  
    C3ex(:,:,kh_tot-k+1)=facm/facp; %YSpCI  
    C4ex(:,:,k)=1/facp/epsr/epsz; ixW@7m  
    C4ex(:,:,kh_tot-k+1)=1/facp/epsr/epsz; :@1eph0  
    D5hz(:,:,k)=facp; smdZxFl  
    D5hz(:,:,kh_tot-k+1)=facp; GiP`dtK
 
    D6hz(:,:,k)=facm; \%/#x V  
    D6hz(:,:,kh_tot-k+1)=facm; (Fynok  
c~{9a_G  
end E Q4KV  
y_*PQZ$c<  
%   PEC walls BYO"u6  
C1ey(:,:,1)=-1; `Ja?fI'H-  
C1ey(:,:,kh_tot)=-1; /CuXa%Ci^  
C2ey(:,:,1)=0; bV edFm  
C2ey(:,:,kh_tot)=0; lY~4'8^  
C3ex(:,:,1)=-1; cE`6uq7p  
C3ex(:,:,kh_tot)=-1;  %ObLWH'  
C4ex(:,:,1)=0; AZzuI*  
C4ex(:,:,kh_tot)=0; )x}l3\s  
"jTKSgv+q5  
%figure )+6v  
%set(gcf,'DoubleBuffer','on') 6'kS_Zu{<  
d)@<W1;  
%*********************************************************************** $e\h}A6  
%     Begin time stepping loop ~/8M 3k/  
%*********************************************************************** S
-7'it!1  
nB%;S  
for n=1:nmax K TsgJ\W  
     G2]4n T  
    % Update magnetic field gXonF'  
    bstore=bx; d0aCY  
    bx(2:ie_tot,:,:)=D1hx(2:ie_tot,:,:).*  bx(2:ie_tot,:,:)-... eh4gQ^l  
                     D2hx(2:ie_tot,:,:).*((ez(2:ie_tot,2:jh_tot,:)-ez(2:ie_tot,1:je_tot,:))-... rj6tZJZ#o0  
                                          (ey(2:ie_tot,:,2:kh_tot)-ey(2:ie_tot,:,1:ke_tot)))./delta; '" X_B0k  
    hx(2:ie_tot,:,:)= D3hx(2:ie_tot,:,:).*hx(2:ie_tot,:,:)+... >$
NDv  
                      D4hx(2:ie_tot,:,:).*(D5hx(2:ie_tot,:,:).*bx(2:ie_tot,:,:)-... ddfs8\  
                                           D6hx(2:ie_tot,:,:).*bstore(2:ie_tot,:,:)); 6ZKsz5:=  
    bstore=by; f6_];]yP  
    by(:,2:je_tot,:)=D1hy(:,2:je_tot,:).*  by(:,2:je_tot,:)-... ?lbH02P{v  
                     D2hy(:,2:je_tot,:).*((ex(:,2:je_tot,2:kh_tot)-ex(:,2:je_tot,1:ke_tot))-... is1's[  
                                          (ez(2:ih_tot,2:je_tot,:)-ez(1:ie_tot,2:je_tot,:)))./delta; L7= Q<D<  
    hy(:,2:je_tot,:)= D3hy(:,2:je_tot,:).*hy(:,2:je_tot,:)+... ! iptT(2  
                      D4hy(:,2:je_tot,:).*(D5hy(:,2:je_tot,:).*by(:,2:je_tot,:)-... s[K^9wz  
                                           D6hy(:,2:je_tot,:).*bstore(:,2:je_tot,:)); -6tgsfEr  
    bstore=bz; 4Ue_Y'LmM  
    bz(:,:,2:ke_tot)=D1hz(:,:,2:ke_tot).*  bz(:,:,2:ke_tot)-... :Xn7Ha[f  
                     D2hz(:,:,2:ke_tot).*((ey(2:ih_tot,:,2:ke_tot)-ey(1:ie_tot,:,2:ke_tot))-... %r-V2)  
                                          (ex(:,2:jh_tot,2:ke_tot)-ex(:,1:je_tot,2:ke_tot)))./delta; M t*6}Cl  
    hz(:,:,2:ke_tot)= D3hz(:,:,2:ke_tot).*hz(:,:,2:ke_tot)+... `((Y
c]:7  
                      D4hz(:,:,2:ke_tot).*(D5hz(:,:,2:ke_tot).*bz(:,:,2:ke_tot)-... e$u4vC~  
                                           D6hz(:,:,2:ke_tot).*bstore(:,:,2:ke_tot)); ?gO8kPg/D  
     +$$$  
    % Update electric field DHw
&+MY  
    dstore=dx; f'<Q.Vh<  
    dx(:,2:je_tot,2:ke_tot)=C1ex(:,2:je_tot,2:ke_tot).*  dx(:,2:je_tot,2:ke_tot)+... L lw&& K  
                            C2ex(:,2:je_tot,2:ke_tot).*((hz(:,2:je_tot,2:ke_tot)-hz(:,1:je_tot-1,2:ke_tot))-... 9Ro6fjjE  
                                                        (hy(:,2:je_tot,2:ke_tot)-hy(:,2:je_tot,1:ke_tot-1)))./delta;
^K7ic,{  
    ex(:,2:je_tot,2:ke_tot)=C3ex(:,2:je_tot,2:ke_tot).*ex(:,2:je_tot,2:ke_tot)+... 6*qL[m.F[o  
                            C4ex(:,2:je_tot,2:ke_tot).*(C5ex(:,2:je_tot,2:ke_tot).*dx(:,2:je_tot,2:ke_tot)-... N0K){  
                                                        C6ex(:,2:je_tot,2:ke_tot).*dstore(:,2:je_tot,2:ke_tot)); y.=/J8->  
    dstore=dy; R*oXmuOsYA  
    dy(2:ie_tot,:,2:ke_tot)=C1ey(2:ie_tot,:,2:ke_tot).*  dy(2:ie_tot,:,2:ke_tot)+... :Gu+m  
                            C2ey(2:ie_tot,:,2:ke_tot).*((hx(2:ie_tot,:,2:ke_tot)-hx(2:ie_tot,:,1:ke_tot-1))-... 21ppSN>  
                                                        (hz(2:ie_tot,:,2:ke_tot)-hz(1:ie_tot-1,:,2:ke_tot)))./delta;  QV h4  
    ey(2:ie_tot,:,2:ke_tot)=C3ey(2:ie_tot,:,2:ke_tot).*ey(2:ie_tot,:,2:ke_tot)+... \S*$UE]uG  
                            C4ey(2:ie_tot,:,2:ke_tot).*(C5ey(2:ie_tot,:,2:ke_tot).*dy(2:ie_tot,:,2:ke_tot)-... i
]=&  
                                                        C6ey(2:ie_tot,:,2:ke_tot).*dstore(2:ie_tot,:,2:ke_tot)); b{d4xU8'  
    dstore=dz; xXY.AoO6  
    dz(2:ie_tot .. d{3@h+zL  
Q~MC7-n>  

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

   你的第一个问题,我当时也看了半天,后来对照着taflove的2D程序发现,PML的设置也是这么写的,后来我就默认了,具体怎么弄我想期待高人的解答 B_ja&) !s1  
  你的第二问题你可以参考herosward师兄翻译的tafloveUPML章节,论坛有你可以去下载,根据UPML的递推公式,c1ex这些都是公式的参数,我们已c1ex为例: Uu"0rUzt  
c1ex=(2*epsr*ky-sigma*dt)/(2*epsr*ky+sigma*dt) \ A%eG&  
ZUp\Ep}  
你可以看到这个参数只跟y方向的设置有关,于x,z方向的网格设置无关,你初始化总场的系数,当然要把y方向范围设置成总场的范围,c1ex(:,jh_bc:jh_tot-upml,:)=c1中jh_bc:jh_tot-upml就是y方向总场的范围 (f_g7B2&y  
希望能对你有所帮助,同时期待高手的补充

关于第一个问题的，采用7.60作为PML的电导系数变量的设计，\sigma (x_i) 不仅仅是在x_i点上衡量7.60.而是对\sigma(x)在一个胞上做积分，在除以胞的边长。如果你对7.60做积分，会得到TAFLOVE用的表达式。

恩，看懂了，谢谢指点~

明白了，可是他用这种做法是为了更精确点吗？

可能是尽量匀化sigma从而减少数值上的反射吧

呵呵，发现一个相同的问题 :D^Y?  
没看清帖子，重新编辑，删掉

kmax=1; 4K% YS   ,*q#qW!!  
kfactor=(kmax-1.0)/delta/(orderbc+1.0)/delbc^orderbc;  5 Ep  我觉得这里好像没有必要啊.不是一直为0吗,请指教!

其实我觉得第一个问题，完全可以这么写， pDLu+}@  
C1ex(:,:,:)=C1;     c+,7Zu!  
C2ex(:,:,:)=C2; $V`KrA~]  
C3ex(:,:,:)=C3; &=+cov(3  
C4ex(:,:,:)=C4; ]Ssw32yn  
C5ex(:,:,:)=C5; +cPE4(d  
C6ex(:,:,:)=C6; PK:o}IWn~x  
C7ex(:,:,:)=C7; 3p?<iVE  
C8ex(:,:,:)=C8; i6!T`Kau  
不会影响结果，这样写还快速，不用去考虑哪里是PML，哪里是主机算区域。因为后面还是要给部分重新赋值，我先赋值了也没关系，后面会覆盖。