[已解决]请牛人指教一下，三维程序中均匀平面波入射时连接边界这样加入有问题吗？

里面的计算过程根据葛老师编的书上公式！（前4个式子计算E0，后面两个坐标转换，不知道这样对不对）请教牛人指导，磁场也按同样方法加入！ ^jE8+h  
ie为x方向格子数+1，je，ke分别为y、z方向 jXg  
left_x(i,1,k)存储左侧面ex值，其他的类似！ An`3Ex[  
pB:$lS  
for i=1:ie G}d-(X
 
        for k=1:ke OO) ~HV4\  
            temp1(1)=i*sin(theta)*cos(psi)+jets*cos(theta)*sin(psi)+k*cos(psi); J([s5:.[  
            temp1(2)=floor(temp1(1));                                %左表面  电场 f9u^R=Ff[  
            temp1(3)=temp1(1)-temp1(2); eU@Cr7@,|  
            temp1(4)=temp1(2)+1; BUZ _)  
            left_x(i,1,k)=(1-temp1(3))*ein(temp1(2))+temp1(3)*ein(temp1(4)); YDJ4c;37  
            left_x(i,1,k)=left_x(i,1,k)*(-sin(psi)*sin(alpha)+cos(theta)*cos(psi)*cos(alpha)); 0&+k.Vg  
            left_z(i,1,k)=-left_x(i,1,k)*sin(theta)*cos(alpha); PmpNAVE'  
            g"VMeW^  
            temp1(1)=i*sin(theta)*cos(psi)+jete*cos(theta)*sin(psi)+k*cos(psi); IM@tN L  
            temp1(2)=floor(temp1(1));                                %右边面  电场 <Zb/  
            temp1(3)=temp1(1)-temp1(2); _
fk#<  
            temp1(4)=temp1(2)+1; o{:xp r=(  
            right_x(i,1,k)=(1-temp1(3))*ein(temp1(2))+temp1(3)*ein(temp1(4)); d3Mva,bw<  
            right_x(i,1,k)=right_x(i,1,k)*(-sin(psi)*sin(alpha)+cos(theta)*cos(psi)*cos(alpha)); <O<LYN+(  
            right_z(i,1,k)=-right_x(i,1,k)*sin(theta)*cos(alpha); _qwQ;!9  
        end / ~%KVe  
    end NpP')m!`}  
        r,1e 'd:  
    for i=1:ie YIRZ+H<Q  
        for j=1:je r=uN9ro  
            temp1(1)=i*sin(theta)*cos(psi)+j*cos(theta)*sin(psi)+kets*cos(psi); 8IQtz2  
            temp1(2)=floor(temp1(1));                                %下边面  电场 pKUP2m`MW  
            temp1(3)=temp1(1)-temp1(2); Uz7oL8  
            temp1(4)=temp1(2)+1; g|X;ahTT  
            bottom_x(i,j,1)=(1-temp1(3))*ein(temp1(2))+temp1(3)*ein(temp1(4)); ?%tMohL  
            bottom_x(i,j,1)=bottom_x(i,j,1)*(-sin(psi)*sin(alpha)+cos(theta)*cos(psi)*cos(alpha)); =wWpP-J&  
            bottom_y(i,j,1)=bottom_x(i,j,1)*(cos(psi)*sin(alpha)+cos(theta)*sin(psi)*cos(alpha)); Dim> 7Wbh  
            YY((#"o;l  
            temp1(1)=i*sin(theta)*cos(psi)+j*cos(theta)*sin(psi)+kete*cos(psi); thlY0XCq,%  
            temp1(2)=floor(temp1(1));                                %上边面  电场 7cDU
2l  
            temp1(3)=temp1(1)-temp1(2); rqPo)AL  
            temp1(4)=temp1(2)+1; op2Of<{h  
            top_x(i,j,1)=(1-temp1(3))*ein(temp1(2))+temp1(3)*ein(temp1(4)); y9H% Xl  
            top_x(i,j,1)=top_x(i,j,1)*(-sin(psi)*sin(alpha)+cos(theta)*cos(psi)*cos(alpha)); OR1DYHHT/1  
            top_y(i,j,1)=top_x(i,j,1)*(cos(psi)*sin(alpha)+cos(theta)*sin(psi)*cos(alpha)); $ ,Ck70_  
        end e}Vw!w  
    end ;*TIM%6#  
+n0r0:z0  
    for j=1:je *|.0Myjo  
        for k=1:ke K)tQ]P  
            temp1(1)=iets*sin(theta)*cos(psi)+j*cos(theta)*sin(psi)+k*cos(psi); ,_.I\
EY[  
            temp1(2)=floor(temp1(1));                                %后边面  电场 =ac_,]z  
            temp1(3)=temp1(1)-temp1(2); 0oZsb
\  
            temp1(4)=temp1(2)+1; (IqZ@->nw  
            back_y(1,j,k)=(1-temp1(3))*ein(temp1(2))+temp1(3)*ein(temp1(4)); 6RO(]5wX  
            back_y(1,j,k)=back_y(1,j,k)*(cos(psi)*sin(alpha)+cos(theta)*sin(psi)*cos(alpha)); {P%9  
            back_z(1,j,k)=-back_y(1,j,k)*sin(theta)*cos(alpha); S{t+
>/  
            #p(h]T32  
            temp1(1)=iete*sin(theta)*cos(psi)+j*cos(theta)*sin(psi)+k*cos(psi); K7knK  
            temp1(2)=floor(temp1(1));                           
.. Qd_6)M-  
6 gL=u-2  

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

从运行结果来看，不正确，不知道是不是这一步有问题？做过这样程序的一定指教一下小弟！

自己顶一下！ e]jH+IR:>  
书中的计算r'的公式是不是错了？其中中间的一项中cos（theta）是不是要改为sin（theta）？又知道的嘛？自己推算的时候感觉书上的错了，有不敢冒昧！望牛人指点！

中间一项本来就是sin（theta）

            left_x(i,1,k)=(1-temp1(3))*ein(temp1(2))+temp1(3)*ein(temp1(4)); N}\[Gr  
            left_x(i,1,k)=left_x(i,1,k)*(-sin(psi)*sin(alpha)+cos(theta)*cos(psi)*cos(alpha)); >J|]moSVA  
            left_z(i,1,k)=-left_x(i,1,k)*sin(theta)*cos(alpha); TYI7<-Mp:[  
  }K8/-d6  
貌似这里也不对，left_x(i,1,k)好像在第二步的时候值已经改变了。不应该用这个迭代吧！

谢了，今天一天确实找了很多毛病，不过你提的这个错误，还没发现呢！茅塞顿开了！不知该说啥了，谢谢您了！o(∩_∩)o...明天再调调试试！

明天就回家了，再来顶顶！

有空看看。

我曾经推导过，我想看看你的程序吧

clear dMx4ykrR  
clc (eCFWmO  
%*************************************** HmK*bZ  
%基本参数 bU 63X={  
%**************************************** a'\By?V]  
cz=2.99792456e8; _B2V "p  
murz=4.0*pi*1.0e-7; 8&"(WuZ@  
epsz=1.0/(cz*cz*murz); vFrt|JC_{  
freq0=1.15*1.0e10;            %激发源的频率 8 6QE/M  
omega0=2*pi*freq0;            %入射波的角速度 U-wLt(Y<  
lambda=cz/freq0;              %激发源中心波长 TaJB4zB  
psi=45*pi/180;                %入射波方向在xoy平面内与x轴的夹角 t"j|nz{m  
theta=45*pi/180;              %入射波方向与z轴夹角 1G~S|,8p  
alpha=45*pi/180;              %入射波为线性极化，alpha为极化角 $V6^G*Q  
incidentstart=1;              %入射波开始位置 zT~B6  
incidentend=1000; zm9TvoC%}  
source=20; :PbDU$x  
t=2*pi/omega0;                %入射波周期 'cDx{?  
z=(murz/epsz)^0.5;            %特征阻抗 Rd+P,PO  
%*********************************************** 3>z[
PPw  
%格子参数 pO<-.,  
%*********************************************** hM@\RPsY  
ie=100;                     *5$&`&,  
je=100;                     wxLXh6|6%_  
ke=100;                         \F,DA"K_  
ib=ie-1;                     %x方向格子大小 OEr:xK2T  
jb=je-1;                     %y方向格子大小 KRsAv^']  
kb=ke-1;                     %z方向格子大小 .'[/|4H  
"%8A:^1  
iet=80;                     %x方向总场大小 gJ2 H=#M  
jet=80;                     %y方向总场大小 NkxCs  
ket=80;                     %z方向总场大小 ^C&+ ~+  
ibt=iet-1; V's:>;  
jbt=iet-1; j*T]HaM  
kbt=ket-1; 0JRD  
M4Z@O3OIE  
iets=(ie-iet)/2+1;          %x轴方向,总场的开始位置 q&'Lbxc>c  
iete=ie-(ie-iet)/2;         %x轴方向,总场的结束位置   !/K8xD$  
jets=(je-jet)/2+1;          %y轴方向,总场的开始位置 lhC6S'vq  
jete=je-(je-jet)/2;         %y轴方向,总场的结束位置 >zY \Llv  
kets=(ke-ket)/2+1;          %z轴方向,总场的开始位置 ]Pn!nSg  
kete=ke-(ke-ket)/2;         %z轴方向,总场的结束位置 C-P06Q]  
09M;}4ev&7  
dx=lambda/100; bAxTLIf  
dt=0.95*dx/(cz*2^0.5); WKA'=,`v  
ca=1; LE?u`i,e=+  
cb=dt/epsz; @D`zKYwX1  
cp=1; -U2mfW  
cq=dt/murz; PM$Ee #62R  
tmax=1000; `29TY&p+"  
%******************************************* s_jBu  
%场的初始场强
V9x8R  
%******************************************* Iy {U'a!  
ex=zeros(ib,je,ke);           %ex场 2^exL h  
ey=zeros(ie,jb,ke); iZ[tHw||  
ez=zeros(ie,je,kb); MrE<vw@he  
hx=zeros(ie,jb,kb); 
NnxM3*  
hy=zeros(ib,je,kb); jj\[
7 O*  
hz=zeros(ib,jb,ke); "C%!8`K{a*  
ein=zeros((incidentend-incidentstart+1),1); guN4-gGDr<  
hin=zeros((incidentend-incidentstart),1); ]0c Pml  
temp=zeros(2,1);                    %处理一维的入射波   aI#4H+/  
temp1=zeros(4,1);                   %处理连接边界投影 b`cYpcs  
back_y=zeros(3,je,ke);              %处理后表面的连接边界，其中1为电场2为磁场3临时时储存r’以下同 %IpSK 0<Sp  
back_z=zeros(2,je,ke); J."{<&  
front_y=zeros(3,je,ke); oL/o*^  
front_z=zeros(2,je,ke); hQJWKAf,/  
left_z=zeros(ie,2,ke); MBk"KF  
left_x=zeros(ie,3,ke); Wf02$c0#K  
right_z=zeros(ie,2,ke); w'Z!;4E0  
right_x=zeros(ie,3,ke); 7x.%hRk  
top_y=zeros(ie,je,2); ZofHic  
top_x=zeros(ie,je,3); o}8{Bh^  
bottom_y=zeros(ie,je,2); 1TqF6`;+  
bottom_x=zeros(ie,je,3); +I|8Q|^SD  
>3;^l/2c  
.")b?#K  
DOD6Liau{Q  
%*********************************************** c|wCKn}`  
% Movie initialization TW`mxj_J2  
%*********************************************** nYv#4*  
tview(:,:)=ez(:,:,ke/2); w`yx=i#  
sview(:,:)=ez(:,je/2,:); twqFs  
<D/al9  
subplot(2,1,1),pcolor(tview'); <Mgf]v.QS  
shading flat; ]S8LY.Az5  
caxis([-1.0 1.0]); &B-[oqC?  
colorbar; yYAnwf  
axis image; snk$^  
title(['Ez(i,j,k=50),timestep=0']); ~ }KzJiL  
xlabel('i coordinate'); Oo%!>!Lt,  
ylabel('j coordinate'); _=`x])mM  
bLG]Wa  
subplot(2,1,2),pcolor(sview'); 1czG55 |  
shading flat; sV0Z  
caxis([-1.0 1.0]); :q2YBa  
colorbar; +h[e0J|v{  
axis image; &R]pw`mTH  
title(['Ez(i,j=50,k),timestep=0']); xi {|  
xlabel('i coordinate'); ='/Z;3jt]x  
ylabel('k coordinate'); i|2$8G3  
y{rn-?`{  
rect=get(gcf,'Position'); I"!'AI-  
rect(1:2)=[0 0]; ":WYcaSi  
M=moviein(tmax/1,gcf,rect); 5ouQQ)vA  
%********************************************** Md(JIlh3  
% 开始时间的 2o{@nN8%  
%********************************************** 5 Rz/Ri\c=  
for n=1:tmax `ENP=kL(+  
    %********************************** +802`eax  
    % 产生入射波 m^$5K's&  
    %********************************** okBE|g  
    hin(incidentstart:incidentend-1)=hin(incidentstart:incidentend-1)-... HY;oy(  
        (ein(incidentstart+1:incidentend)-ein(incidentstart:incidentend-1))*dt/murz/dx; jW5iqU"{*  
    ein(incidentstart+1:incidentend-1)=ein(incidentstart+1:incidentend-1)-... O
.:I,D&]  
        (hin(incidentstart+1:incidentend-1)-hin(incidentstart:incidentend-2))*dt/epsz/dx; {*=E?oF@  
    ein(source,1)=1*sin(omega0*n*dt); hjY0w  
    ein(incidentstart)=temp(1)+(ein(incidentstart+1)-ein(incidentstart))*(cz*dt-dx)/(cz*dt+dx); ]UUI~sFE  
    temp(1)=ein(incidentstart+1); 9G:TW|)L[Q  
    ein(incidentend)=temp(2)+(ein(incidentend-1)-ein(incidentend))*(cz*dt-dx)/(cz*dt+dx); <*@~n- R$  
    temp(2)=ein(incidentend-1); nlfPg-78B+  
    %********************************** kJ8vKcc  
    % 将连接边界投影到一维的入射波上 CV^0.  
    %**********************************   9={N4}<  
    for i=1:ie hYvNcOSks  
        for k=1:ke n85r^W  
            temp1(1)=i*sin(theta)*cos(psi)+jets*sin(theta)*sin(psi)+k*cos(theta); Jirct,k  
            temp1(2)=floor(temp1(1));                                %左表面入射  电场 QaMDGD  
            temp1(3)=temp1(1)-temp1(2); #4y,a_)  
            temp1(4)=temp1(2)+1; !P|5#.eC  
            left_x(i,3,k)=(1-temp1(3))*ein(temp1(2))+temp1(3)*ein(temp1(4)); yKDZ+3xK]  
            left_x(i,1,k)=left_x(i,3,k)*(-sin(psi)*sin(alpha)+cos(theta)*cos(psi)*cos(alpha)); fcAIg(vW  
            left_z(i,1,k)=-left_x(i,3,k)*sin(theta)*cos(alpha); 7Jx%JgF  
             vj3isI4lU  
            temp1(1)=i*sin(theta)*cos(psi)+jete*sin(theta)*sin(psi)+k*cos(theta); |QYZRz  
            temp1(2)=floor(temp1(1));                                %右边面入射  电场 M'u=H  
            temp1(3)=temp1(1)-temp1(2); ^pcRW44K  
            temp1(4)=temp1(2)+1; F!CAitxd  
            right_x(i,3,k)=(1-temp1(3))*ein(temp1(2))+temp1(3)*ein(temp1(4)); /_OOPt=G  
            right_x(i,1,k)=right_x(i,3,k)*(-sin(psi)*sin(alpha)+cos(theta)*cos(psi)*cos(alpha)); ZO7bSxAN-  
            right_z(i,1,k)=-right_x(i,3,k)*sin(theta)*cos(alpha); QFzFL-H~N  
        end nUqy1(  
    end ,+-?Zv 2  
         UD*+"~
 
    for i=1:ie Pf8u/?/  
        for j=1:je XV2=8#R  
            temp1(1)=i*sin(theta)*cos(psi)+j*sin(theta)*sin(psi)+kets*cos(theta); 7Dl%UG]  
            temp1(2)=floor(temp1(1));                                %下边面入射  电场 w`#fH  
            temp1(3)=temp1(1)-temp1(2); 7zI5PGWw  
            temp1(4)=temp1(2)+1; {,f[r*{Y  
            bottom_x(i,j,3)=(1-temp1(3))*ein(temp1(2))+temp1(3)*ein(temp1(4)); UvD-C?u'  
            bottom_x(i,j,1)=bottom_x(i,j,3)*(-sin(psi)*sin(alpha)+cos(theta)*cos(psi)*cos(alpha)); ;QidDi_s>  
            bottom_y(i,j,1)=bottom_x(i,j,3)*(cos(psi)*sin(alpha)+cos(theta)*sin(psi)*cos(alpha)); zUQe0Gc.b^  
             ^gm>!-Gx  
            temp1(1)=i*sin(theta)*cos(psi)+j*sin(theta)*sin(psi)+kete*cos(theta); b7'F|h^  
            temp1(2)=floor(temp1(1));                                %上边面入射  电场 =h\E<dw  
            temp1(3)=temp1(1)-temp1(2); "-U3
=+  
            temp1(4)=temp1(2)+1; vXubY@k2  
            top_x(i,j,3)=(1-temp1(3))*ein(temp1(2))+temp1(3)*ein(temp1(4)); _<u;4RO(s  
            top_x(i,j,1)=top_x(i,j,3)*(-sin(psi)*sin(alpha)+cos(theta)*cos(psi)*cos(alpha)); t{ H1u  
            top_y(i,j,1)=top_x(i,j,3)*(cos(psi)*sin(alpha)+cos(theta)*sin(psi)*cos(alpha)); >ITEd  
        end Ygx,t|?7  
    end 4g!7 4a  
5+FLSk  
    for j=1:je 5Bd(>'ig_  
        for k=1:ke 9r8D*PvS  
            temp1(1)=iets*sin(theta)*cos(psi)+j*sin(theta)*sin(psi)+k*cos(theta); !Zj#
.6c9  
            temp1(2)=floor(temp1(1));                                %后边面入射  电场 1 ;Ju
]  
            temp1(3)=temp1(1)-temp1(2); 0#`)Prop6  
            temp1(4)=temp1(2)+1; Q [:<S/w  
            back_y(3,j,k)=(1-temp1(3))*ein(temp1(2))+temp1(3)*ein(temp1(4)); n!?r }n8  
            back_y(1,j,k)=back_y(3,j,k)*(cos(psi)*sin(alpha)+cos(theta)*sin(psi)*cos(alpha)); _?Ckq
 
            back_z(1,j,k)=-back_y(3,j,k)*sin(theta)*cos(alpha); Vi,Y@+4  
             ?waebuj>  
            temp1(1)=iete*sin(theta)*cos(psi)+j*sin(theta)*sin(psi)+k*cos(theta); ;.0LRWcJ  
            temp1(2)=floor(temp1(1));                                %前边面入射  电场 #11RLvDQd  
            temp1(3)=temp1(1)-temp1(2); hNVMz`r  
            temp1(4)=temp1(2)+1; WY.
5K =}  
            front_y(3,j,k)=(1-temp1(3))*ein(temp1(2))+temp1(3)*ein(temp1(4));  QT_^M1%  
            front_y(1,j,k)=front_y(3,j,k)*(cos(psi)*sin(alpha)+cos(theta)*sin(psi)*cos(alpha)); a>(~C
'(<  
            front_z(1,j,k)=-front_y(3,j,k)*sin(theta)*cos(alpha); L<E/,IdE  
        end @Z=wE3T@  
    end ;Xh5oB\)W  
     sy.:T]ZH  
     for i=1:ie \P@S"QO  
        for k=1:ke -$ali[  
            temp1(1)=i*sin(theta)*cos(psi)+(jets-1)*sin(theta)*sin(psi)+k*cos(theta); YE
_6OLW  
            temp1(2)=floor(temp1(1));                                %左表面入射  磁场 &E]"c]i+  
            temp1(3)=temp1(1)-temp1(2); <{ #<5 8  
            temp1(4)=temp1(2)+1; -R74/GBg  
            left_x(i,3,k)=(1-temp1(3))*hin(temp1(2))+temp1(3)*hin(temp1(4)); EOQaY  
            left_x(i,2,k)=left_x(i,3,k)*(-sin(psi)*cos(alpha)-cos(theta)*cos(psi)*sin(alpha));      %hx  j0-1 !=knppY  
            left_z(i,2,k)=left_x(i,3,k)*sin(theta)*sin(alpha);                                      %hz  j0-1                                       ~ a>S#S  
             t[j9R#02?  
            temp1(1)=i*sin(theta)*cos(psi)+(jete+1)*sin(theta)*sin(psi)+k*cos(theta); h&$
Py  
            temp1(2)=floor(temp1(1));                                %右边面入射  磁场 ~,G]glu8  
            temp1(3)=temp1(1)-temp1(2); LT&/0  
            temp1(4)=temp1(2)+1; }#|2z}!  
            right_x(i,3,k)=(1-temp1(3))*hin(temp1(2))+temp1(3)*hin(temp1(4)); Cg*kN"8q  
            right_x(i,2,k)=right_x(i,3,k)*(-sin(psi)*cos(alpha)-cos(theta)*cos(psi)*sin(alpha));   %hx   jb+1 }+JLn%H)  
            right_z(i,2,k)=right_x(i,3,k)*sin(theta)*sin(alpha);                                   %hz   jb+1 l]u
7.~b  
        end ]gHLcr3  
     end +WdL  
     h{H]xe[Q  
    for i=1:ie %Y 2G  
        for j=1:je [U']kt  
            temp1(1)=i*sin(theta)*cos(psi)+j*sin(theta)*sin(psi)+(kets-1)*cos(theta); /'"R Mq  
            temp1(2)=floor(temp1(1));                                %下边面入射  磁场 wK
LN:aRF2  
            temp1(3)=temp1(1)-temp1(2); /gX%ABmS  
            temp1(4)=temp1(2)+1; 43F^J
%G  
            bottom_x(i,j,3)=(1-temp1(3))*hin(temp1(2))+temp1(3)*hin(temp1(4));
d1lH[r!Z  
            bottom_x(i,j,2)=bottom_x(i,j,3)*(-sin(psi)*cos(alpha)-cos(theta)*cos(psi)*sin(alpha));  %hx  k0-1 'rh\CA/}D  
            bottom_y(i,j,2)=bottom_x(i,j,3)*(cos(psi)*cos(alpha)-cos(theta)*sin(psi)*sin(alpha));   %hy  k0-1 Q6|@N~UeZ  
             [[$Mh_MD  
            temp1(1)=i*sin(theta)*cos(psi)+j*sin(theta)*sin(psi)+(kete+1)*cos(theta); =ty2_6&>  
            temp1(2)=floor(temp1(1));                                %上边面入射  磁场 _;VYFs  
            temp1(3)=temp1(1)-temp1(2); sOC| B  
            temp1(4)=temp1(2)+1; oo'iwq-\  
            top_x(i,j,3)=(1-temp1(3))*hin(temp1(2))+temp1(3)*hin(temp1(4)); a^x  0 l  
            top_x(i,j,2)=top_x(i,j,3)*(-sin(psi)*cos(alpha)-cos(theta)*cos(psi)*sin(alpha));      %hx  kc+1 O E]~@eU  
            top_y(i,j,2)=top_x(i,j,3)*(cos(psi)*cos(alpha)-cos(theta)*sin(psi)*sin(alpha));       %hy  kc+1 YOlH*cZtg  
        end klo^K9!  
    end iI}nW  
     ^Wk0*.wg  
    for j=1:je R1~7F{FW  
        for k=1:ke BMF3XcH~G  
            temp1(1)=(iets-1)*sin(theta)*cos(psi)+j*sin(theta)*sin(psi)+k*cos(theta); X)y*#U  
            temp1(2)=floor(temp1(1));                                %后边面入射  磁场 lkyJ;}_**  
            temp1(3)=temp1(1)-temp1(2); "|\94  
            temp1(4)=temp1(2)+1; 3} l;  
            back_y(3,j,k)=(1-temp1(3))*hin(temp1(2))+temp1(3)*hin(temp1(4));
?Cc$]  
            back_y(2,j,k)=back_y(3,j,k)*(cos(psi)*cos(alpha)-cos(theta)*sin(psi)*sin(alpha));  %hy  i0-1 Upu%.[7  
            back_z(2,j,k)=back_y(3,j,k)*sin(theta)*sin(alpha);                                 %hz  i0-1 wV?[3bEhM  
             Hj1k-Bs&'w  
            temp1(1)=(iete+1)*sin(theta)*cos(psi)+j*sin(theta)*sin(psi)+k*cos(theta); 2t.fD@  
            temp1(2)=floor(temp1(1));                                %前边面入射  磁场 F-i&M1\_  
            temp1(3)=temp1(1)-temp1(2); L% zuI& q  
            temp1(4)=temp1(2)+1; " _mmR
M  
            front_y(3,j,k)=(1-temp1(3))*hin(temp1(2))+temp1(3)*hin(temp1(4)); -/1

d&  
            front_y(2,j,k)=front_y(3,j,k)*(cos(psi)*cos(alpha)-cos(theta)*sin(psi)*sin(alpha));   %hy  ia+1 <%(f9j  
            front_z(2,j,k)=front_y(3,j,k)*sin(theta)*sin(alpha);                                  %hz  ia+1 *eMLbU7  
        end iAg}pwU  
    end ?SB5b,  
    %*************************************************** %va[jJ  
    % 更新电场        (个别地方吸收边界无需特殊考虑，所以在设置吸收边界时，循环范围还需调整) ,qYf#fU#7  
    %*************************************************** Zq9>VqGe  
    ex(1:ib,2:je-1,2:ke-1)=ca*ex(1:ib,2:je-1,2:ke-1)+cb*(hz(1:ib,2:je-1,2:ke-1)-hz(1:ib,1:je-2,2:ke-1)-... br10ptEx  
        hy(1:ib,2:je-1,2:ke-1)+hy(1:ib,2:je-1,1:ke-2))/dx;       4bWfx_0W  
    ey(2:ie-1,1:jb,2:ke-1)=ca*ey(2:ie-1,1:jb,2:ke-1)+cb*(hx(2:ie-1,1:jb,2:ke-1)-hx(2:ie-1,1:jb,1:ke-2)-... FmR\`yY_,  
        hz(2:ie-1,1:jb,2:ke-1)+hz(1:ie-2,1:jb,2:ke-1))/dx; &9k"9  
    ez(2:ie-1,2:je-1,1:kb)=ca*ez(2:ie-1,2:je-1,1:kb)+cb*(hy(2:ie-1,2:je-1,1:kb)-hx(1:ie-2,2:je-1,1:kb)-... cx<h_  
        hx(2:ie-1,2:je-1,1:kb)+hx(2:ie-1,1:je-2,1:kb))/dx;   Us*
Vn  
    %*************************************************** jw/wcP  
    % 电场连接边界条件   边界面 ^=3 ^HQ'Zm  
    %*************************************************** [10$a(g\x  
    for j=jets-1:jete-1 OfW%&LAMQ  
        for k=kets-1:kete-1 $F<%Jl7_Z  
            ey(iets,j,k)=ey(iets,j,k)+back_z(2,j,k)*dt/epsz/dx; p+!f(H  
            ez(iets,j,k)=ez(iets,j,k)-back_y(2,j,k)*dt/epsz/dx; E=3#TBd  
            ey(iete,j,k)=ey(iete,j,k)-front_z(2,j,k)*dt/epsz/dx; :E}6S  
            ez(iete,j,k)=ez(iete,j,k)+front_y(2,j,k)*dt/epsz/dx; C?bXrG\  
        end 0;'j!`l9  
    end BR%{bY^ 5p  
    for i=iets-1:iete-1 La@\q[U{@  
        for k=kets-1:kete-1 _%XbxP6rH  
            ez(i,jets,k)=ez(i,jets,k)+left_x(i,2,k)*dt/epsz/dx; g+VRT,r  
            ex(i,jets,k)=ex(i,jets,k)-left_z(i,2,k)*dt/epsz/dx; ;k-g
_{M  
            ez(i,jete,k)=ez(i,jete,k)-right_x(i,2,k)*dt/epsz/dx; OrzM hQaf  
            ex(i,jete,k)=ex(i,jete,k)+right_z(i,2,k)*dt/epsz/dx; kK08W3@&t  
        end qNhH%tYQ
 
    end x!Y(Y=i>  
    for j=jets-1:jete-1 JHHb|  
        for k=kets-1:kete-1 9j9YQ2  
            ex(i,j,kets)=ex(i,j,kets)+bottom_y(i,j,2)*dt/epsz/dx; 9$8X>T^  
            ey(i,j,kets)=ey(i,j,kets)-bottom_x(i,j,2)*dt/epsz/dx; {P,>Q4N  
            ex(i,j,kete)=ex(i,j,kete)-top_y(i,j,2)*dt/epsz/dx; pFG]IM7o/u  
            ey(i,j,kete)=ey(i,j,kete)+top_x(i,j,2)*dt/epsz/dx; tvv[$b&  
        end 1fmSk$ y.9  
    end 3{I=.mUUm  
    %*************************************************** L)@`58Eil  
    % 电场连接边界条件  棱边 oM-b96  
    %***************************************************     lrq>TJEcx  
    ex(iets:iete,jets,kets)=ex(iets:iete,jets,kets)-left_z(iets:iete,2,kets)*dt/epsz/dx+... 3{6ps : w  
        bottom_y(iets:iete,jets,2)*dt/epsz/dx; <d3PDO@w/  
    ex(iets:iete,jete,kets)=ex(iets:iete,jete,kets)+right_z(iets:iete,2,kets)*dt/epsz/dx+... &%@/Dwr  
        bottom_y(iets:iete,jete,2)*dt/epsz/dx; Bi %Z2/  
    ex(iets:iete,jete,kete)=ex(iets:iete,jete,kete)+right_z(iets:iete,2,kete)*dt/epsz/dx-... 1*TXDo_T  
        top_y(iets:iete,jete,2)*dt/epsz/dx; 9T?~$XlX  
    ex(iets:iete,jets,kete)=ex(iets:iete,jets,kete)-left_z(iets:iete,2,kete)*dt/epsz/dx-... JvT%R`i  
        top_y(iets:iete,jets,2)*dt/epsz/dx; 6oPUYn-  
     ?=TL2"L  
    ey(iets,jets:jete,kets)=ey(iets,jets:jete,kets)-bottom_x(iets,jets:jete,2)*dt/epsz/dx+... '3IkPy1Uz  
        back_z(2,jets:jete,kets)*dt/epsz/dx; I=Gr^\x=  
    ey(iets,jets:jete,kete)=ey(iets,jets:jete,kete)+top_x(iets,jets:jete,2)*dt/epsz/dx+... PV5-^Y"v  
        back_z(2,jets:jete,kete)*dt/epsz/dx; &IIJKn|_  
    ey(iete,jets:jete,kete)=ey(iete,jets:jete,kete)+top_x(iete,jets:jete,2)*dt/epsz/dx-... 3}v0{c  
        front_z(2,jets:jete,kete)*dt/epsz/dx; uv?8V@x2  
    ey(iete,jets:jete,kets)=ey(iete,jets:jete,kets)-bottom_x(iete,jets:jete,2)*dt/epsz/dx-... xzuPie\  
        front_z(2,jets:jete,kets)*dt/epsz/dx;  >cC Gx  
     Ka[Sm|-q  
    ez(iets,jets,kets:kete)=ez(iets,jets,kets:kete)-back_y(2,jets,kets:kete)*dt/epsz/dx+... >|y>e{P  
        left_x(iets,2,kets:kete)*dt/epsz/dx; ))8Emk^Q{  
    ez(iete,jets,kets:kete)=ez(iete,jets,kets:kete)+front_y(2,jets,kets:kete)*dt/epsz/dx+... v#{G8'+%  
        left_x(iete,2,kets:kete)*dt/epsz/dx; [P(rY  
    ez(iete,jete,kets:kete)=ez(iete,jete,kets:kete)+front_y(2,jete,kets:kete)*dt/epsz/dx-... y^5T/M  
        right_x(iete,2,kets:kete)*dt/epsz/dx; gNG0k$nP  
    ez(iets,jete,kets:kete)=ez(iets,jete,kets:kete)-back_y(2,jete,kets:kete)*dt/epsz/dx-... IS3e|o*]MP  
        right_x(iets,2,kets:kete)*dt/epsz/dx; oUnq"]  
    %**************************************************** ~5x4?2  
    % 电场吸收边界    边界面 "6 uTo0  
    %**************************************************** }(8D!XgWa  
     $1:}(nO,  
     U&tfl/  
     m T>b;  
    %**************************************************** @&/s~3  
    % 电场吸收边界    棱边 f$~ _FX  
    %**************************************************** 7*R{u*/e  
     o6`Y7,]  
     J r=REa0  
     u<BHf@AI  
    %*************************************************** UwvGr h  
    % 更新磁场           (个别地方吸收边界无需特殊考虑，所以在设置吸收边界时，循环范围还需调整) [i2A{(x  
    %*************************************************** 0[_O+u  
    hx(2:ie-1,1:jb,1:kb)=cp*hx(2:ie-1,1:jb,1:kb)-cq*(ez(2:ie-1,2:jb+1,1:kb)-ez(2:ie-1,1:jb,1:kb)-... 1jR=h7^=  
        ey(2:ie-1,1:jb,2:kb+1)+ey(2:ie-1,1:jb,1:kb))/dx; "H wVK  
    hy(1:jb,2:je-1,1:kb)=cp*hy(1:jb,2:je-1,1:kb)-cq*(ex(1:jb,2:je-1,2:kb+1)-ex(1:jb,2:je-1,1:kb)-... GLbc/qs  
        ez(2:jb+1,2:je-1,1:kb)+ez(1:jb,2:je-1,1:kb))/dx; ,+swH;=7#r  
    hz(1:ib,1:jb,2:ke-1)=cp*hz(1:ib,1:jb,2:ke-1)-cq*(ey(2:ib+1,1:jb,2:ke-1)-ey(1:ib,1:jb,2:ke-1)-... |?4~T:  
        ex(1:ib,2:jb+1,2:ke-1)+ex(1:ib,1:jb,2:ke-1))/dx;   ~xsb5M5  
    %*************************************************** 8#NIs@DJ  
    % 磁场连接边界条件   边界面 !y[}|  
    %*************************************************** t?pIE cl  
    for j=jets:jete we6+2  
        for k=kets:kete [,.[gWA  
            hy(iets,j,k)=hy(iets,j,k)-back_z(1,j,k)*dt/murz/dx; }[XB]Xf  
            hz(iets,j,k)=hz(iets,j,k)+back_y(1,j,k)*dt/murz/dx; 5P5A,K  
            hy(iete,j,k)=hy(iete,j,k)+front_z(1,j,k)*dt/murz/dx; FN{H\W1cf  
            hz(iete,j,k)=hz(iete,j,k)-front_y(1,j,k)*dt/murz/dx; [a#?}((  
        end ?uNTUU,  
    end 7C,<iY  
    for i=iets:iete r{;VTQ  
        for k=kets:kete a>;3 j  
            hz(i,jets,k)=hz(i,jets,k)-left_x(i,1,k)*dt/murz/dx; ]{q-Y<{"  
            hx(i,jets,k)=hx(i,jets,k)+left_z(i,1,k)*dt/murz/dx; A52LH,  
            hz(i,jete,k)=hz(i,jete,k)+right_x(i,1,k)*dt/murz/dx; c+)36/; X  
            hx(i,jete,k)=hx(i,jete,k)-right_z(i,1,k)*dt/murz/dx; /h.:br?M#P  
        end "t3uW6&  
    end
<-Kb@V3  
    for i=iets:iete y\r^\ S9%  
        for j=jets:jete o(v"?Y6  
            hx(i,j,kets)=hx(i,j,kets)-bottom_y(i,j,1)*dt/murz/dx; ^+.+IcH  
            hy(i,j,kets)=hy(i,j,kets)+bottom_x(i,j,1)*dt/murz/dx; yoq\9* ?u^  
            hx(i,j,kete)=hx(i,j,kete)+top_y(i,j,1)*dt/murz/dx; =rf)yp-D  
            hy(i,j,kete)=hy(i,j,kete)-top_x(i,j,1)*dt/murz/dx; (1saof*p%  
        end F)[XIY&2/  
    end llTQ\7zP  
    %**************************************************** Zi2Eu4p l{  
    % 磁场吸收边界    边界面 1[fkXO{  
    %**************************************************** ,6y-.m7>  
     [="g|/M)  
     Y&1!Z*OL;  
     E-BOIy,  
     3[00-~&U  
     wI1M0@}PV  
    %**************************************************** e{qp!N1!  
    % Visualize fields /h}wM6pg  
    %**************************************************** [N0/">c  
    if mod(n,10)==0; Y%n{`9=  
        timestep=int2str(n); xr7-[)3Q$  
        tview(:,:)=ez(:,:,ke/2); t wa(M?  
        sview(:,:)=ez(:,je/2,:); ~#IWM+I  
u`L!za7fi  
        subplot(2,1,1),pcolor(tview'); i4{ /  
        shading flat; &CxyP_  
        caxis([-1.0 1.0]); z<3}TD  
        colorbar; {Kq*5Aq8  
        axis image; [FeJ8P>z  
        title(['Ez(i,j,k=50),timestep=',timestep]); 
pUCEYR  
        xlabel('i coordinate'); =Ov;'MC  
        ylabel('j coordinate'); #2ZrdD"5kQ  
LXTipWKz  
        subplot(2,1,2),pcolor(sview'); n;S0fg  
        shading flat; 'AAF/9  
        caxis([-1.0 1.0]); !DY2{Wb  
        colorbar; JWUv H  
        axis image; YKs^%GO+  
        title(['Ez(i,j=50,k),timestep=',timestep]); WNF=NNO-R  
        xlabel('i coordinate'); uJ)\P  
        ylabel('k coordinate'); wHo#%Y,Nmi  
         f//j{P[  
        nn=n/10; ?h>%Ix  
        M(:,nn)=getframe(gcf,rect); z$8e6*  
    end 9:}RlL+cOk  
end 5W:Gl?$S}  
movie(gcf,M,0,10,rect);