[求助]三维时谐场——介质球的RCS计算

       下面是我整的三维时谐场——用FDTD计算介质球的RCS的程序，但是在近远场外推的时候出现了问题，怎么都跑不出图来， ~nb(e$?N  
希望大家能帮我找出错误，谢谢！ o+.ySSBl+  
       IE=81; r`\@Fv,&#  
JE=81; &;~?\>?I  
KE=81; Zl{9G?abCT  
ia=9; N.0g%0A.D  
ja=9; !l]_c5  
ka=9; bDm7$ (  
ic=(IE+1)/2; OL%}C*Zq  
jc=(JE+1)/2; MiR$N  
kc=(KE+1)/2; ~FQHT?DAo  
ib=IE-ia; P)#h4|xZ  
jb=JE-ja; ]s!id[j
 
kb=KE-ka; dS\!tdHP-Q  
epsz=8.8e-12; Y`#6MhFT7  
muz=4*pi*1e-7; -+M360  
lamda=0.032; =XK}eQ_d  
freq=lamda/3e+8; =
z]rZSq*o  
%freq=30e+7; jgS
3#  
ddx=lamda/42; KMK8jJ  
dt=ddx/6e+8; ve:Oe{Ie{  
arg=2*pi*freq*dt; y/:%S2za>  
dx=zeros(IE,JE,KE); wph8ln"C-  
dy=zeros(IE,JE,KE); 9 )!
}  
dz=zeros(IE,JE,KE); kDXQpe  
ex=zeros(IE,JE,KE); (> _Lb  
ey=zeros(IE,JE,KE); uE
<8L(*B  
ez=zeros(IE,JE,KE); |>[qC O  
hx=zeros(IE,JE,KE); <*WGvCh%w  
hy=zeros(IE,JE,KE); w/"vf3}(9  
hz=zeros(IE,JE,KE); {bR2S&=OmK  
ix=zeros(IE,JE,KE); IUDH"~f  
iy=zeros(IE,JE,KE); U{/fY/kq  
iz=zeros(IE,JE,KE); _RbM'_y+E  
gax=zeros(IE,JE,KE); SM%/pu;  
gay=zeros(IE,JE,KE); {}rnn$HQe  
gaz=zeros(IE,JE,KE); 6yO5{._M  
gbx=zeros(IE,JE,KE); G.^)5!By  
gby=zeros(IE,JE,KE); l($8HAJ  
gbz=zeros(IE,JE,KE); '7/
F]S0K  
END=200; >*opEI+  
e_inc=zeros(1,END); 0kNKt(_  
h_inc=zeros(1,END); TLp2a<Iy  
e_temp=zeros(1,4); =mF"D:s*  
idxl=zeros(ia,JE,KE); ^<;CIXo  
idxh=zeros(ia,JE,KE); Tl'wA^~H  
ihxl=zeros(ia,JE,KE); '=%`;?j  
ihxh=zeros(ia,JE,KE); j*[P\Cm  
idyl=zeros(IE,ja,KE); *^Ges;5$"  
idyh=zeros(IE,ja,KE); /-i m g^^  
ihyl=zeros(IE,ja,KE); Qe\vx1GRLH  
ihyh=zeros(IE,ja,KE); W$2\GPJt  
idzl=zeros(IE,JE,ka); MnZljB  
idzh=zeros(IE,JE,ka); 1G.?Y3DC<  
ihzl=zeros(IE,JE,ka); ^1vKhO+p$  
ihzh=zeros(IE,JE,ka); dvqg H  
gi1=zeros(1,IE);gi2=ones(1,IE);gi3=ones(1,IE); :CK`v6 Qs  
gj1=zeros(1,JE);gj2=ones(1,JE);gj3=ones(1,JE); Shn=Q  
gk1=zeros(1,KE);gk2=ones(1,KE);gk3=ones(1,KE); #ACT&J  
fi1=zeros(1,IE);fi2=ones(1,IE);fi3=ones(1,IE); (_8.gS[  
fj1=zeros(1,JE);fj2=ones(1,JE);fj3=ones(1,JE); dP(.l}O  
fk1=zeros(1,KE);fk2=ones(1,KE);fk3=ones(1,KE); U\B9Ab  
sita=0*pi/180;    %pointing angle (QL:7  
fai=0*pi/180;     %azimuth angle CLk,]kA'r  
arfa=0*pi/180;    %polarizing angle RgUQ
:  
sp1=sin(sita);cp1=cos(sita); }( CYok  
sp2=sin(fai);cp2=cos(fai); &P>& T  
sp3=sin(arfa);cp3=cos(arfa); 8)=(
eI$  
ex_left=zeros(IE,KE);ex_right=zeros(IE,KE); c"-X:m"  
ex_up=zeros(IE,JE);ex_down=zeros(IE,JE); |JF,n~n  
ey_front=zeros(JE,KE);ey_back=zeros(JE,KE); 6W{Nw<  
ey_up=zeros(IE,JE);ey_down=zeros(IE,JE); W]DGt|JP  
ez_left=zeros(IE,KE);ez_right=zeros(IE,KE); [@5cYeW3.  
ez_front=zeros(JE,KE);ez_back=zeros(JE,KE); 8Iu6r}k?~`  
hx_left=zeros(IE,KE);hx_right=zeros(IE,KE); *~shvtq  
hx_up=zeros(IE,JE);hx_down=zeros(IE,JE); dKCl#~LAI'  
hy_front=zeros(JE,KE);hy_back=zeros(JE,KE);
4x(m.u@  
hy_up=zeros(IE,JE);hy_down=zeros(IE,JE); ?jUgDwc(w  
hz_left=zeros(IE,KE);hz_right=zeros(IE,KE); h@\-]zN{  
hz_front=zeros(JE,KE);hz_back=zeros(JE,KE); D!! B4zt  
npml=4;        %pml层的厚度 t!SxJB e  
xxn=1:-1/npml:1/npml; TWGn:mi  
xn=0.3333*xxn.^3; yn$1n
t4  
fi1(1:npml)=xn; nW"O+s3  
fi1(IE:-1:IE-npml+1)=xn; Xwdcy J!  
gi2(1:npml)=1./(1+xn); >l><d!
hw  
gi2(IE:-1:IE-npml+1)=1./(1+xn); T4Vp0i  
gi3(1:npml)=(1-xn)./(1+xn); o$l8"Uv  
gi3(IE:-1:IE-npml+1)=(1-xn)./(1+xn); /;d 5p  
fj1(1:npml)=xn; |9\i+)C  
fj1(JE:-1:JE-npml+1)=xn; m$^5{qpg  
gj2(1:npml)=1./(1+xn); /<n7iIK)  
gj2(JE:-1:JE-npml+1)=1./(1+xn); F-rhxJd  
gj3(1:npml)=(1-xn)./(1+xn); [j5+PV  
gj3(JE:-1:JE-npml+1)=(1-xn)./(1+xn); }4!}vkVx  
fk1(1:npml)=xn; xvl
{o  
fk1(KE:-1:KE-npml+1)=xn; fs;\_E
[)  
gk2(1:npml)=1./(1+xn); S|)atJJ0G"  
gk2(KE:-1:KE-npml+1)=1./(1+xn); fdX|t"oz  
gk3(1:npml)=(1-xn)./(1+xn); $Wj{B@k  
gk3(KE:-1:KE-npml+1)=(1-xn)./(1+xn); gC(S(osF  
xxn=(1-0.5/npml):-1/npml:0.5/npml; -dO8Uis$  
xn=0.3333*xxn.^3; b@8z+,_  
gi1(1:npml)=xn; MD;Z UAX<  
gi1(IE-1:-1:IE-npml)=xn; `ovMfL.u  
fi2(1:npml)=1./(1+xn); "qF/7`e[  
fi2(IE-1:-1:IE-npml)=1./(1+xn); du$M
  
fi3(1:npml)=(1-xn)./(1+xn); GiX3c^V"1  
fi3(IE-1:-1:IE-npml)=(1-xn)./(1+xn); } % Ie  
gj1(1:npml)=xn; FXDB> }8  
gj1(JE-1:-1:JE-npml)=xn; }xt^}:D  
fj2(1:npml)=1./(1+xn); )1B?<4  
fj2(JE-1:-1:JE-npml)=1./(1+xn); ym%` l!
 
fj3(1:npml)=(1-xn)./(1+xn); 9xz@2b@  
fj3(JE-1:-1:JE-npml)=(1-xn)./(1+xn); ^pd7nr~Y  
gk1(1:npml)=xn; 4gC(zJ  
gk1(KE-1:-1:KE-npml)=xn; _q#pEv  
fk2(1:npml)=1./(1+xn); @@U  
fk2(KE-1:-1:KE-npml)=1./(1+xn); ]s0wJD=  
fk3(1:npml)=(1-xn)./(1+xn); poW%Fzj  
fk3(KE-1:-1:KE-npml)=(1-xn)./(1+xn); ><D2of|  
WL=42;
=E]tEi  
dT=WL/2; 8Cp@k=  
omiga=pi/WL; A :KZyd"Z  
Z=sqrt(muz/epsz); RHn3\N  
odist=2; T=u"y;&L  
iao=ia+1-odist; J/K~8sc  
jao=ja+1-odist; qQ^CSn98J  
kao=ka+1-odist; @@K/0:],  
ibo=ib+odist; gAorb\
iJ  
jbo=jb+odist; ru2M"]T  
kbo=kb+odist; ,LxZbo!  
exi_oup=zeros(IE,JE);exr_oup=zeros(IE,JE);ex_oup=zeros(IE,JE); FCEmg0qdjD  
exi_odown=zeros(IE,JE);exr_odown=zeros(IE,JE);ex_odown=zeros(IE,JE); !KOa'Ic$V  
eyi_oup=zeros(IE,JE);eyr_oup=zeros(IE,JE);ey_oup=zeros(IE,JE); |[iO./zP  
eyi_odown=zeros(IE,JE);eyr_odown=zeros(IE,JE);ey_odown=zeros(IE,JE); aY?VP?BL  
exi_oleft=zeros(IE,KE);exr_oleft=zeros(IE,KE);ex_oleft=zeros(IE,KE); =cS5f#0  
ezi_oleft=zeros(IE,KE);ezr_oleft=zeros(IE,KE);ez_oleft=zeros(IE,KE); rZpsC}C'  
exi_oright=zeros(IE,KE);exr_oright=zeros(IE,KE);ex_oright=zeros(IE,KE); i&lW&]  
ezi_oright=zeros(IE,KE);ezr_oright=zeros(IE,KE);ez_oright=zeros(IE,KE); +@!\3a4!  
eyi_ofront=zeros(JE,KE);eyr_ofront=zeros(JE,KE);ey_ofront=zeros(JE,KE); i"iy 0?  
ezi_ofront=zeros(JE,KE);ezr_ofront=zeros(JE,KE);ez_ofront=zeros(JE,KE); frPQi{u$  
eyi_oback=zeros(JE,KE);eyr_oback=zeros(JE,KE);ey_oback=zeros(JE,KE); { "Cu)AFy  
ezi_oback=zeros(JE,KE);ezr_oback=zeros(JE,KE);ez_oback=zeros(JE,KE); 5j.@)XXe  
hxi_oup=zeros(IE,JE);hxr_oup=zeros(IE,JE);hx_oup=zeros(IE,JE); (nq""kO6'  
hxi_odown=zeros(IE,JE);hxr_odown=zeros(IE,JE);hx_odown=zeros(IE,JE); C=r`\W  
hyi_oup=zeros(IE,JE);hyr_oup=zeros(IE,JE);hy_oup=zeros(IE,JE); [zN*P$U]  
hyi_odown=zeros(IE,JE);hyr_odown=zeros(IE,JE);hy_odown=zeros(IE,JE); 0(Y,Q(JTo&  
hxi_oleft=zeros(IE,KE);hxr_oleft=zeros(IE,KE);hx_oleft=zeros(IE,KE); Zl&ED{k<  
hzi_oleft=zeros(IE,KE);hzr_oleft=zeros(IE,KE);hz_oleft=zeros(IE,KE); ,~38IIS>_  
hxi_oright=zeros(IE,KE);hxr_oright=zeros(IE,KE);hx_oright=zeros(IE,KE); 2IW!EUR  
hzi_oright=zeros(IE,KE);hzr_oright=zeros(IE,KE);hz_oright=zeros(IE,KE); 9y&;6V.'  
hyi_ofront=zeros(JE,KE);hyr_ofront=zeros(JE,KE);hy_ofront=zeros(JE,KE); }JST(d&  
hzi_ofront=zeros(JE,KE);hzr_ofront=zeros(JE,KE);hz_ofront=zeros(JE,KE); CKZEX*mPC  
hyi_oback=zeros(JE,KE);hyr_oback=zeros(JE,KE);hy_oback=zeros(JE,KE); T^#d;A  
hzi_oback=zeros(JE,KE);hzr_oback=zeros(JE,KE);hz_oback=zeros(JE,KE); F*#!hWtb  
%%%介质球的设置 e#k rr  
epsilon=[1 4]; Gr&e]M[l  
%epsilon=[1 81]; aW dI  
sigma=[0 0]; "]G'^  
%sigma=[0 5]; :3R3>o6m  
%sigma=[0 0]; cq?,v?m  
radius=20; `ql8y'  
for i=ia+1:ib 5EQ)pH+  
    for j=ja+1:jb *&IvEu  
        for k=ka+1:kb BNjMq  
            eps=epsilon(1); ]!^wB 3j  
            cond=sigma(1); qS! Lt3+  
            xdist=(ic-i-0.5); 5<`83;R9  
            ydist=jc-j; ktynI
N  
            zdist=kc-k; iR9duP+  
            dist=sqrt(xdist^2+ydist^2+zdist^2); +U J~/XV  
            if dist<=radius Hs8JJGXWB  
                eps=epsilon(2); $rk=#;6]v;  
                cond=sigma(2); Q.eD:
@%iE  
            end ]?T^tJ  
            gax(i,j,k)=1/(eps+cond*dt/epsz); m%})H"5  
            gbx(i,j,k)=cond*dt/epsz; Tim/7*vx  
        end PPV T2;9  
    end PR!0=E*}  
end Ig?9"{9p  
for i=ia+1:ib cy!P!t,@  
    for j=ja+1:jb 9]
ZfSn)  
        for k=ka+1:kb b49h @G  
            eps=epsilon(1); 6p&2A  
            cond=sigma(1); g
|| q 3  
            xdist=ic-i; b|k^   
            ydist=jc-j-0.5; CfS;F  
            zdist=kc-k; +RM!j9Rq  
            dist=sqrt(xdist^2+ydist^2+zdist^2); +924_,zF  
            if dist<=radius 4@\$k+v  
                eps=epsilon(2); 0[d*Z  
                cond=sigma(2); /^jl||'H,:  
            end ._j?1Fw`  
            gay(i,j,k)=1/(eps+cond*dt/epsz); 1>\V>g
9  
            gby(i,j,k)=cond*dt/epsz; QA^FP8!j  
        end V\L%*6O  
    end Q6URaw#Yt`  
end GY@:[u.&  
for i=ia+1:ib G?LC!9MB  
    for j=ja+1:jb #+_=(J  
        for k=ka+1:kb j~.tyxOq#  
            eps=epsilon(1); 'J0I$-QYk  
            cond=sigma(1); wsQuJrG  
            xdist=ic-i; !PeSnO  
            ydist=jc-j; 6A=k;do  
            zdist=kc-k-0.5; R#y"SxD()  
            dist=sqrt(xdist^2+ydist^2+zdist^2); JQej$=*  
            if dist<=radius
P"}"q ![  
                eps=epsilon(2); PU%f`)  
                cond=sigma(2); >0iCQKq  
            end tKuJ &I~  
            gaz(i,j,k)=1/(eps+cond*dt/epsz); ~&<vAgy,  
            gbz(i,j,k)=cond*dt/epsz; Zw{?^6;cS  
        end 3:#6/@wQ  
    end S?JGg.)  
end f>Ua7!b  
real_in=0; V'dw=W17V  
imag_in=0; x%ccNP0  
%t0=40; G `3{Q7k  
%spread=10; -\B*reC  
T=1000; h:G>w`X  
%N=0; 6!itr"  
%isource=floor(sqrt((ib-ia)^2+(jb-ja)^2+(kb-ka)^2)); "WtYqXyd  
isource=4; T+RC#&>  
is=ja+1; '`<Fys&:  
for N=1:T dP_bFUzg  
    Yd~J(  
    for j=2:END `bV&n!Y_  
        e_inc(j)=e_inc(j)+0.5*(h_inc(j-1)-h_inc(j)); EBL-+%J8  
    end 3'i(wI~
<[  
    NySa%7@CD  
    %real_in=real_in+cos(arg*T)*e_inc(ja-1); \JR^uJ{Y  
    %imag_in=imag_in-sin(arg*T)*e_inc(ja-1); [742s]j  
    x3U>5F@  
    %pulse=exp(-0.5*((t0-T)/spread)^2); +03/A`PKrB  
    pulse=sin(arg*N); >/ A'G  
    e_inc(isource)=pulse; =w`uZ;l$Q  
    R:Ih#2R  
    e_inc(1)=e_temp(1); #l- 0$  
    e_temp(1)=e_temp(2); T fIOS]  
    e_temp(2)=e_inc(2); za 7+xF  
    e_inc(END)=e_temp(3); [MbbL  
    e_temp(3)=e_temp(4); ="vg/@.>i  
    e_temp(4)=e_inc(END-1); <z#Fj`2{  
    KkpbZ7\@  
    %%%  Calculate Incident component on connected edge Eld[z{n"  
    %%%  up/down connect edge 88S:E7 $  
    for i=ia+1:ib i0+e3!QU  
        for j=ja+1:jb [kxOv7a  
            k=kb; a(?)r[=  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; yw7
(!1j=  
            w=dr-floor(dr); kc=Z6(=  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); ~G=E Q]a
  
            ex_up(i,j)=ein*(sp2*cp3-cp1*cp2*sp3); O8" t.W  
            3>MILEY^  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; 3< 6h~ek)  
            w=dr-floor(dr); 9v-Y*\!w.  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); Q}<QE:-&E  
            ey_up(i,j)=ein*(-cp2*cp3-cp1*sp2*sp3); uHmvHA~/c8  
            nsVLgTbx  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; H-u SdT  
            w=dr-floor(dr); r1vS~ 4Z  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); &&LB0vH
!J  
            hx_up(i,j)=hin*(sp2*sp3+cp1*cp2*cp3); {dRZ2U3  
            Oi^cs=}  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; 5cU:wc  
            w=dr-floor(dr); kn$_X
4^?  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); q
I~*G3  
            hy_up(i,j)=hin*(-cp2*sp3+cp1*sp2*cp3); -
Hw3rv3o  
        end ]yqE6Lf9  
    end } d8\ Jg  
    QZ;DZMP  
    for i=ia+1:ib P,-5af*;  
        for j=ja+1:jb ]e"!ZR?XJ  
            k=ka+1; &)||~  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; srO>l ;Vf/  
            w=dr-floor(dr); #Y   
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); \~Z%}$ =  
            ex_down(i,j)=ein*(sp2*cp3-cp1*cp2*sp3); >35w"a7S  
            ?)k]Vg.  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; <pHm=q/U  
            w=dr-floor(dr); eu_ZsseZ  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); @^-Y&N!b=  
            ey_down(i,j)=ein*(-cp2*cp3-cp1*sp2*sp3); f`/('}t  
            |%F4`gz8KP  
            k=ka; <`; {gX1  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; "cBqZzkk9j  
            w=dr-floor(dr); hp9LV2_5  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); <BZC5b6  
            hx_down(i,j)=hin*(sp2*sp3+cp1*cp2*cp3); Nz`v+sp  
            U{pg y#/  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; r`EjD}2d  
            w=dr-floor(dr); 39P55B/o%  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); U{[YCs fk  
            hy_down(i,j)=hin*(-cp2*sp3+cp1*sp2*cp3); u1#(~[.  
        end tAJ}36aG  
    end uG6.(A1LM  
    U2WHs3  
    %%%  left/right connect edge tQNrDp+  
    for i=ia+1:ib rrj.]^E_~  
        for k=ka+1:kb o'(BL:8s  
            j=ja+1; m9:ah<  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; \**j\m   
            w=dr-floor(dr); } -;)G~h/"  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); eQ8t.~5;-  
            ex_left(i,k)=ein*(sp2*cp3-cp1*cp2*sp3); S`FIb'J  
            &<uLr *+*  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; 5n}<V-yJ*m  
            w=dr-floor(dr); .n YlYY'   
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); 6XU p$Pd(  
            ez_left(i,k)=ein*(sp1*sp3); o}/|"(K  
            6G"UXNa,  
            j=ja; "qR, V9\  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; .RFH@''  
            w=dr-floor(dr); 5<v1v
&  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); +ls`;f  
            hx_left(i,k)=hin*(sp2*sp3+cp1*cp2*cp3); QQV8Vlv"  
            =G^'wwpv(  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; ku]?"{Xx  
            w=dr-floor(dr); `\\s%}vZ*T  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); *xsBFCRU  
            hz_left(i,k)=hin*(-sp1*cp3); ysIhUpd  
        end y'4Qt.1ukN  
    end "uIa
K
b  
       YMK ![ q-  
    for i=ia+1:ib '
=Lpch2J  
        for k=ka+1:kb Ow4(1eE_  
            j=jb; (y.
N-I,  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; {CBb^BP  
            w=dr-floor(dr); LOfw #+]d  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); _[Imwu}  
            ex_right(i,k)=ein*(sp2*cp3-cp1*cp2*sp3); _~\} fY  
            <n#X~}i)  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; `m<O!I"A  
            w=dr-floor(dr); jED.0,+K!  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); YkB@fTTS
 
            ez_right(i,k)=ein*(sp1*sp3); J-dB  
            v7./u4S|V  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; xt,Qn460;  
            w=dr-floor(dr); 1@KiP`DA  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); H~Vf;k>  
            hx_right(i,k)=hin*(sp2*sp3+cp1*cp2*cp3); 9.M'FCd~M  
            ug2W{D  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; pUqC88*j  
            w=dr-floor(dr); LAf#Rco4  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); ^&1O:G*"  
            hz_right(i,k)=hin*(-sp1*cp3); =VuSi(d;e{  
        end 9+N%Io?!  
    end I(pq3_9$  
    h L [eA  
    %%%  front/back connect edge q: FhuOP  
    for j=ja+1:jb miQ*enZi  
        for k=ka+1:kb -NN=(p!<  
            i=ib; &Q?@VNi  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; h`1<+1J9  
            w=dr-floor(dr); \b(&-=(  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); [F+W]Jk,  
            ey_front(j,k)=ein*(-cp2*cp3-cp1*sp2*sp3); r\ft{Z<P  
            xLoQ0rt 6  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; Nv36#^Z  
            w=dr-floor(dr); =ejU(1 g  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); =cjO]  
            ez_front(j,k)=ein*(sp1*sp3); pl&nr7\  
            /^nIOAeE  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; {P~rf&Ee  
            w=dr-floor(dr); IV. })8
 
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1);
qNj?Rwc  
            hy_front(j,k)=hin*(-cp2*sp3+cp1*sp2*cp3); H2R3I<j  
            #lvt4a"P"  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; a,RCK~GR  
            w=dr-floor(dr); z6E =%-`  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); 8ex;g^e  
            hz_front(j,k)=hin*(-sp1*cp3); l|gi2~ %Y  
        end zQY ,}a  
    end o$.#A]Flb  
    PJN9[Y{^3  
    for j=ja+1:jb n Ab~  
        for k=ka+1:kb %w65)BFQ  
            i=ia+1; #'s$6gT=  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; -\?-  
            w=dr-floor(dr); 8Zsaq1S  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); }BlyEcw'aN  
            ey_back(j,k)=ein*(-cp2*cp3-cp1*sp2*sp3); wX]$xZ!s  
            Ju47}t%HB  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; pPRX#3
 
            w=dr-floor(dr); ]}rNxT4<  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); { %X2K  
            ez_back(j,k)=ein*(sp1*sp3); zG='U  
            4DCh+|r  
            i=ia; diJpbR^JP  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; IXa~,a H71  
            w=dr-floor(dr); OmWEa  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); ~-7/9$ay5  
            hy_back(j,k)=hin*(-cp2*sp3+cp1*sp2*cp3); ?{f6su@rW  
            08nh y[  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; VR>!Ch  
            w=dr-floor(dr); uKk#V6
t#  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); bOr11?
 
            hz_back(j,k)=hin*(-sp1*cp3); St%x\[D  
        end "crR{OjE"  
    end LQPQ !):;  
    I &iyj99n  
    %%% Dx H;nzo3x  
    for j=2:JE mgx|5Otg  
        for k=2:KE MZh.Xo  
            for i=2:ia u5}:[4N%I  
                curl_h=hz(i,j,k)-hz(i,j-1,k)-hy(i,j,k)+hy(i,j,k-1); NzQvciJ@"  
                idxl(i,j,k)=idxl(i,j,k)+curl_h; f~mwDkf?L  
                dx(i,j,k)=gj3(j)*gk3(k)*dx(i,j,k)+gj2(j)*gk2(k)*0.5*(curl_h+gi1(i)*idxl(i,j,k)); jJiuq#;T3  
            end %;:![?M  
            for i=ia+1:ib M,H8ZO:R  
                curl_h=hz(i,j,k)-hz(i,j-1,k)-hy(i,j,k)+hy(i,j,k-1); Ljz)%y[s  
                dx(i,j,k)=gj3(j)*gk3(k)*dx(i,j,k)+gj2(j)*gk2(k)*0.5*curl_h; ]w6F%d  
            end PDD2ouv4  
            for i=ib+1:IE @9 S ::  
                ixh=i-ib; vm+3!s:u  
                curl_h=hz(i,j,k)-hz(i,j-1,k)-hy(i,j,k)+hy(i,j,k-1); (]'wQ4iQ  
                idxh(ixh,j,k)=idxh(ixh,j,k)+curl_h; L1RD`qXu.  
                dx(i,j,k)=gj3(j)*gk3(k)*dx(i,j,k)+gj2(j)*gk2(k)*0.5*(curl_h+gi1(i)*idxh(ixh,j,k)); |9S8sfw  
            end Q*#Lr4cm{  
        end )m7%cyfC  
    end Cu#n5SF*  
   j.Uy>ol  
    for i=ia+1:ib-1 vf3)T;X>  
        for k=ka+1:kb j13-?fQ&  
            dx(i,ja+1,k)=dx(i,ja+1,k)-0.5*hz_left(i,k); M8WjqTq  
            dx(i,jb,k)=dx(i,jb,k)+0.5*hz_right(i,k); #VX]tr
h,  
        end EX{%CPp7}  
    end ck]I?  
    a8T9=KY^  
    for i=ia+1:ib-1 #(614-r/  
        for j=ja+1:jb GqCBD-@4v.  
            dx(i,j,ka+1)=dx(i,j,ka+1)+0.5*hy_down(i,j); AJi+JO-  
            dx(i,j,kb)=dx(i,j,kb)-0.5*hy_up(i,j); D*-  
        end  ?pEPwc  
    end E6~VHQa2?  
    8X`DFeJ  
    %%% Dy 6Z#Nh@!+C  
    for i=2:IE 4utwcXL  
        for k=2:KE ^x O](,H  
            for j=2:ja o i'iZX  
                curl_h=hx(i,j,k)-hx(i,j,k-1)-hz(i,j,k)+hz(i-1,j,k); ;?HP/dZLz  
                idyl(i,j,k)=idyl(i,j,k)+curl_h; 4Y59^  
                dy(i,j,k)=gi3(i)*gk3(k)*dy(i,j,k)+gi2(i)*gk2(k)*0.5*(curl_h+gj1(j)*idyl(i,j,k)); KU$,{Sn6@  
            end .c]>*/(+  
            for j=ja+1:jb a+LK~mC*  
                curl_h=hx(i,j,k)-hx(i,j,k-1)-hz(i,j,k)+hz(i-1,j,k); h623)C;  
                dy(i,j,k)=gi3(i)*gk3(k)*dy(i,j,k)+gi2(i)*gk2(k)*0.5*curl_h; L-?ty@-i  
            end `"CA$Se8  
            for j=jb+1:JE H1U$ApD  
            jyh=j-jb; %l&oRBC  
            curl_h=hx(i,j,k)-hx(i,j,k-1)-hz(i,j,k)+hz(i-1,j,k); Ne<S_u2nT  
            idyh(i,jyh,k)=idyh(i,jyh,k)+curl_h; y$7Ys:R~  
            dy(i,j,k)=gi3(i)*gk3(k)*dy(i,j,k)+gi2(i)*gk2(k)*0.5*(curl_h+gj1(j)*idyh(i,jyh,k)); >A{Dpsi\  
            end cL#-vW<s3  
        end > .NLmzUX  
    end bxh-#x &  
    }F{s\qUt  
    for i=ia+1:ib O)&W0`VY  
        for j=ja+1:jb-1 _$$.5?4  
            dy(i,j,ka+1)=dy(i,j,ka+1)-0.5*hx_down(i,j); VCc=dME  
            dy(i,j,kb)=dy(i,j,kb)+0.5*hx_up(i,j); u=nd7:bv  
        end E}9wzPs  
    end !~C%0{9+u@  
    ( xooU 8d  
    for j=ja+1:jb-1 )e0kr46  
        for k=ka+1:kb k vZw4Pk  
            dy(ia+1,j,k)=dy(ia+1,j,k)+0.5*hz_back(j,k); P.Bwfa  
            dy(ib,j,k)=dy(ib,j,k)-0.5*hz_front(j,k); n32"cFPpT  
        end #:BkDidt2v  
    end npzp/mcIe)  
   u4FD}nV  
    %%%  Dz W6>t!1oO+  
    for i=2:IE JqO1 a?H  
        for j=2:JE
tm5{h{AM  
            for k=2:ka c?CfM>  
                curl_h=hy(i,j,k)-hy(i-1,j,k)-hx(i,j,k)+hx(i,j-1,k); V. i
{IW  
                idzl(i,j,k)=idzl(i,j,k)+curl_h; ^dLu#,;  
                dz(i,j,k)=gi3(i)*gj3(j)*dz(i,j,k)+gi2(i)*gj2(j)*0.5*(curl_h+gk1(k)*idzl(i,j,k)); AmIW
$(Ce  
            end %]7 6u7b/  
            for k=ka+1:kb OcV,pJ  
                curl_h=hy(i,j,k)-hy(i-1,j,k)-hx(i,j,k)+hx(i,j-1,k); ?0:]%t18  
                dz(i,j,k)=gi3(i)*gj3(j)*dz(i,j,k)+gi2(i)*gj2(j)*0.5*curl_h; u_NLgM7*  
            end ~3M4F^  
            for k=kb+1:KE 1LS1 ZY  
                kzh=k-kb; 4Lg ,J9  
                curl_h=hy(i,j,k)-hy(i-1,j,k)-hx(i,j,k)+hx(i,j-1,k); QBGm)h?=  
                idzh(i,j,kzh)=idzh(i,j,kzh)+curl_h; m
mJnE  
                dz(i,j,k)=gi3(i)*gj3(j)*dz(i,j,k)+gi2(i)*gj2(j)*0.5*(curl_h+gk1(k)*idzh(i,j,kzh)); M *w{PjU  
            end ?`e@ o?  
        end zB0*KgAn{  
    end s.7=!JQ#]p  
    |Io:D:  
    for i=ia+1:ib N,lr~6)  
        for k=ka+1:kb-1 =| T^)J  
            dz(i,ja+1,k)=dz(i,ja+1,k)+0.5*hx_left(i,k); :y7K3:d3  
            dz(i,jb,k)=dz(i,jb,k)-0.5*hx_right(i,k); [c=P)t7 V  
        end ,u>LAo0  
    end (rhlK} C  
    'i$._Tx  
    for j=ja+1:jb DqWy@7 a  
        for k=ka+1:kb-1 2F*>&n&Db7  
            dz(ia+1,j,k)=dz(ia+1,j,k)-0.5*hy_back(j,k); |oU I2<"  
            dz(ib,j,k)=dz(ib,j,k)+0.5*hy_front(j,k);  t* Ct*  
        end :z P:4NW  
    end vrb@::sy0T  
    rzHBop-8  
    %%%  Calculate E from D v3cMPN  
    for i=2:IE-1 A^FkU  
        for j=2:JE-1 0<v5_pB  
            for k=2:KE-1 (bv{17K  
                ex(i,j,k)=gax(i,j,k)*(dx(i,j,k)-ix(i,j,k)); ./;uhj  
                ix(i,j,k)=ix(i,j,k)+gbx(i,j,k)*ex(i,j,k); qvu1u GCc  
                ?&9=f\/P  
                ey(i,j,k)=gay(i,j,k)*(dy(i,j,k)-iy(i,j,k)); uQc("F  
                iy(i,j,k)=iy(i,j,k)+gby(i,j,k)*ey(i,j,k); b?OA|JqX  
                Gn4b*Y&M]3  
                ez(i,j,k)=gaz(i,j,k)*(dz(i,j,k)-iz(i,j,k)); ^vM6_=g2E%  
                iz(i,j,k)=iz(i,j,k)+gbz(i,j,k)*ez(i,j,k); l4i51S"  
            end ppn  8  
        end '8zd]U  
    end }v}F8}4  
    )nf%S+KV  
    %%% Calculate the Fourier transform of Ex ]/Nt  
    40MKf/9  
    %%% Calculate the Incident field s"#N;  
    for j=1:END-1 @2. :fK  
        h_inc(j)=h_inc(j)+0.5*(e_inc(j)-e_inc(j+1)); v`
QDms,{  
    end ;b65s9n^b  
    @~s5{4  
    %%% Hx nG3SDL#(k  
    for j=1:JE-1 75p9_)>96  
        for k=1:KE-1 sXEIC#rq  
            for i=1:ia i[9gcL"  
                curl_e=ey(i,j,k+1)-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k); M2ex 3
m  
                ihxl(i,j,k)=ihxl(i,j,k)+curl_e; WqefH{PB  
                hx(i,j,k)=fj3(j)*fk3(k)*hx(i,j,k)+fj2(j)*fk2(k)*0.5*(curl_e+fi1(i)*ihxl(i,j,k)); Xj+_
"0 #  
            end !m:WoQ
/  
            for i=ia+1:ib iCpm^XT  
                curl_e=ey(i,j,k+1)-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k); zqt<[=O  
                hx(i,j,k)=fj3(j)*fk3(k)*hx(i,j,k)+fj2(j)*fk2(k)*0.5*curl_e;  ij:a+T  
            end ;taZi
xOH  
            for i=ib+1:IE FJH>P\+  
                ixh=i-ib; vkJyD/;=  
                curl_e=ey(i,j,k+1)-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k);
*LhwIY  
                ihxh(ixh,j,k)=ihxh(ixh,j,k)+curl_e; k-3;3Mq  
                hx(i,j,k)=fj3(j)*fk3(k)*hx(i,j,k)+fj2(j)*fk2(k)*0.5*(curl_e+fi1(i)*ihxh(ixh,j,k)); X=-=z5  
            end BMO,eQcB  
        end Z*w(
{k7]  
    end ]\K?%z
  
    N[O .p]8  
    for i=ia+1:ib |)To 0Z  
        for k=ka+1:kb-1 |-vyhr0  
            hx(i,ja,k)=hx(i,ja,k)+0.5*ez_left(i,k); A@|Z^T:  
            hx(i,jb,k)=hx(i,jb,k)-0.5*ez_right(i,k); :I7qw0?  
        end i_m&qy<v  
    end 9\?&u_ U"  
    "Cxj_V@\  
    for i=ia+1:ib 2|\mBP`ok  
        for j=ja+1:jb-1 >F^$ ' b]  
            hx(i,j,ka)=hx(i,j,ka)-0.5*ey_down(i,j); yn ofDGAf  
            hx(i,j,kb)=hx(i,j,kb)+0.5*ey_up(i,j); !{V`N|0  
        end KL?<lp"  
    end i#YDdz  
        :T$}@& -  
    %%% Hy 5K9W5hA:D  
    for i=1:IE-1 -SD:G]un  
        for k=1:KE-1 Ay6T*Nu`  
            for j=1:ja 0y<9JvN$9  
                curl_e=ez(i+1,j,k)-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k); =.S2gO >  
                ihyl(i,j,k)=ihyl(i,j,k)+curl_e; hoBFC1  
                hy(i,j,k)=fi3(i)*fk3(k)*hy(i,j,k)+fi2(i)*fk2(k)*0.5*(curl_e+fj1(j)*ihyl(i,j,k)); q*R~gEi#yk  
            end )x[=}0C  
            for j=ja+1:jb =CFg~8W  
                curl_e=ez(i+1,j,k)-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k); o/
,%rA4  
                hy(i,j,k)=fi3(i)*fk3(k)*hy(i,j,k)+fi2(i)*fk2(k)*0.5*curl_e; GB `n  
            end rkS'OC  
            for j=jb+1:JE ]b}3f<  
                jyh=j-jb; m8L %!6o  
                curl_e=ez(i+1,j,k)-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k); RgFpc*.T  
                ihyh(i,jyh,k)=ihyh(i,jyh,k)+curl_e; YdvXp/P:|  
                hy(i,j,k)=fi3(i)*fk3(k)*hy(i,j,k)+fi2(i)*fk2(k)*0.5*(curl_e+fj1(j)*ihyh(i,jyh,k)); "Ke_dM  
            end B> i^w1  
        end BYB9M  
    end W'k&DKhTqF  
    7>7n|N  
    for i=ia+1:ib-1 |1ry*~  
        for j=ja+1:jb *tZ3?X[b  
            hy(i,j,ka)=hy(i,j,ka)+0.5*ex_down(i,j); .Sw4{m[g  
            hy(i,j,kb)=hy(i,j,kb)-0.5*ex_up(i,j); k(>J?\iNW  
        end u]B b^[  
    end vr47PM2al  
    P/~dY  
    for j=ja+1:jb jyg>'"W  
        for k=ka+1:kb-1 gjT`<CW  
            hy(ia,j,k)=hy(ia,j,k)-0.5*ez_back(j,k); eWYet2!Q  
            hy(ib,j,k)=hy(ib,j,k)+0.5*ez_front(j,k); $(J)F-DB i  
        end ~O./A-l  
    end aAoAjVNkK  
    nC/T$ #G  
   %%% Hz ~:!&}e5  
    for i=1:IE-1 |OF3J,q  
        for j=1:JE-1 3H"F~_H  
            for k=1:ka b,HXD~=  
                curl_e=ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k); ,trh)ZZYW|  
                ihzl(i,j,k)=ihzl(i,j,k)+curl_e; rf$[8d  
                hz(i,j,k)=fi3(i)*fj3(j)*hz(i,j,k)+fi2(i)*fj2(j)*0.5*(curl_e+fk1(k)*ihzl(i,j,k)); )`5-rm~*  
            end 79Q,XRWh|  
            for k=ka+1:kb xop\W4s_  
                curl_e=ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k); hhjsg?4uL  
                hz(i,j,k)=fi3(i)*fj3(j)*hz(i,j,k)+fi2(i)*fj2(j)*0.5*curl_e; 2s 9U&  
            end 5.0;xz}#y  
            for k=kb+1:KE zSsogAx  
                kzh=k-kb; . VI #  
                curl_e=ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k); qy\SOAh  
                ihzh(i,j,kzh)=ihzh(i,j,kzh)+curl_e; e!i.u
'z  
                hz(i,j,k)=fi3(i)*fj3(j)*hz(i,j,k)+fi2(i)*fj2(j)*0.5*(curl_e+fk1(k)*ihzh(i,j,kzh));  ;XYfw)  
            end }?{. 'Hv0  
        end 3_VWtGQ  
    end *=E4|>Ul,  
    `~=Is.V[  
    for i=ia+1:ib-1 < R@&<E6  
        for k=ka+1:kb EjVB\6,  
            hz(i,ja,k)=hz(i,ja,k)-0.5*ex_left(i,k); Z.>?Dt  
            hz(i,jb,k)=hz(i,jb,k)+0.5*ex_right(i,k); . 55aY~We  
        end _x#r,1V+D  
    end yCg>]6B  
    &E k\  
    for j=ja+1:jb-1 1nLFtiki  
        for k=ka+1:kb kr1^`>O5  
            hz(ia,j,k)=hz(ia,j,k)+0.5*ey_back(j,k); lgAE`Os  
            hz(ib,j,k)=hz(ib,j,k)-0.5*ey_front(j,k); 1^;&?E  
        end ;rF:$37^  
    end e8]mdU{)  
    ER$~kFE2yP  
    %%%  Visualize fields 1 gRR  
    %timestep=int2str(N); ~ YZi"u  
    %tview(:,:)=ez(:,:,kc); oU 8o;zk0  
    %pcolor(tview); H @E-=Ly  
    %shading interp; R?%|RCht1  
    %caxis([-1,1]); D3 E!jQ1  
    %colorbar; b;|^62  
    %axis image; Xg"Mjmr  
    %title(['Ez(i,j,k=kc), time step = ',timestep]); |@)ij c4i  
    %xlabel('i'); CV7.hF<  
    %ylabel('j'); `+@%l*TQ  
    %pause(0.05);   }fxH>79g  
        Ipmr@%~  
if N==T-dT ?at~il$z'  
    eio=e_inc(is);  6Ue6b$xE  
    hio=e_inc(is); :Ws3+OI'm3  
    for i=iao:ibo
3?/}  
        for j=jao:jbo            54LCoG/  
            exi_oup(i,j)=ex(i,j,kbo); 8|Y^z_C  
            eyi_oup(i,j)=ey(i,j,kbo); 4]8PF  
            {k~$\J?.  
            exi_odown(i,j)=ex(i,j,kao); a. 5`Q2  
            eyi_odown(i,j)=ey(i,j,kao); 7,LT4wYH  
            oq<#  
            hxi_oup(i,j)=0.5*(hx(i,j,kbo)+hx(i,j,kbo-1)); J5rR?[i{  
            hyi_oup(i,j)=0.5*(hy(i,j,kbo)+hy(i,j,kbo-1)); vqxTf)ys  
            S"Zs'7dy`  
            hxi_odown(i,j)=0.5*(hx(i,j,kao)+hx(i,j,kao-1)); ]X _&  
            hyi_odown(i,j)=0.5*(hy(i,j,kao)+hy(i,j,kao-1)); o_$r*Z|HG  
        end +Q_Gm3^  
    end D{z=)'/F  
    for i=iao:ibo +l3 vIN
 
        for k=kao:kbo IX3U
\_I#  
            exi_oleft(i,k)=ex(i,jao,k); Yc5$915  
            ezi_oleft(i,k)=ez(i,jao,k); ;[zx'e?!  
            B'0Il"g'  
            exi_oright(i,k)=ex(i,jbo,k); 4Y/!V[  
            ezi_oright(i,k)=ez(i,jbo,k); nocH~bAf2  
            W]*wxzf!5z  
            hxi_oleft(i,k)=0.5*(hx(i,jao,k)+hx(i,jao-1,k)); YGn:_9  
            hzi_oleft(i,k)=0.5*(hz(i,jao,k)+hz(i,jao-1,k)); Hm^p^,}_x  
            (c*Dvpo1  
            hxi_oright(i,k)=0.5*(hx(i,jbo,k)+hx(i,jbo-1,k)); I5"=b}V5  
            hzi_oright(i,k)=0.5*(hz(i,jbo,k)+hz(i,jbo-1,k)); kI;^V  
        end Zd%\x[f9ck  
    end Ua3ERBX{  
    for j=jao:jbo JT-Zo OZ  
        for k=kao:kbo ]V]o%onW  
            eyi_ofront(j,k)=ey(ibo,j,k); G0&w#j  
            ezi_ofront(j,k)=ez(ibo,j,k); BzUx@,  
            (jyJ-qe  
            eyi_oback(j,k)=ey(iao,j,k); BZ T%+s;u9  
            ezi_oback(j,k)=ez(iao,j,k); LOYv%9$0*p  
            (6+0U1[Iz  
            hyi_ofront(j,k)=0.5*(hy(ibo,j,k)+hy(ibo-1,j,k)); Tuy*Df  
            hzi_ofront(j,k)=0.5*(hz(ibo,j,k)+hz(ibo-1,j,k)); +%~g$#tlJo  
            FxT [4  
            hyi_oback(j,k)=0.5*(hy(iao,j,k)+hy(iao-1,j,k)); \%VoX`B  
            hzi_oback(j,k)=0.5*(hz(iao,j,k)+hz(iao-1,j,k)); NV;tsuA|  
        end ^b.J z}  
    end 8ksDXf`.  
end 8UZEC-K  
if N==T *Ee# x!O  
    ero=e_inc(is); s`2o\]  
    hro=h_inc(is); 0nd<6S+fs  
    for i=iao:ibo U\z+{]<<  
        for j=jao:jbo            n(+:l'#H
J  
            exr_oup(i,j)=ex(i,j,kbo); [6tQv<}^
 
            eyr_oup(i,j)=ey(i,j,kbo); Aws TDM  
            ouI0"R&@  
            exr_odown(i,j)=ex(i,j,kao); Q&LkST-i  
            eyr_odown(i,j)=ey(i,j,kao); +Snjb0  
            i\'N1S<D  
            hxr_oup(i,j)=0.5*(hx(i,j,kbo)+hx(i,j,kbo-1)); yiUdUw/  
            hyr_oup(i,j)=0.5*(hy(i,j,kbo)+hy(i,j,kbo-1)); DQ}]'*@?  
            (TQXG^n$gY  
            hxr_odown(i,j)=0.5*(hx(i,j,kao)+hx(i,j,kao-1)); KKzvoc?Bt  
            hyr_odown(i,j)=0.5*(hy(i,j,kao)+hy(i,j,kao-1)); J.d `tiN  
        end ;vx9xs?6  
    end Pg36'aTe%j  
    for i=iao:ibo NlMx!f>b%/  
        for k=kao:kbo zYNJF>^<  
            exr_oleft(i,k)=ex(i,jao,k); 5%$#3LT|  
            ezr_oleft(i,k)=ez(i,jao,k); BkF[nL*|  
            a`uT'g[*  
            exr_oright(i,k)=ex(i,jbo,k); "8I4]
'  
            ezr_oright(i,k)=ez(i,jbo,k); '$VR_N\  
            g/J!U8W"  
            hxr_oleft(i,k)=0.5*(hx(i,jao,k)+hx(i,jao-1,k)); YGO@X(ej,  
            hzr_oleft(i,k)=0.5*(hz(i,jao,k)+hz(i,jao-1,k)); ;sf/tX  
            \4[Ta,;t  
            hxr_oright(i,k)=0.5*(hx(i,jbo,k)+hx(i,jbo-1,k)); hdB[H8Q  
            hzr_oright(i,k)=0.5*(hz(i,jbo,k)+hz(i,jbo-1,k)); H#y"3E<s  
        end GXl?Zg  
    end A0ToX) |C  
    for j=jao:jbo BSY#xe V  
        for k=kao:kbo yiq#p"Hs  
            eyr_ofront(j,k)=ey(ibo,j,k); .%A2  
            ezr_ofront(j,k)=ez(ibo,j,k); Fqv5WoYVf  
            GdU W$.  
            eyr_oback(j,k)=ey(iao,j,k); JNz0!wi  
            ezr_oback(j,k)=ez(iao,j,k); [C6?:'}FA  
            DPOPRi~  
            hyr_ofront(j,k)=0.5*(hy(ibo,j,k)+hy(ibo-1,j,k)); aBnbu vp  
            hzr_ofront(j,k)=0.5*(hz(ibo,j,k)+hz(ibo-1,j,k)); NsHveOK1.  
            O#\>j  
            hyr_oback(j,k)=0.5*(hy(iao,j,k)+hy(iao-1,j,k)); g]#Wve  
            hzr_oback(j,k)=0.5*(hz(iao,j,k)+hz(iao-1,j,k)); :6+~"7T  
        end K;l'IN"N  
    end K;PpS*!  
end S8
VR#  
    CmOb+:4@K  
end V]}b3Y
!(  
eo=ero+sqrt(-1)*eio; bRK9Qt#3  
eo=abs(eo); [)A#9L~s=  
ho=(hro+sqrt(-1)*hio)*exp(-sqrt(-1)*0.5*omiga); %K7}yy&9C  
for i=iao:ibo (Nahtx!/9  
    for j=jao:jbo xJ
hbGK  
        ex_oup(i,j)=(exr_oup(i,j)+sqrt(-1)*exi_oup(i,j))/eo; D,/9rH  
        ey_oup(i,j)=(eyr_oup(i,j)+sqrt(-1)*eyi_oup(i,j))/eo; I\~V0<"jI
 
        [4j;FN Fa  
        ex_odown(i,j)=(exr_odown(i,j)+sqrt(-1)*exi_odown(i,j))/eo; o{p_s0IX;S  
        ey_odown(i,j)=(eyr_odown(i,j)+sqrt(-1)*eyi_odown(i,j))/eo; ise}> A!t  
        val<N293L>  
        hx_oup(i,j)=(hxr_oup(i,j)+sqrt(-1)*hxi_oup(i,j))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; ]r6bJ2  
        hy_oup(i,j)=(hyr_oup(i,j)+sqrt(-1)*hyi_oup(i,j))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; {)qP34rM  
        cG:`Zj~4  
        hx_odown(i,j)=(hxr_odown(i,j)+sqrt(-1)*hxi_odown(i,j))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; v_0!uT5~NE  
        hy_odown(i,j)=(hyr_odown(i,j)+sqrt(-1)*hyi_odown(i,j))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; X;:xGZ-oY  
    end _O%p{t'q<  
end ?xZmm%JF  
for i=iao:ibo sBnPS[Oo  
    for k=kao:kbo 8B/9{8  
        ex_oleft(i,k)=(exr_oleft(i,k)+sqrt(-1)*exi_oleft(i,k))/eo; D$E9%'i
r  
        ez_oleft(i,k)=(ezr_oleft(i,k)+sqrt(-1)*ezi_oleft(i,k))/eo; rR~X>+K  
        ~x:]ch|  
        ex_oright(i,k)=(exr_oright(i,k)+sqrt(-1)*exi_oright(i,k))/eo; #+0R!Y  
        ez_oright(i,k)=(ezr_oright(i,k)+sqrt(-1)*ezi_oright(i,k))/eo; 6,1|y%(f  
        [!mjUsut*  
        hx_oleft(i,k)=(hxr_oleft(i,k)+sqrt(-1)*hxi_oleft(i,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; )g_zPt  
        hz_oleft(i,k)=(hzr_oleft(i,k)+sqrt(-1)*hzi_oleft(i,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; UAZ&*{MM^  
        ?y)X$D^  
        hx_oright(i,k)=(hxr_oright(i,k)+sqrt(-1)*hxi_oright(i,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; %LHV0u  
        hz_oright(i,k)=(hzr_oright(i,k)+sqrt(-1)*hzi_oright(i,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; QOk"UP  
    end ':]a.yA\1  
end 1,'^BgI,  
for j=jao:jbo &NeYKh?  
    for k=kao:kbo "lf3hWGw  
        ey_ofront(j,k)=(eyr_ofront(j,k)+sqrt(-1)*eyi_ofront(j,k))/eo; ]H~,K]@.  
        ez_ofront(j,k)=(ezr_ofront(j,k)+sqrt(-1)*ezi_ofront(j,k))/eo; YWrY{6M  
        o q)"1  
        ey_oback(j,k)=(eyr_oback(j,k)+sqrt(-1)*eyi_oback(j,k))/eo; Els=:4  
        ez_oback(j,k)=(ezr_oback(j,k)+sqrt(-1)*ezi_oback(j,k))/eo; 0\@|M@X=  
        l$_rA~Mo  
        hy_ofront(j,k)=(hyr_ofront(j,k)+sqrt(-1)*hyi_ofront(j,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; :w%bw\}  
        hz_ofront(j,k)=(hzr_ofront(j,k)+sqrt(-1)*hzi_ofront(j,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z;  LB7I`W  
        8{%9%{  
        hy_oback(j,k)=(hyr_oback(j,k)+sqrt(-1)*hyi_oback(j,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; M) XQi/  
        hz_oback(j,k)=(hzr_oback(j,k)+sqrt(-1)*hzi_oback(j,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; %0]b5u  
    end `|Z@UPHzG  
end Eo&qc 17)`  
%%%%%%%%%  Far Field Calculate  %%%%%%% -U[`pUY?f  
rk=2*pi/WL; 5O ;^Mk|  
fx=zeros(1,181); )%w8>1}c  
fy=zeros(1,181); Lu.tRZ`$38  
fz=zeros(1,181); !13 /+ u  
fmx=zeros(1,181); .WyX/E$I^!  
fmy=zeros(1,181); y
4rJ-  
fmz=zeros(1,181); d ovwB`5  
rcs=zeros(1,181); B-'BJ|*4I  
est=zeros(1,181); [(EH  
efa=zeros(1,181); (Gc`3jJ  
fai_s=0*pi/180; _6zP]|VBr  
cps2=cos(fai_s); 6RQCKN)  
sps2=sin(fai_s); eAKQR  
ss=0; BV?N_/DXp  
for sita_s=0:pi/180:pi sNmC#,  
    ss=ss+1; @i-@mxk6<  
    sps1=sin(sita_s); =2]
rA  
    cps1=cos(sita_s); ]V^.!=gh$  
    rk=2*pi/WL; +gCy@_2;  
    %%% fx、fmx——left/right output boundry {}&f\6OI%  
    for i=iao+1:ibo 9~%]|_(  
        for k=kao:kbo (A fbS=[  
            temp1=exp(sqrt(-1)*rk*((i-0.5)*sps1*cps2+jao*sps1*sps2+k*cps1)); L4*fF  
            fx(ss)=fx(ss)-(1-0.5*(k==kao)|(k==kbo))*ddx*ddx*hz_oleft(i,k)*temp1; x~Ly$A2p  
            fmx(ss)=fmx(ss)+(1-0.5*(k==kao)|(k==kbo))*ddx*ddx*ez_oleft(i,k)*temp1; edk9Qd9  
            >@

"3Q`  
            temp2=exp(sqrt(-1)*rk*((i-0.5)*sps1*cps2+jbo*sps1*sps2+k*cps1)); yi-"hT`  
            fx(ss)=fx(ss)+(1-0.5*(k==kao)|(k==kbo))*ddx*ddx*hz_ori .. ,xe@G)a  
RdvTtXg  

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

大段的代码肯定难看，你可以尝试以下改正方法： 5!'R'x5e  
1，低级失误，i,j,k混淆,x,y,z混淆 ]g;^w?9h  
2，看动态显示图，这是最有效的方法，哪里出问题直接定位到处理这部分网格的程序中。 \Y!T>nWn)I  
3，运行一次，观察MATLAB中的各变量，按经验判断错误的计算值

好的，谢谢！ 31WZJm^  
q1H=/[a  
根据你的提醒调试，希望能找出错误。 $fj])>=H