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

       下面是我整的三维时谐场——用FDTD计算介质球的RCS的程序，但是在近远场外推的时候出现了问题，怎么都跑不出图来， Muarryh}  
希望大家能帮我找出错误，谢谢！ )hn,rmn (P  
       IE=81; ,@<-h* m  
JE=81; M>0~Ek%3  
KE=81; Tvk=NJ  
ia=9; TsR20P@  
ja=9; 0wOgQ
n  
ka=9; [T
NYPA>{  
ic=(IE+1)/2; Qi[T!1  
jc=(JE+1)/2;
SH5k^EJ  
kc=(KE+1)/2; p4t(xm2T  
ib=IE-ia; |=7%Edkd  
jb=JE-ja; AjZT- Q0L  
kb=KE-ka; t);5Cw
_  
epsz=8.8e-12; BURiLEYZl  
muz=4*pi*1e-7; RI%*5lM8;  
lamda=0.032; ?lbX.+  
freq=lamda/3e+8; ~L?p/3m   
%freq=30e+7; u_mm*o~)g  
ddx=lamda/42; =r GkM.^  
dt=ddx/6e+8;
d*H-l3N  
arg=2*pi*freq*dt; Q>X ;7nt0  
dx=zeros(IE,JE,KE); h; {?z  
dy=zeros(IE,JE,KE); X_=oJi|:  
dz=zeros(IE,JE,KE); !4 6^}3  
ex=zeros(IE,JE,KE); 3?fya8W<  
ey=zeros(IE,JE,KE); E{j6OX\  
ez=zeros(IE,JE,KE); S:DB%V3  
hx=zeros(IE,JE,KE); \Z)''
:},C  
hy=zeros(IE,JE,KE); LxMOs Nv  
hz=zeros(IE,JE,KE); ^ 41p+  
ix=zeros(IE,JE,KE); @UvjJ  
iy=zeros(IE,JE,KE); J :,  
iz=zeros(IE,JE,KE); C8(sH@  
gax=zeros(IE,JE,KE); h7o{l7`)  
gay=zeros(IE,JE,KE); B:9Z;g@&  
gaz=zeros(IE,JE,KE); OY?uqP}c  
gbx=zeros(IE,JE,KE); @ cv`}k  
gby=zeros(IE,JE,KE); q.MM|;_u`  
gbz=zeros(IE,JE,KE); FmnA+fA  
END=200; J\2F%kBej?  
e_inc=zeros(1,END); 9WBDSx_(Q  
h_inc=zeros(1,END); |z5olu$gVc  
e_temp=zeros(1,4); dedi6Brl  
idxl=zeros(ia,JE,KE); ujwI4oj"c  
idxh=zeros(ia,JE,KE); |QHWX^pO  
ihxl=zeros(ia,JE,KE); W>+<r9Rt4  
ihxh=zeros(ia,JE,KE); }De)_E\~  
idyl=zeros(IE,ja,KE); h#;yA"j1&  
idyh=zeros(IE,ja,KE); \2~.r/`1  
ihyl=zeros(IE,ja,KE); hf%W grO.  
ihyh=zeros(IE,ja,KE); fV[xv4D.  
idzl=zeros(IE,JE,ka); p-ry{"XA  
idzh=zeros(IE,JE,ka); SqEO ]~  
ihzl=zeros(IE,JE,ka); 4N*^%  
ihzh=zeros(IE,JE,ka); A~h8 >zz*  
gi1=zeros(1,IE);gi2=ones(1,IE);gi3=ones(1,IE); )TXn7{M:  
gj1=zeros(1,JE);gj2=ones(1,JE);gj3=ones(1,JE); slw^BK3t  
gk1=zeros(1,KE);gk2=ones(1,KE);gk3=ones(1,KE); hI/p9 `w  
fi1=zeros(1,IE);fi2=ones(1,IE);fi3=ones(1,IE); W&rjJZY6  
fj1=zeros(1,JE);fj2=ones(1,JE);fj3=ones(1,JE); GhQ.}@*  
fk1=zeros(1,KE);fk2=ones(1,KE);fk3=ones(1,KE); L9G=+T9  
sita=0*pi/180;    %pointing angle !M}&dW2  
fai=0*pi/180;     %azimuth angle f!1KGP  
arfa=0*pi/180;    %polarizing angle 23Eg|Xk  
sp1=sin(sita);cp1=cos(sita); v^KJU +  
sp2=sin(fai);cp2=cos(fai); kV-a'"W5  
sp3=sin(arfa);cp3=cos(arfa); D4`7,JC}<  
ex_left=zeros(IE,KE);ex_right=zeros(IE,KE); W/J3sAYv  
ex_up=zeros(IE,JE);ex_down=zeros(IE,JE); .]8 Jeb  
ey_front=zeros(JE,KE);ey_back=zeros(JE,KE); xXLKL6F(\  
ey_up=zeros(IE,JE);ey_down=zeros(IE,JE); @ '<lD*W  
ez_left=zeros(IE,KE);ez_right=zeros(IE,KE); Ph-3,cC  
ez_front=zeros(JE,KE);ez_back=zeros(JE,KE); >o(*jZ  
hx_left=zeros(IE,KE);hx_right=zeros(IE,KE); c L84}1QD  
hx_up=zeros(IE,JE);hx_down=zeros(IE,JE); R|tjvp-[}  
hy_front=zeros(JE,KE);hy_back=zeros(JE,KE); q!NwfXJM  
hy_up=zeros(IE,JE);hy_down=zeros(IE,JE); w}7`Vas9  
hz_left=zeros(IE,KE);hz_right=zeros(IE,KE); t6LTGWs/_o  
hz_front=zeros(JE,KE);hz_back=zeros(JE,KE); ;r8,Wx@f1C  
npml=4;        %pml层的厚度 ysvn*9h+&  
xxn=1:-1/npml:1/npml; viKN:n! Ev  
xn=0.3333*xxn.^3; d{DlW |_  
fi1(1:npml)=xn; {%~Sbcq4F  
fi1(IE:-1:IE-npml+1)=xn; b
}[{'  
gi2(1:npml)=1./(1+xn); ]\Ez{MdAT  
gi2(IE:-1:IE-npml+1)=1./(1+xn); .i`+}@iA  
gi3(1:npml)=(1-xn)./(1+xn); .H9!UQ&It  
gi3(IE:-1:IE-npml+1)=(1-xn)./(1+xn); $z` jR*  
fj1(1:npml)=xn; tXfB.[U  
fj1(JE:-1:JE-npml+1)=xn; z#O{rwnl  
gj2(1:npml)=1./(1+xn); egKYlfe"  
gj2(JE:-1:JE-npml+1)=1./(1+xn); _{T`
ka  
gj3(1:npml)=(1-xn)./(1+xn); hXYVi6(k  
gj3(JE:-1:JE-npml+1)=(1-xn)./(1+xn); qB"y'UW8  
fk1(1:npml)=xn; R(x%<I  
fk1(KE:-1:KE-npml+1)=xn; |[xi"E\  
gk2(1:npml)=1./(1+xn); en#W<"_"
  
gk2(KE:-1:KE-npml+1)=1./(1+xn); r?H {Y3,  
gk3(1:npml)=(1-xn)./(1+xn); EEF}Wf$f  
gk3(KE:-1:KE-npml+1)=(1-xn)./(1+xn); g2F~0%HY  
xxn=(1-0.5/npml):-1/npml:0.5/npml; nEs l  
xn=0.3333*xxn.^3; 60QElJ9D  
gi1(1:npml)=xn; m@"!=CTKd  
gi1(IE-1:-1:IE-npml)=xn; w6Ue5Ix,!  
fi2(1:npml)=1./(1+xn); ox)/*c<  
fi2(IE-1:-1:IE-npml)=1./(1+xn); y:pypuwt;  
fi3(1:npml)=(1-xn)./(1+xn); 'O2{0  
fi3(IE-1:-1:IE-npml)=(1-xn)./(1+xn); $^`hu%s,~  
gj1(1:npml)=xn; qOkw6jfluh  
gj1(JE-1:-1:JE-npml)=xn; I7]45pF  
fj2(1:npml)=1./(1+xn); drF"kTD"7  
fj2(JE-1:-1:JE-npml)=1./(1+xn); r`6XF  
fj3(1:npml)=(1-xn)./(1+xn); D|UDLaz~  
fj3(JE-1:-1:JE-npml)=(1-xn)./(1+xn); ,{YC|uB  
gk1(1:npml)=xn; g"<kj"  
gk1(KE-1:-1:KE-npml)=xn; Ip?Ueaei  
fk2(1:npml)=1./(1+xn); W5sVQ`S-  
fk2(KE-1:-1:KE-npml)=1./(1+xn); d6~wJMFl  
fk3(1:npml)=(1-xn)./(1+xn); !L>
'g  
fk3(KE-1:-1:KE-npml)=(1-xn)./(1+xn); w=O:|Xu#*  
WL=42; >oLM2VJ  
dT=WL/2; $8X?|fV)  
omiga=pi/WL; .\< \J|3  
Z=sqrt(muz/epsz); 3",6 E(  
odist=2; O+(Z`,^  
iao=ia+1-odist; hLO nX<%a  
jao=ja+1-odist; %lEPFp  
kao=ka+1-odist; JQI`9$asuC  
ibo=ib+odist; \5X34'7  
jbo=jb+odist; ./.E=,j  
kbo=kb+odist; 6V!yfps)  
exi_oup=zeros(IE,JE);exr_oup=zeros(IE,JE);ex_oup=zeros(IE,JE); ca$D|3  
exi_odown=zeros(IE,JE);exr_odown=zeros(IE,JE);ex_odown=zeros(IE,JE); CYG'WFvZZ  
eyi_oup=zeros(IE,JE);eyr_oup=zeros(IE,JE);ey_oup=zeros(IE,JE); zoO>N'b3)  
eyi_odown=zeros(IE,JE);eyr_odown=zeros(IE,JE);ey_odown=zeros(IE,JE); >uxak2nM-  
exi_oleft=zeros(IE,KE);exr_oleft=zeros(IE,KE);ex_oleft=zeros(IE,KE); *:)#'cenI  
ezi_oleft=zeros(IE,KE);ezr_oleft=zeros(IE,KE);ez_oleft=zeros(IE,KE); &a1ag
i7M  
exi_oright=zeros(IE,KE);exr_oright=zeros(IE,KE);ex_oright=zeros(IE,KE); Cc9<ABv?  
ezi_oright=zeros(IE,KE);ezr_oright=zeros(IE,KE);ez_oright=zeros(IE,KE); c1b@3  
eyi_ofront=zeros(JE,KE);eyr_ofront=zeros(JE,KE);ey_ofront=zeros(JE,KE); qCIZW  
ezi_ofront=zeros(JE,KE);ezr_ofront=zeros(JE,KE);ez_ofront=zeros(JE,KE); OB5(4TY  
eyi_oback=zeros(JE,KE);eyr_oback=zeros(JE,KE);ey_oback=zeros(JE,KE); Cf8(Jk`v|  
ezi_oback=zeros(JE,KE);ezr_oback=zeros(JE,KE);ez_oback=zeros(JE,KE); &Gxk
~p<  
hxi_oup=zeros(IE,JE);hxr_oup=zeros(IE,JE);hx_oup=zeros(IE,JE); vFy/  
hxi_odown=zeros(IE,JE);hxr_odown=zeros(IE,JE);hx_odown=zeros(IE,JE); .LV=Z0ja  
hyi_oup=zeros(IE,JE);hyr_oup=zeros(IE,JE);hy_oup=zeros(IE,JE); MI<hShc\  
hyi_odown=zeros(IE,JE);hyr_odown=zeros(IE,JE);hy_odown=zeros(IE,JE); ]ujH7T  
hxi_oleft=zeros(IE,KE);hxr_oleft=zeros(IE,KE);hx_oleft=zeros(IE,KE); iZ)7%R?5  
hzi_oleft=zeros(IE,KE);hzr_oleft=zeros(IE,KE);hz_oleft=zeros(IE,KE); lD\lFN(:
 
hxi_oright=zeros(IE,KE);hxr_oright=zeros(IE,KE);hx_oright=zeros(IE,KE); wgfn:LR  
hzi_oright=zeros(IE,KE);hzr_oright=zeros(IE,KE);hz_oright=zeros(IE,KE); *}3~8fu{  
hyi_ofront=zeros(JE,KE);hyr_ofront=zeros(JE,KE);hy_ofront=zeros(JE,KE); @4hxGk=
 
hzi_ofront=zeros(JE,KE);hzr_ofront=zeros(JE,KE);hz_ofront=zeros(JE,KE); )FE'#\  
hyi_oback=zeros(JE,KE);hyr_oback=zeros(JE,KE);hy_oback=zeros(JE,KE); Fx)]AJ~[t  
hzi_oback=zeros(JE,KE);hzr_oback=zeros(JE,KE);hz_oback=zeros(JE,KE); ]+ug:E{7  
%%%介质球的设置 awo=%vJ&  
epsilon=[1 4]; kU4Zij-O  
%epsilon=[1 81]; l{Xsh;%=  
sigma=[0 0]; IRXpk6|  
%sigma=[0 5]; hnH:G`[F  
%sigma=[0 0]; 6lsU/`.  
radius=20; V?%>Ex$  
for i=ia+1:ib Y`\zLX"_m  
    for j=ja+1:jb &`t-[5O\  
        for k=ka+1:kb l&
5| =  
            eps=epsilon(1); N4'b]:`n  
            cond=sigma(1); 2VMau.eQ  
            xdist=(ic-i-0.5); '$l*FWOEal  
            ydist=jc-j; r=8]Ub[  
            zdist=kc-k; EGU?54  
            dist=sqrt(xdist^2+ydist^2+zdist^2); R3%T}^;f  
            if dist<=radius Yb414K  
                eps=epsilon(2); &<k)W  
                cond=sigma(2); \-
`L}$  
            end ENjrv  
            gax(i,j,k)=1/(eps+cond*dt/epsz); sK}Ru?a)  
            gbx(i,j,k)=cond*dt/epsz; m,qU})  
        end Xs*~[k'  
    end 7{/qQGL  
end m6aoh^I  
for i=ia+1:ib AaKILIIQZ  
    for j=ja+1:jb ,tmo6D62  
        for k=ka+1:kb b&V
=X{V4  
            eps=epsilon(1); TkmN.@w_C  
            cond=sigma(1); ,_P(!7Z8  
            xdist=ic-i; O:pQf/Xn  
            ydist=jc-j-0.5; !X*L<)=nh  
            zdist=kc-k; vAU^<$D27  
            dist=sqrt(xdist^2+ydist^2+zdist^2); Qy$8!(  
            if dist<=radius $kkdB,y  
                eps=epsilon(2); &UQKZ.  
                cond=sigma(2); .r[J} O"  
            end Cg8{NNeD  
            gay(i,j,k)=1/(eps+cond*dt/epsz); 6'|J ;  
            gby(i,j,k)=cond*dt/epsz; V|<qO-#.  
        end G8f7N;D  
    end T\. 8og  
end $F;$-2  
for i=ia+1:ib /vSFQ}W  
    for j=ja+1:jb Y==# yNwM  
        for k=ka+1:kb :IbrV@gN{@  
            eps=epsilon(1); *s;$`8fM<  
            cond=sigma(1); yu3EPT!~  
            xdist=ic-i; goWD~'\  
            ydist=jc-j; 0Rrz   
            zdist=kc-k-0.5; hq(3%- 7&  
            dist=sqrt(xdist^2+ydist^2+zdist^2); l.fNkLC#  
            if dist<=radius HwM:bY N  
                eps=epsilon(2); eAm7*2  
                cond=sigma(2); -&@[]/  
            end 5#q ^lL
  
            gaz(i,j,k)=1/(eps+cond*dt/epsz); L.ndLd  
            gbz(i,j,k)=cond*dt/epsz; v>7tJ[s  
        end />.&  
    end 7u o4F=%  
end >-EoE;s  
real_in=0; DEqk9Exk`  
imag_in=0; 9`-ofwr'|  
%t0=40; &> }MoB  
%spread=10; nolTvqMT  
T=1000; z[IG+2  
%N=0; OlMCF.W#3  
%isource=floor(sqrt((ib-ia)^2+(jb-ja)^2+(kb-ka)^2)); `[57U,v  
isource=4; z}OY'}sk8  
is=ja+1; &=@R,  
for N=1:T 9CUMqaY2  
    8I NVn'G  
    for j=2:END /.7x[Yc  
        e_inc(j)=e_inc(j)+0.5*(h_inc(j-1)-h_inc(j)); C'*1
w  
    end ",Ek| z  
    z&cfFx#h)  
    %real_in=real_in+cos(arg*T)*e_inc(ja-1); ;bkS0Vmg  
    %imag_in=imag_in-sin(arg*T)*e_inc(ja-1); 7A3e-51>  
    {%b>/r  
    %pulse=exp(-0.5*((t0-T)/spread)^2); I`DdhMi7  
    pulse=sin(arg*N); PY)C=={p  
    e_inc(isource)=pulse; ,7>_Lp_v  
    KDg!Y(m{  
    e_inc(1)=e_temp(1); 2zArAch  
    e_temp(1)=e_temp(2); y)Ip\.KV\  
    e_temp(2)=e_inc(2); FSc730rM  
    e_inc(END)=e_temp(3); lT1*e(I  
    e_temp(3)=e_temp(4); ^ chlAQz(  
    e_temp(4)=e_inc(END-1); ax7ub  
    $5%tGFh  
    %%%  Calculate Incident component on connected edge %D e<H*  
    %%%  up/down connect edge LPewoAXO  
    for i=ia+1:ib C@hnT<e  
        for j=ja+1:jb "R4~ 8r  
            k=kb; $N:m 9R  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; 8Bo'0  
            w=dr-floor(dr); _S@s  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); dpGaI  
            ex_up(i,j)=ein*(sp2*cp3-cp1*cp2*sp3); =xai 7iM  
            [ivJ&'vB  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; D-~
HJ  
            w=dr-floor(dr); 0`,a@Q4  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); &'T7 ~M:  
            ey_up(i,j)=ein*(-cp2*cp3-cp1*sp2*sp3); G4=%<+
  
            pd>EUdbrp&  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; ~EE*/vX  
            w=dr-floor(dr); }*IX3
4  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); [<8<+lH=P  
            hx_up(i,j)=hin*(sp2*sp3+cp1*cp2*cp3); 1B,RRHXn6  
            >:KPvq!0  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1;  QqtFNG  
            w=dr-floor(dr); ~cU,3g  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); WoesE:NiR  
            hy_up(i,j)=hin*(-cp2*sp3+cp1*sp2*cp3); W53i5u(  
        end "AUY+ LN  
    end oA_AnD?G+  
    2$\Du9+  
    for i=ia+1:ib H)t8d_^|j  
        for j=ja+1:jb m'z<d  
            k=ka+1; @iao"&  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; A-c3B+  
            w=dr-floor(dr); VEE:Z^U!  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); EC?!%iO`  
            ex_down(i,j)=ein*(sp2*cp3-cp1*cp2*sp3); wyvs#T  
            pz.<5  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; ,i,f1XJ|  
            w=dr-floor(dr); /of,4aaK7  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); [eG- &u  
            ey_down(i,j)=ein*(-cp2*cp3-cp1*sp2*sp3); evu@uq  
            R"=G?d)  
            k=ka; y4P mL  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; YNk?1#k?i  
            w=dr-floor(dr); PNg,
bcl  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); J|,| *t  
            hx_down(i,j)=hin*(sp2*sp3+cp1*cp2*cp3); aHwrFkn  
            is#?O5:2  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; ^J@Y?CQl\  
            w=dr-floor(dr); VGq2ITg9eE  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); -L1{0{Z  
            hy_down(i,j)=hin*(-cp2*sp3+cp1*sp2*cp3);  1=W>zC  
        end ?}B_'NZ%  
    end xxiEL2"`>  
    ]5uCs[  
    %%%  left/right connect edge cz&FOP+!  
    for i=ia+1:ib 0^:O:X  
        for k=ka+1:kb jc)D*Cf  
            j=ja+1; .NwHr6/s*  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; $;$_N43  
            w=dr-floor(dr); |/YT.c%  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); :N%]<Mq  
            ex_left(i,k)=ein*(sp2*cp3-cp1*cp2*sp3); 7NoB   
            nyT
fTn  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; F0Rk[GM
 
            w=dr-floor(dr); Jw"'ZW#W  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); 'rq [P",  
            ez_left(i,k)=ein*(sp1*sp3); khW3z*e#  
            J}%&;uv  
            j=ja; wQ4/eQ*  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; 9K FWa0G  
            w=dr-floor(dr); AGPZd9  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); H[!by)H  
            hx_left(i,k)=hin*(sp2*sp3+cp1*cp2*cp3); o`hF1*yp  
            ,v>P05  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; 'OkF.bs  
            w=dr-floor(dr); j YO#  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); T9N /;3  
            hz_left(i,k)=hin*(-sp1*cp3); E]_sl/`{od  
        end pbloL3d.;+  
    end J"I{0>@  
       S(9fGh  
    for i=ia+1:ib OW1[Y-o[  
        for k=ka+1:kb a?\ Au  
            j=jb; A!goR-J]  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; ;Fp"]z!Qh+  
            w=dr-floor(dr); '.d el7s  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); q$>At}4  
            ex_right(i,k)=ein*(sp2*cp3-cp1*cp2*sp3); =,,!a/U  
            viY _Y.Yjy  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1;  A5Y z|  
            w=dr-floor(dr); ,%>/8*  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); S%g`X  
            ez_right(i,k)=ein*(sp1*sp3); bBC3% H^  
            <:>a51HBX  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; LI"ghz=F  
            w=dr-floor(dr); /,GDG=ra  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); :|N5fkhN  
            hx_right(i,k)=hin*(sp2*sp3+cp1*cp2*cp3); 4}fG{B
k  
            e9N"{kDs6  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; 5BTQJa  
            w=dr-floor(dr); gnx!_H\h<  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); L[Tr"BW  
            hz_right(i,k)=hin*(-sp1*cp3); b^6Ooc/-k  
        end ~hYG%  
    end sq-[<ryk  
    H-^>Co_  
    %%%  front/back connect edge u@[D*c1!H  
    for j=ja+1:jb 3 LoB-4u?  
        for k=ka+1:kb om;jXf}A  
            i=ib; v3
4XcA  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; BEifUgCh  
            w=dr-floor(dr); dhsQfWg#}  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); ]NG`MZ  
            ey_front(j,k)=ein*(-cp2*cp3-cp1*sp2*sp3); <E!M<!h  
            <_ddGg~  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; J?u@' "u  
            w=dr-floor(dr); mqw&SxU9  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); &*[T
 
            ez_front(j,k)=ein*(sp1*sp3); l9#M`x9  
            nwHi3ojD:  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; |c dQJW  
            w=dr-floor(dr); n\wO[l)  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); F0+u#/#  
            hy_front(j,k)=hin*(-cp2*sp3+cp1*sp2*cp3); @ 5|F:J  
            ` *h-j/M  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; fh}\#WE"  
            w=dr-floor(dr); w*\)]bTs  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); iI&J_Y{1a_  
            hz_front(j,k)=hin*(-sp1*cp3); f!0*^d  
        end |HhUU1!  
    end `NyvJt^<  
    yH@W6'.  
    for j=ja+1:jb /&cb`^"U^  
        for k=ka+1:kb b":cj:mxL  
            i=ia+1; #R'm|En'  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; 0rjx
WPc  
            w=dr-floor(dr); 7 45Uo'  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); [bH6>{3u  
            ey_back(j,k)=ein*(-cp2*cp3-cp1*sp2*sp3); &37QUdp+p  
            "4oY F:h  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; 8L6!CP_!  
            w=dr-floor(dr); W@JmG`Sy  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); N3 07lGb  
            ez_back(j,k)=ein*(sp1*sp3); |*i0h`a  
            3dQV5E.  
            i=ia; $iupzVrro  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; '-S^z"ZrI  
            w=dr-floor(dr); UIovv%7zZ  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); *<"xF'C  
            hy_back(j,k)=hin*(-cp2*sp3+cp1*sp2*cp3); IF$f^$  
            #3\F<AJ<VB  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; O: :X$O7  
            w=dr-floor(dr); ixE72bX  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); NuLQkf)  
            hz_back(j,k)=hin*(-sp1*cp3); 3-Xum*)Y  
        end YeptYW@xfw  
    end H'&x
4[J:  
    j<Lj1P3  
    %%% Dx /mo
4Q?^  
    for j=2:JE ,*}SfCon  
        for k=2:KE D CcM~  
            for i=2:ia 31Y+bxQ  
                curl_h=hz(i,j,k)-hz(i,j-1,k)-hy(i,j,k)+hy(i,j,k-1); aOA;"jR1  
                idxl(i,j,k)=idxl(i,j,k)+curl_h; uiPfAPZ  
                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)); q.g<gu]  
            end 3?.3Z!H/  
            for i=ia+1:ib -[".km  
                curl_h=hz(i,j,k)-hz(i,j-1,k)-hy(i,j,k)+hy(i,j,k-1); `N]!-=o  
                dx(i,j,k)=gj3(j)*gk3(k)*dx(i,j,k)+gj2(j)*gk2(k)*0.5*curl_h; XGbtmmQG  
            end E5/-?(N  
            for i=ib+1:IE TB1E1  
                ixh=i-ib; p4*VE5[?_+  
                curl_h=hz(i,j,k)-hz(i,j-1,k)-hy(i,j,k)+hy(i,j,k-1); yL7a*C&  
                idxh(ixh,j,k)=idxh(ixh,j,k)+curl_h; 1|q$Wn:*  
                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)); %q`_vtUT  
            end CES^ c-. k  
        end Ky|dRbK,  
    end +F]X  
   K8ecSs}}J  
    for i=ia+1:ib-1 q6%jCt2'  
        for k=ka+1:kb D42Bm&JocO  
            dx(i,ja+1,k)=dx(i,ja+1,k)-0.5*hz_left(i,k); 5BN!uUkm+  
            dx(i,jb,k)=dx(i,jb,k)+0.5*hz_right(i,k); cAC]%~orx  
        end "z1\I\ ^  
    end A<ur20  
    {I&>`?7.  
    for i=ia+1:ib-1 wtu WzHrF  
        for j=ja+1:jb wm=!tx\`k  
            dx(i,j,ka+1)=dx(i,j,ka+1)+0.5*hy_down(i,j); WIa4!\Ky!  
            dx(i,j,kb)=dx(i,j,kb)-0.5*hy_up(i,j); 4 B"tz!  
        end ,mx>)}l95  
    end 7D@O:yO  
    [x kbzJ  
    %%% Dy %J5zfNe)&  
    for i=2:IE am]$`7R5d  
        for k=2:KE /608P:U  
            for j=2:ja Uw8O"}U8  
                curl_h=hx(i,j,k)-hx(i,j,k-1)-hz(i,j,k)+hz(i-1,j,k); ts<5%{M(  
                idyl(i,j,k)=idyl(i,j,k)+curl_h; J/:
9;{R  
                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)); jn&[=Y-  
            end Ugp[Ugr  
            for j=ja+1:jb =EA*h_"q9  
                curl_h=hx(i,j,k)-hx(i,j,k-1)-hz(i,j,k)+hz(i-1,j,k); y9cW&rDH  
                dy(i,j,k)=gi3(i)*gk3(k)*dy(i,j,k)+gi2(i)*gk2(k)*0.5*curl_h; U^-:qT;CX  
            end oe=^CeW"  
            for j=jb+1:JE
4. 7m*  
            jyh=j-jb; E=_M=5]  
            curl_h=hx(i,j,k)-hx(i,j,k-1)-hz(i,j,k)+hz(i-1,j,k); %V%*0S|U  
            idyh(i,jyh,k)=idyh(i,jyh,k)+curl_h; )[.FUx  
            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)); \25Rq/&w  
            end T<=Ci?C v  
        end C4H M  
    end K.Nun)<  
    Gg]>S#^3  
    for i=ia+1:ib 5A g4o  
        for j=ja+1:jb-1 "7(2m  
            dy(i,j,ka+1)=dy(i,j,ka+1)-0.5*hx_down(i,j); PO%Z.ol9  
            dy(i,j,kb)=dy(i,j,kb)+0.5*hx_up(i,j); n,!PyJ  
        end 6}qp;mR E]  
    end 8$xd;+`y'  
    F.<sKQ&A  
    for j=ja+1:jb-1 F4xYfbwY"]  
        for k=ka+1:kb R^.E";/h  
            dy(ia+1,j,k)=dy(ia+1,j,k)+0.5*hz_back(j,k); -NHA{?6r  
            dy(ib,j,k)=dy(ib,j,k)-0.5*hz_front(j,k); UA>UW!I  
        end P:eY>~m<;  
    end ,_,Z<X/
 
   f02<u  
    %%%  Dz sgW*0o  
    for i=2:IE s"7wG!yf  
        for j=2:JE 5 ?~ ?8Hi  
            for k=2:ka ae`*0wbv  
                curl_h=hy(i,j,k)-hy(i-1,j,k)-hx(i,j,k)+hx(i,j-1,k); 0,a\vs%@X  
                idzl(i,j,k)=idzl(i,j,k)+curl_h; TRW{`b[  
                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)); h:\WW;s[B  
            end $p#)xx7  
            for k=ka+1:kb s.d }*H-o  
                curl_h=hy(i,j,k)-hy(i-1,j,k)-hx(i,j,k)+hx(i,j-1,k); 6 9_etv  
                dz(i,j,k)=gi3(i)*gj3(j)*dz(i,j,k)+gi2(i)*gj2(j)*0.5*curl_h; .bE+dA6:v  
            end pTT7#b(t  
            for k=kb+1:KE >`R}ulz)  
                kzh=k-kb; ..`J-k  
                curl_h=hy(i,j,k)-hy(i-1,j,k)-hx(i,j,k)+hx(i,j-1,k); vq-Tq>  
                idzh(i,j,kzh)=idzh(i,j,kzh)+curl_h; 3J%(2}{y  
                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)); 56T<s+X>  
            end xE`uFHuS}  
        end 1S/KT4  
    end =~}\g;K1Q  
    *p0n{F9  
    for i=ia+1:ib b+bgGLo  
        for k=ka+1:kb-1 ZCsL%(  
            dz(i,ja+1,k)=dz(i,ja+1,k)+0.5*hx_left(i,k); `irz'/"p  
            dz(i,jb,k)=dz(i,jb,k)-0.5*hx_right(i,k); r%\(5H f  
        end \bNN]=  
    end dtD)VNkBZ  
    wl|cipy"  
    for j=ja+1:jb ByqB4Hv2  
        for k=ka+1:kb-1 96x$Xl;  
            dz(ia+1,j,k)=dz(ia+1,j,k)-0.5*hy_back(j,k); 8,d<&3D  
            dz(ib,j,k)=dz(ib,j,k)+0.5*hy_front(j,k); *{5p/}p  
        end s tvI  
    end 29p
`G1n  
    Eh/B[u7T[  
    %%%  Calculate E from D =|_:H$94  
    for i=2:IE-1 fG_.&!P  
        for j=2:JE-1 ^NwXvp>7-  
            for k=2:KE-1 5i `q  
                ex(i,j,k)=gax(i,j,k)*(dx(i,j,k)-ix(i,j,k)); Z.TYi~d/9D  
                ix(i,j,k)=ix(i,j,k)+gbx(i,j,k)*ex(i,j,k); 5$9g4  
                }q7rR:g  
                ey(i,j,k)=gay(i,j,k)*(dy(i,j,k)-iy(i,j,k)); T
JjcX?:(  
                iy(i,j,k)=iy(i,j,k)+gby(i,j,k)*ey(i,j,k); I<
yd=#:n  
                {=};<;_F  
                ez(i,j,k)=gaz(i,j,k)*(dz(i,j,k)-iz(i,j,k)); 53=5xE= `D  
                iz(i,j,k)=iz(i,j,k)+gbz(i,j,k)*ez(i,j,k); @Cg%7AF  
            end nQbF~  
        end Z3X/SQ'0  
    end ?"-%>y@w  
    {/d<Jm:  
    %%% Calculate the Fourier transform of Ex &Z3g$R 9  
    g$$i WC!S<  
    %%% Calculate the Incident field [*^`rQ  
    for j=1:END-1 H <7r  
        h_inc(j)=h_inc(j)+0.5*(e_inc(j)-e_inc(j+1)); k@vN_U
n  
    end =pSuyM'  
    ,v(K|P@  
    %%% Hx Z&BJ/qk \-  
    for j=1:JE-1 /E8{:>2  
        for k=1:KE-1 HQqnJ;ns<  
            for i=1:ia GIG\bQSv2  
                curl_e=ey(i,j,k+1)-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k); 7L2$(
d4  
                ihxl(i,j,k)=ihxl(i,j,k)+curl_e; BSU%.tmI  
                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)); ^~\cx75D  
            end |-t>_+. J'  
            for i=ia+1:ib Vm\ly;v'R  
                curl_e=ey(i,j,k+1)-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k); |e49F  
                hx(i,j,k)=fj3(j)*fk3(k)*hx(i,j,k)+fj2(j)*fk2(k)*0.5*curl_e; =*<Cw?Gc  
            end 'e @`HG  
            for i=ib+1:IE r: Ij\YQ  
                ixh=i-ib; TrA&yXXL  
                curl_e=ey(i,j,k+1)-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k); t6m&+N  
                ihxh(ixh,j,k)=ihxh(ixh,j,k)+curl_e; ixc~DV+@[  
                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)); K`}8fU  
            end P|c[EUT  
        end www#.D%'U  
    end vw(X9xa  
    ffDh0mDN  
    for i=ia+1:ib _ zM/>Qa  
        for k=ka+1:kb-1 #2AKO/  
            hx(i,ja,k)=hx(i,ja,k)+0.5*ez_left(i,k); CRXIVver  
            hx(i,jb,k)=hx(i,jb,k)-0.5*ez_right(i,k); +$_.${uwV  
        end gT3i{iU  
    end J" :R,w`  
    wNQhz.>y  
    for i=ia+1:ib jFAnhbbCE  
        for j=ja+1:jb-1 OAx5 LTd  
            hx(i,j,ka)=hx(i,j,ka)-0.5*ey_down(i,j); Am>^{qh9  
            hx(i,j,kb)=hx(i,j,kb)+0.5*ey_up(i,j); e\o>(is  
        end MC=pN(l  
    end N_<wiwI<  
        nKxu8YAJe  
    %%% Hy 5MU-Eu|*>  
    for i=1:IE-1 l}\q }7\)  
        for k=1:KE-1 fH{ _X  
            for j=1:ja o!bIaeEaU  
                curl_e=ez(i+1,j,k)-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k); Y'n+,g  
                ihyl(i,j,k)=ihyl(i,j,k)+curl_e; SbrKNADH%  
                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)); =Rv!c+?  
            end l6kqP  
            for j=ja+1:jb 0_\@!#-sml  
                curl_e=ez(i+1,j,k)-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k); Z]VmTB  
                hy(i,j,k)=fi3(i)*fk3(k)*hy(i,j,k)+fi2(i)*fk2(k)*0.5*curl_e; /4K ^-  
            end /O@TqH  
            for j=jb+1:JE &?[uY5Mk  
                jyh=j-jb; {k:W?`  
                curl_e=ez(i+1,j,k)-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k); MLVrL r t  
                ihyh(i,jyh,k)=ihyh(i,jyh,k)+curl_e; :&59N^So|  
                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)); <<#j?%  
            end E+'P|~>oX  
        end v>S[}du  
    end @ ={Hx$zL  
    (tY0/s  
    for i=ia+1:ib-1 xb6y=L  
        for j=ja+1:jb Xx)PyO  
            hy(i,j,ka)=hy(i,j,ka)+0.5*ex_down(i,j); Tw`n3y?  
            hy(i,j,kb)=hy(i,j,kb)-0.5*ex_up(i,j); kF,_o/Jc  
        end |/!3N  
    end ]!% p21e  
    A!x
x#+M  
    for j=ja+1:jb V@%:y tDf  
        for k=ka+1:kb-1
Wycood*  
            hy(ia,j,k)=hy(ia,j,k)-0.5*ez_back(j,k); Xo,}S\wcn  
            hy(ib,j,k)=hy(ib,j,k)+0.5*ez_front(j,k); e#{,M8  
        end Lt|'("($*  
    end ]|6)'L&]*s  
    4n}tDHvd  
   %%% Hz 5}v<?<l9\  
    for i=1:IE-1 TDqH"q0  
        for j=1:JE-1 qhE1 7Hf  
            for k=1:ka 816O
V  
                curl_e=ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k); jPum2U_  
                ihzl(i,j,k)=ihzl(i,j,k)+curl_e; x~;1CB  
                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)); o=%pR|  
            end Uxll<z,  
            for k=ka+1:kb c}$C=s5 h}  
                curl_e=ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k); ()cqax4  
                hz(i,j,k)=fi3(i)*fj3(j)*hz(i,j,k)+fi2(i)*fj2(j)*0.5*curl_e; ]x12_+  
            end {`KRr:w  
            for k=kb+1:KE ^*-6PV#Z  
                kzh=k-kb; J/T$.*X  
                curl_e=ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k); hv4om+  
                ihzh(i,j,kzh)=ihzh(i,j,kzh)+curl_e; o$4xi
nK  
                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)); L'A>IBrz  
            end * |dz.Tr  
        end NA`qC.K   
    end 8Pgw_ 21N1  
    A3su!I2S  
    for i=ia+1:ib-1 3_5]0:?]-  
        for k=ka+1:kb 4ju=5D];  
            hz(i,ja,k)=hz(i,ja,k)-0.5*ex_left(i,k); !=/w
psH  
            hz(i,jb,k)=hz(i,jb,k)+0.5*ex_right(i,k); R}T8cVxc  
        end
=PQMd  
    end N?Nu'  
    k/lU]~PE  
    for j=ja+1:jb-1
Tl/!Dn  
        for k=ka+1:kb agI"Kh]j?  
            hz(ia,j,k)=hz(ia,j,k)+0.5*ey_back(j,k); EklcnM|6  
            hz(ib,j,k)=hz(ib,j,k)-0.5*ey_front(j,k); _{k-&I  
        end d[( 
}  
    end 3=@lJ?Ym  
    -G|?Kl  
    %%%  Visualize fields +xL*`fn  
    %timestep=int2str(N); %k+G-oT5  
    %tview(:,:)=ez(:,:,kc); _Qh:*j!  
    %pcolor(tview); /[<1D|f%  
    %shading interp; eI#b%h  
    %caxis([-1,1]); ,JU3w  
    %colorbar; *IGCFZbp41  
    %axis image; @Q%9b)\\  
    %title(['Ez(i,j,k=kc), time step = ',timestep]); X`s6lV%\  
    %xlabel('i'); p#QR^|7"  
    %ylabel('j'); / %
9DO  
    %pause(0.05);   |mS-<e8LY4  
        !XtZI3Xu  
if N==T-dT Kn#CIFbBN  
    eio=e_inc(is); CW+]Jv]"  
    hio=e_inc(is); VnU/_#n  
    for i=iao:ibo 7kT&}`g.  
        for j=jao:jbo            _G62E$=  
            exi_oup(i,j)=ex(i,j,kbo); g##yR/L  
            eyi_oup(i,j)=ey(i,j,kbo); v9_7OMl/x  
            &%=]lP]  
            exi_odown(i,j)=ex(i,j,kao); *y
dh.R<hb  
            eyi_odown(i,j)=ey(i,j,kao); GxynLXWo>  
            2d60o~E  
            hxi_oup(i,j)=0.5*(hx(i,j,kbo)+hx(i,j,kbo-1)); l.iT+T  
            hyi_oup(i,j)=0.5*(hy(i,j,kbo)+hy(i,j,kbo-1)); B3';Tcs  
            zn3]vU!  
            hxi_odown(i,j)=0.5*(hx(i,j,kao)+hx(i,j,kao-1)); nD
5+&M0  
            hyi_odown(i,j)=0.5*(hy(i,j,kao)+hy(i,j,kao-1)); 9!>Ks8'.d  
        end ,qz:(Nr  
    end <Na .6P  
    for i=iao:ibo uH.1'bR?a  
        for k=kao:kbo @ =XJ<  
            exi_oleft(i,k)=ex(i,jao,k); P/ XO5`  
            ezi_oleft(i,k)=ez(i,jao,k); Z!81\5  
            <0h,{28  
            exi_oright(i,k)=ex(i,jbo,k); K'@lXA:  
            ezi_oright(i,k)=ez(i,jbo,k); <_8p6{=  
            lK2=[%,~  
            hxi_oleft(i,k)=0.5*(hx(i,jao,k)+hx(i,jao-1,k)); r:~q{  
            hzi_oleft(i,k)=0.5*(hz(i,jao,k)+hz(i,jao-1,k)); J=`2{ 'l  
            .Xcf*$.;s  
            hxi_oright(i,k)=0.5*(hx(i,jbo,k)+hx(i,jbo-1,k)); :T9 P9<  
            hzi_oright(i,k)=0.5*(hz(i,jbo,k)+hz(i,jbo-1,k)); `P43O gA  
        end Hgk@I;  
    end 1Z_2s2`p  
    for j=jao:jbo ;i>(r;ZM  
        for k=kao:kbo `z!6zo2d  
            eyi_ofront(j,k)=ey(ibo,j,k); {k15!(:i~a  
            ezi_ofront(j,k)=ez(ibo,j,k); x^XP<R{D  
            ~:xR0dqx  
            eyi_oback(j,k)=ey(iao,j,k); 04E S>'@  
            ezi_oback(j,k)=ez(iao,j,k); C&/_mm
5  
            $p)e.ZMgE  
            hyi_ofront(j,k)=0.5*(hy(ibo,j,k)+hy(ibo-1,j,k)); oa"_5kn,  
            hzi_ofront(j,k)=0.5*(hz(ibo,j,k)+hz(ibo-1,j,k)); <
t"KNKI  
            ny'wS  
            hyi_oback(j,k)=0.5*(hy(iao,j,k)+hy(iao-1,j,k)); mr[+\ 5  
            hzi_oback(j,k)=0.5*(hz(iao,j,k)+hz(iao-1,j,k)); N2
Qb+  
        end t4H*&U  
    end ?`+G0
VT  
end PFSh_9.q  
if N==T TOge!Q>a  
    ero=e_inc(is); G5T(  
    hro=h_inc(is); p?H2W-  
    for i=iao:ibo XlaGR2-%  
        for j=jao:jbo            &VdKL2  
            exr_oup(i,j)=ex(i,j,kbo); wpK[;  
            eyr_oup(i,j)=ey(i,j,kbo); i[O{M`Z%  
            /V"6Q'D  
            exr_odown(i,j)=ex(i,j,kao); "IA:,j.#g  
            eyr_odown(i,j)=ey(i,j,kao); }J0HEpn4  
            Wi"3kps q  
            hxr_oup(i,j)=0.5*(hx(i,j,kbo)+hx(i,j,kbo-1)); z0-[ RGg  
            hyr_oup(i,j)=0.5*(hy(i,j,kbo)+hy(i,j,kbo-1)); !;U;5e=0  
            " Z2D@l  
            hxr_odown(i,j)=0.5*(hx(i,j,kao)+hx(i,j,kao-1)); O
BEHUJ5  
            hyr_odown(i,j)=0.5*(hy(i,j,kao)+hy(i,j,kao-1)); y8Oz4|  
        end .:(T}\]R  
    end r=4vN=:  
    for i=iao:ibo sD{d8s[(  
        for k=kao:kbo 6BY(Y(z  
            exr_oleft(i,k)=ex(i,jao,k); 9.^2CM6l  
            ezr_oleft(i,k)=ez(i,jao,k); HU47S  
            (p!w`MSv  
            exr_oright(i,k)=ex(i,jbo,k); mhv ;pM6  
            ezr_oright(i,k)=ez(i,jbo,k); mrGfu:r  
            XbYST%|.  
            hxr_oleft(i,k)=0.5*(hx(i,jao,k)+hx(i,jao-1,k)); F['%?+<3  
            hzr_oleft(i,k)=0.5*(hz(i,jao,k)+hz(i,jao-1,k)); h}n?4B~Gi  
            M=t;t0  
            hxr_oright(i,k)=0.5*(hx(i,jbo,k)+hx(i,jbo-1,k)); r-'CB  
            hzr_oright(i,k)=0.5*(hz(i,jbo,k)+hz(i,jbo-1,k)); GWCU9n  
        end }c@duf-l  
    end dUc([&  
    for j=jao:jbo gaC[%M  
        for k=kao:kbo iu.$P-s  
            eyr_ofront(j,k)=ey(ibo,j,k); #8CeTR23cw  
            ezr_ofront(j,k)=ez(ibo,j,k); go!jx6~;x  
            uMb[0-5  
            eyr_oback(j,k)=ey(iao,j,k); _(J4  
            ezr_oback(j,k)=ez(iao,j,k); pdq5EUdS  
            3ahbv%y  
            hyr_ofront(j,k)=0.5*(hy(ibo,j,k)+hy(ibo-1,j,k)); /`+ubF
Xc  
            hzr_ofront(j,k)=0.5*(hz(ibo,j,k)+hz(ibo-1,j,k)); l",X  
            *Q!I^]CR  
            hyr_oback(j,k)=0.5*(hy(iao,j,k)+hy(iao-1,j,k)); t+t
D  
            hzr_oback(j,k)=0.5*(hz(iao,j,k)+hz(iao-1,j,k)); Z1N=tL  
        end Prqr,  
    end Rl[SqmnI)@  
end )n{9*{Ch  
    X ApSKJ  
end ^EZ?wdL  
eo=ero+sqrt(-1)*eio; bS&XlgnKi  
eo=abs(eo); 
{D`_q
|  
ho=(hro+sqrt(-1)*hio)*exp(-sqrt(-1)*0.5*omiga); 6|B;C  
for i=iao:ibo 3,dIW*<**  
    for j=jao:jbo =`/GBT$  
        ex_oup(i,j)=(exr_oup(i,j)+sqrt(-1)*exi_oup(i,j))/eo; fb_q2p} G  
        ey_oup(i,j)=(eyr_oup(i,j)+sqrt(-1)*eyi_oup(i,j))/eo; qE@H~&  
        v&3 Oc  
        ex_odown(i,j)=(exr_odown(i,j)+sqrt(-1)*exi_odown(i,j))/eo; L 4Sa,ZL  
        ey_odown(i,j)=(eyr_odown(i,j)+sqrt(-1)*eyi_odown(i,j))/eo; V[^+lR  
        W+'f|J=  
        hx_oup(i,j)=(hxr_oup(i,j)+sqrt(-1)*hxi_oup(i,j))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; Y5R|)x  
        hy_oup(i,j)=(hyr_oup(i,j)+sqrt(-1)*hyi_oup(i,j))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; rvG0aqO`  
        +T8]R7b9  
        hx_odown(i,j)=(hxr_odown(i,j)+sqrt(-1)*hxi_odown(i,j))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; #%{x*y:Ms  
        hy_odown(i,j)=(hyr_odown(i,j)+sqrt(-1)*hyi_odown(i,j))*exp(-sqrt(-1)*0.5*omiga)/eo/Z;
/=gU  
    end * am
Z  
end 4#&w-W  
for i=iao:ibo abtYa  
    for k=kao:kbo LLU]KZhtY|  
        ex_oleft(i,k)=(exr_oleft(i,k)+sqrt(-1)*exi_oleft(i,k))/eo; }_,={<g  
        ez_oleft(i,k)=(ezr_oleft(i,k)+sqrt(-1)*ezi_oleft(i,k))/eo; 6SP!J*F  
        .-?Txkwb  
        ex_oright(i,k)=(exr_oright(i,k)+sqrt(-1)*exi_oright(i,k))/eo; DP08$Iq  
        ez_oright(i,k)=(ezr_oright(i,k)+sqrt(-1)*ezi_oright(i,k))/eo; 4O>0gK{w  
        *f k3IvAXu  
        hx_oleft(i,k)=(hxr_oleft(i,k)+sqrt(-1)*hxi_oleft(i,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; txFcV  
        hz_oleft(i,k)=(hzr_oleft(i,k)+sqrt(-1)*hzi_oleft(i,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; &2//\Qz  
        0RGqpJxk  
        hx_oright(i,k)=(hxr_oright(i,k)+sqrt(-1)*hxi_oright(i,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; FKU)# Eo  
        hz_oright(i,k)=(hzr_oright(i,k)+sqrt(-1)*hzi_oright(i,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; m*rw?nLZ  
    end A(T=  
end U`kO<ztk  
for j=jao:jbo gI{56Z  
    for k=kao:kbo CJf4b:SY@  
        ey_ofront(j,k)=(eyr_ofront(j,k)+sqrt(-1)*eyi_ofront(j,k))/eo; +O?`uV  
        ez_ofront(j,k)=(ezr_ofront(j,k)+sqrt(-1)*ezi_ofront(j,k))/eo; KF}_|~~T  
        Gi Max  
        ey_oback(j,k)=(eyr_oback(j,k)+sqrt(-1)*eyi_oback(j,k))/eo; aSH =|Jnc  
        ez_oback(j,k)=(ezr_oback(j,k)+sqrt(-1)*ezi_oback(j,k))/eo; <qjolMO`  
        -5u. Ix3  
        hy_ofront(j,k)=(hyr_ofront(j,k)+sqrt(-1)*hyi_ofront(j,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; #|^yWw^  
        hz_ofront(j,k)=(hzr_ofront(j,k)+sqrt(-1)*hzi_ofront(j,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; 3bZ:*6W.6  
        Hq0O!Zv  
        hy_oback(j,k)=(hyr_oback(j,k)+sqrt(-1)*hyi_oback(j,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; _64<[2  
        hz_oback(j,k)=(hzr_oback(j,k)+sqrt(-1)*hzi_oback(j,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; ds:&{~7L<T  
    end UdK+,k~m/  
end 3}x6IM2  
%%%%%%%%%  Far Field Calculate  %%%%%%% I7?s+vyds  
rk=2*pi/WL; k5
6*eE
c  
fx=zeros(1,181); t6! B  
fy=zeros(1,181); Y8Z-m (OQ  
fz=zeros(1,181); B/gI~e0  
fmx=zeros(1,181); F,h}HlU  
fmy=zeros(1,181); ^.u J]k0  
fmz=zeros(1,181); %Zi}sm1t  
rcs=zeros(1,181); e?\34F  
est=zeros(1,181); )vy_m_f&  
efa=zeros(1,181); r`A|2(h5B  
fai_s=0*pi/180; ;s!GpO7+  
cps2=cos(fai_s); N51WY7  
sps2=sin(fai_s); N7#,x9+E  
ss=0; f@
H>by N  
for sita_s=0:pi/180:pi n{E9p3
i  
    ss=ss+1; :,pdR>q%(y  
    sps1=sin(sita_s); }ZkGH}K_}  
    cps1=cos(sita_s); F0tx.]uS  
    rk=2*pi/WL; Ie'iAY  
    %%% fx、fmx——left/right output boundry )O$T; U  
    for i=iao+1:ibo ~NZ}@J{00_  
        for k=kao:kbo 0khAi|PY  
            temp1=exp(sqrt(-1)*rk*((i-0.5)*sps1*cps2+jao*sps1*sps2+k*cps1));  Dac ,yW  
            fx(ss)=fx(ss)-(1-0.5*(k==kao)|(k==kbo))*ddx*ddx*hz_oleft(i,k)*temp1; >+F +"NAN  
            fmx(ss)=fmx(ss)+(1-0.5*(k==kao)|(k==kbo))*ddx*ddx*ez_oleft(i,k)*temp1; 8 _|"+Ze  
            $x;(C[  
            temp2=exp(sqrt(-1)*rk*((i-0.5)*sps1*cps2+jbo*sps1*sps2+k*cps1)); G4QsR7  
            fx(ss)=fx(ss)+(1-0.5*(k==kao)|(k==kbo))*ddx*ddx*hz_ori .. 1pZ[rM'}  
j38>5DM6L  

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

大段的代码肯定难看，你可以尝试以下改正方法： k=JT%  
1，低级失误，i,j,k混淆,x,y,z混淆 `]m/za%7  
2，看动态显示图，这是最有效的方法，哪里出问题直接定位到处理这部分网格的程序中。 }I ^e:,{  
3，运行一次，观察MATLAB中的各变量，按经验判断错误的计算值

好的，谢谢！ k7ODQ(*v  
l#G }j^Q  
根据你的提醒调试，希望能找出错误。 #3o]Qo[Sc