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

       下面是我整的三维时谐场——用FDTD计算介质球的RCS的程序，但是在近远场外推的时候出现了问题，怎么都跑不出图来， `78:TU~5S  
希望大家能帮我找出错误，谢谢！ SZ[,(h  
       IE=81; Fs,#d%4@%  
JE=81; qz .{[l  
KE=81; mlmp'f  
ia=9; ^k7`:@ z0U  
ja=9; wZ_k]{J  
ka=9; Vo2{aK;  
ic=(IE+1)/2; [9S
?  
jc=(JE+1)/2; o2C{V1nB  
kc=(KE+1)/2; 1\'zq;I~  
ib=IE-ia; hH|m
oj]  
jb=JE-ja; O!7v&$]1  
kb=KE-ka; e MT5bn  
epsz=8.8e-12; AQH\ ;L  
muz=4*pi*1e-7; D[5Qd)PIL  
lamda=0.032; r(#]Z  
freq=lamda/3e+8; x&SG gl
 
%freq=30e+7; *$eMM*4  
ddx=lamda/42; ogrh"  
dt=ddx/6e+8; #o]/&T=N=  
arg=2*pi*freq*dt; YAvOV-L  
dx=zeros(IE,JE,KE); RZ#~^5DiO  
dy=zeros(IE,JE,KE);  @{|vW  
dz=zeros(IE,JE,KE); xy$agt>j>  
ex=zeros(IE,JE,KE); ^mA^7jB  
ey=zeros(IE,JE,KE); quPNwNy  
ez=zeros(IE,JE,KE); (` N@4w=  
hx=zeros(IE,JE,KE); ~T_4M  
hy=zeros(IE,JE,KE); cU^Z=B  
hz=zeros(IE,JE,KE); muc>4!Q  
ix=zeros(IE,JE,KE); l6wN&JHTh  
iy=zeros(IE,JE,KE);  XAb!hc   
iz=zeros(IE,JE,KE); cn~M:LW23  
gax=zeros(IE,JE,KE); zFn-VEJ)  
gay=zeros(IE,JE,KE); Xj^Hy"HC^~  
gaz=zeros(IE,JE,KE); Hmi]qK[F  
gbx=zeros(IE,JE,KE); tXgsWG?v[H  
gby=zeros(IE,JE,KE); @6N$!Q?  
gbz=zeros(IE,JE,KE); O_u2V'jy9  
END=200; bvK fxAih  
e_inc=zeros(1,END); y)B>g/Hoh  
h_inc=zeros(1,END); TC%ENxDR  
e_temp=zeros(1,4); e{"r3*  
idxl=zeros(ia,JE,KE); LR5X=&k  
idxh=zeros(ia,JE,KE); o'8`>rb  
ihxl=zeros(ia,JE,KE); 6/Pw'4H9$  
ihxh=zeros(ia,JE,KE); <^APq8>  
idyl=zeros(IE,ja,KE); "?<$>\@; q  
idyh=zeros(IE,ja,KE); CQ`$' oy?W  
ihyl=zeros(IE,ja,KE); wInJ!1  
ihyh=zeros(IE,ja,KE); G{ 9p.Q  
idzl=zeros(IE,JE,ka); t%`
GXJb  
idzh=zeros(IE,JE,ka); :Rq>a@Rp  
ihzl=zeros(IE,JE,ka); =/J{>S>(i  
ihzh=zeros(IE,JE,ka); m1mA:R\zM  
gi1=zeros(1,IE);gi2=ones(1,IE);gi3=ones(1,IE); ,o3{?o]s  
gj1=zeros(1,JE);gj2=ones(1,JE);gj3=ones(1,JE); KG!W,tB  
gk1=zeros(1,KE);gk2=ones(1,KE);gk3=ones(1,KE); 0"*!0s~  
fi1=zeros(1,IE);fi2=ones(1,IE);fi3=ones(1,IE); iIe\mV  
fj1=zeros(1,JE);fj2=ones(1,JE);fj3=ones(1,JE); v
e [*t`  
fk1=zeros(1,KE);fk2=ones(1,KE);fk3=ones(1,KE); f""+jc1  
sita=0*pi/180;    %pointing angle NR*s7>  
fai=0*pi/180;     %azimuth angle tS2Orzc>,  
arfa=0*pi/180;    %polarizing angle j~IX  
sp1=sin(sita);cp1=cos(sita);  9z9EK'g  
sp2=sin(fai);cp2=cos(fai); aWe?n;  
sp3=sin(arfa);cp3=cos(arfa); 0KQDw  
ex_left=zeros(IE,KE);ex_right=zeros(IE,KE); rX-V0  
ex_up=zeros(IE,JE);ex_down=zeros(IE,JE); yv@td+-"D  
ey_front=zeros(JE,KE);ey_back=zeros(JE,KE); -f+U:/'.>v  
ey_up=zeros(IE,JE);ey_down=zeros(IE,JE); d1"%sI  
ez_left=zeros(IE,KE);ez_right=zeros(IE,KE); &ZmHR^Flz  
ez_front=zeros(JE,KE);ez_back=zeros(JE,KE); t=IpVl!  
hx_left=zeros(IE,KE);hx_right=zeros(IE,KE); oY#62&wk4  
hx_up=zeros(IE,JE);hx_down=zeros(IE,JE); U49#?^?  
hy_front=zeros(JE,KE);hy_back=zeros(JE,KE); xTy[X"sJ  
hy_up=zeros(IE,JE);hy_down=zeros(IE,JE); L1'#wH  
hz_left=zeros(IE,KE);hz_right=zeros(IE,KE); HFr#Ql>g  
hz_front=zeros(JE,KE);hz_back=zeros(JE,KE); m,,FNYW  
npml=4;        %pml层的厚度 @.PVUP  
xxn=1:-1/npml:1/npml; `H:5D5]  
xn=0.3333*xxn.^3; )_vE"ryThA  
fi1(1:npml)=xn; Z uh!{_x;  
fi1(IE:-1:IE-npml+1)=xn; 'DB4po.   
gi2(1:npml)=1./(1+xn); `H_.<``>  
gi2(IE:-1:IE-npml+1)=1./(1+xn); phT|w H  
gi3(1:npml)=(1-xn)./(1+xn); v7i5R !  
gi3(IE:-1:IE-npml+1)=(1-xn)./(1+xn); ? ^EB"{  
fj1(1:npml)=xn; ^,?dk![1Cv  
fj1(JE:-1:JE-npml+1)=xn; km)5?  
gj2(1:npml)=1./(1+xn); 1Rrl59}5  
gj2(JE:-1:JE-npml+1)=1./(1+xn); k|Hxd^^I  
gj3(1:npml)=(1-xn)./(1+xn); \sUk71L`j  
gj3(JE:-1:JE-npml+1)=(1-xn)./(1+xn); -t<8)9q(  
fk1(1:npml)=xn; OJkiTs{  
fk1(KE:-1:KE-npml+1)=xn; )@X `B d  
gk2(1:npml)=1./(1+xn); VeJM=s.y7  
gk2(KE:-1:KE-npml+1)=1./(1+xn); %(
c5T
)B9  
gk3(1:npml)=(1-xn)./(1+xn); 80p?qe  
gk3(KE:-1:KE-npml+1)=(1-xn)./(1+xn); LV
^V`m0#  
xxn=(1-0.5/npml):-1/npml:0.5/npml; eC9nOwp]xH  
xn=0.3333*xxn.^3; !OVTs3}  
gi1(1:npml)=xn; 8op,;Z7Y  
gi1(IE-1:-1:IE-npml)=xn; M7\; Y  
fi2(1:npml)=1./(1+xn); j"8f,er  
fi2(IE-1:-1:IE-npml)=1./(1+xn); cVg!"  
fi3(1:npml)=(1-xn)./(1+xn); `Y
ZK$ -,  
fi3(IE-1:-1:IE-npml)=(1-xn)./(1+xn); J^Dkx"1GD  
gj1(1:npml)=xn; `qNhB\  
gj1(JE-1:-1:JE-npml)=xn; |\1!*Qp  
fj2(1:npml)=1./(1+xn); dKOW5\H'  
fj2(JE-1:-1:JE-npml)=1./(1+xn); F|eKt/>e  
fj3(1:npml)=(1-xn)./(1+xn); @+9x8*~S'  
fj3(JE-1:-1:JE-npml)=(1-xn)./(1+xn); #~ :j< =o  
gk1(1:npml)=xn; &I&:  
gk1(KE-1:-1:KE-npml)=xn; 63J_u-o  
fk2(1:npml)=1./(1+xn); t3~ZGOn  
fk2(KE-1:-1:KE-npml)=1./(1+xn);  KKfC^g  
fk3(1:npml)=(1-xn)./(1+xn); )CYm
/dk  
fk3(KE-1:-1:KE-npml)=(1-xn)./(1+xn); 4
4uM:;  
WL=42; k,~I>qg  
dT=WL/2; Auq)  
omiga=pi/WL; YPjjSi:#  
Z=sqrt(muz/epsz); @+gr>a1K#  
odist=2; wTpjM@F?J|  
iao=ia+1-odist; 1YGj^7V)|Z  
jao=ja+1-odist; ] Q 'Ed  
kao=ka+1-odist; n^svRM]eQ  
ibo=ib+odist; ^x! N]  
jbo=jb+odist; y5td o'Ex  
kbo=kb+odist; /_N*6a
~  
exi_oup=zeros(IE,JE);exr_oup=zeros(IE,JE);ex_oup=zeros(IE,JE); #s'UA!)  
exi_odown=zeros(IE,JE);exr_odown=zeros(IE,JE);ex_odown=zeros(IE,JE); k63]Qf=5?N  
eyi_oup=zeros(IE,JE);eyr_oup=zeros(IE,JE);ey_oup=zeros(IE,JE); (5^bU<  
eyi_odown=zeros(IE,JE);eyr_odown=zeros(IE,JE);ey_odown=zeros(IE,JE); '^TQ Ubw  
exi_oleft=zeros(IE,KE);exr_oleft=zeros(IE,KE);ex_oleft=zeros(IE,KE); EXdx$I=X  
ezi_oleft=zeros(IE,KE);ezr_oleft=zeros(IE,KE);ez_oleft=zeros(IE,KE); G lz0`z  
exi_oright=zeros(IE,KE);exr_oright=zeros(IE,KE);ex_oright=zeros(IE,KE); lw.4O^  
ezi_oright=zeros(IE,KE);ezr_oright=zeros(IE,KE);ez_oright=zeros(IE,KE); Po%+:0oX  
eyi_ofront=zeros(JE,KE);eyr_ofront=zeros(JE,KE);ey_ofront=zeros(JE,KE); f{b$Y3  
ezi_ofront=zeros(JE,KE);ezr_ofront=zeros(JE,KE);ez_ofront=zeros(JE,KE); nX@lR~g%F  
eyi_oback=zeros(JE,KE);eyr_oback=zeros(JE,KE);ey_oback=zeros(JE,KE); our$Ka31  
ezi_oback=zeros(JE,KE);ezr_oback=zeros(JE,KE);ez_oback=zeros(JE,KE); A]z~Dw3  
hxi_oup=zeros(IE,JE);hxr_oup=zeros(IE,JE);hx_oup=zeros(IE,JE); aR iD}P*V  
hxi_odown=zeros(IE,JE);hxr_odown=zeros(IE,JE);hx_odown=zeros(IE,JE); DNP%]{J  
hyi_oup=zeros(IE,JE);hyr_oup=zeros(IE,JE);hy_oup=zeros(IE,JE); *1}UK9X;  
hyi_odown=zeros(IE,JE);hyr_odown=zeros(IE,JE);hy_odown=zeros(IE,JE); PRs[!EB6  
hxi_oleft=zeros(IE,KE);hxr_oleft=zeros(IE,KE);hx_oleft=zeros(IE,KE); ST#OO!  
hzi_oleft=zeros(IE,KE);hzr_oleft=zeros(IE,KE);hz_oleft=zeros(IE,KE); %s+H& vfQs  
hxi_oright=zeros(IE,KE);hxr_oright=zeros(IE,KE);hx_oright=zeros(IE,KE); "kLu]M<  
hzi_oright=zeros(IE,KE);hzr_oright=zeros(IE,KE);hz_oright=zeros(IE,KE); ileqI/40f  
hyi_ofront=zeros(JE,KE);hyr_ofront=zeros(JE,KE);hy_ofront=zeros(JE,KE); MOiTzL*  
hzi_ofront=zeros(JE,KE);hzr_ofront=zeros(JE,KE);hz_ofront=zeros(JE,KE); MH"{N "|  
hyi_oback=zeros(JE,KE);hyr_oback=zeros(JE,KE);hy_oback=zeros(JE,KE); o3_dHbdI  
hzi_oback=zeros(JE,KE);hzr_oback=zeros(JE,KE);hz_oback=zeros(JE,KE); (&25 8i,  
%%%介质球的设置 duCso M/  
epsilon=[1 4]; (8[etm  
%epsilon=[1 81]; 8{m5P8w'  
sigma=[0 0]; WBo|0(#  
%sigma=[0 5]; FC'v= *  
%sigma=[0 0]; X3z$f(lF%)  
radius=20; L(p{>Ykcc  
for i=ia+1:ib M^i^_}~S;  
    for j=ja+1:jb j89C~xP6  
        for k=ka+1:kb h";0i:  
            eps=epsilon(1);  H %Cb  
            cond=sigma(1); l3d^V&Sk  
            xdist=(ic-i-0.5); `^rN"\  
            ydist=jc-j; 7e{w)m:A  
            zdist=kc-k; e|`QW|9 .  
            dist=sqrt(xdist^2+ydist^2+zdist^2); )1PZ#  
            if dist<=radius %gF; A*  
                eps=epsilon(2); Km5#$IiP;  
                cond=sigma(2); E>/kNl  
            end c$cb2V7,  
            gax(i,j,k)=1/(eps+cond*dt/epsz); $VWeo#b  
            gbx(i,j,k)=cond*dt/epsz; WUVRwJ 5  
        end 20 j9~+  
    end QKj-"y[  
end :YL`GSl  
for i=ia+1:ib [k"@n+%  
    for j=ja+1:jb 9t 3mU:  
        for k=ka+1:kb >dnH  
            eps=epsilon(1); RXM
zwk  
            cond=sigma(1); ?w{lC,  
            xdist=ic-i; E]^wsS>=  
            ydist=jc-j-0.5; :# 1d;jx  
            zdist=kc-k; +=XDNSw  
            dist=sqrt(xdist^2+ydist^2+zdist^2); ?2Q9z-$  
            if dist<=radius HFJna2B`  
                eps=epsilon(2); PT9,R^2T!  
                cond=sigma(2); QB<9Be@e  
            end (+@ Lnz\  
            gay(i,j,k)=1/(eps+cond*dt/epsz); 6Un61s  
            gby(i,j,k)=cond*dt/epsz; rf2+~B{$,  
        end !xU1[,9  
    end mSn>  
end y"|QY!fK  
for i=ia+1:ib 6x1!!X+)+  
    for j=ja+1:jb $=7'Cm
?  
        for k=ka+1:kb ~0m
O<0~  
            eps=epsilon(1); {Vc%ga|E  
            cond=sigma(1); @yBg)1AL  
            xdist=ic-i; -P;_j,~U  
            ydist=jc-j; F %OA  
            zdist=kc-k-0.5; *Q?
ZJS~  
            dist=sqrt(xdist^2+ydist^2+zdist^2); V3<baxdE  
            if dist<=radius o"O=Epg  
                eps=epsilon(2); c:  /Wk  
                cond=sigma(2); N5 BC<pu  
            end K`=O!;  
            gaz(i,j,k)=1/(eps+cond*dt/epsz); 9uV'#sR  
            gbz(i,j,k)=cond*dt/epsz; '=|2, H]  
        end zJp}
JO  
    end ubC(%Y_k  
end 8PQn=k9  
real_in=0; {-(}p+;z  
imag_in=0; @a AR99M  
%t0=40; BPp`r_m8w}  
%spread=10; 8>Y  
T=1000; /Iwnl   
%N=0; |5uvmK  
%isource=floor(sqrt((ib-ia)^2+(jb-ja)^2+(kb-ka)^2)); d$>TC(E=t  
isource=4; kB|jN~  
is=ja+1; pYtG%<  
for N=1:T vvJ{fi  
    d(.e%[`  
    for j=2:END UL81x72O  
        e_inc(j)=e_inc(j)+0.5*(h_inc(j-1)-h_inc(j)); mv7><C  
    end zEG6T*  
    D&Xh|}2A  
    %real_in=real_in+cos(arg*T)*e_inc(ja-1); \a:#e%]qz9  
    %imag_in=imag_in-sin(arg*T)*e_inc(ja-1); %SKp<>;9  
    ,o$F~KPu  
    %pulse=exp(-0.5*((t0-T)/spread)^2); (=v :@\r  
    pulse=sin(arg*N); L5%t.7B  
    e_inc(isource)=pulse; XcOfQs  
    P8tpbdZE-  
    e_inc(1)=e_temp(1); "}_b,5lkGK  
    e_temp(1)=e_temp(2); Eei"baw/  
    e_temp(2)=e_inc(2); o:H^ L,<Tl  
    e_inc(END)=e_temp(3); hdL/zW7]  
    e_temp(3)=e_temp(4); <qG4[W,
[  
    e_temp(4)=e_inc(END-1);
)E--E+j  
    C- Aiv@@<=  
    %%%  Calculate Incident component on connected edge pwg$% lv  
    %%%  up/down connect edge .A;e`cKb  
    for i=ia+1:ib [>5<&[A  
        for j=ja+1:jb d|GQZAEJEt  
            k=kb; =x9SvIm/tH  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; c/l%:!A  
            w=dr-floor(dr); dP>~ExYtm  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); ~H[  
            ex_up(i,j)=ein*(sp2*cp3-cp1*cp2*sp3); <eU1E}BDQ  
            /}G+PUk7  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1;  %2 A-u  
            w=dr-floor(dr); K>=KsG  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); a;=)`  
            ey_up(i,j)=ein*(-cp2*cp3-cp1*sp2*sp3);  yN9k-IPI  
            z/IA @  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; i/ED_<_Vg  
            w=dr-floor(dr); CqMm'6;$a}  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); {!? @u?M  
            hx_up(i,j)=hin*(sp2*sp3+cp1*cp2*cp3); )tH.P: 1~,  
            ~U] "
dbQ  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; R+'$V$g\X  
            w=dr-floor(dr); :}QBrd  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); F`/-Q>Q  
            hy_up(i,j)=hin*(-cp2*sp3+cp1*sp2*cp3); 3\x
@G
)1  
        end j]-0m4QF  
    end grvm2`u  
    ]V|rOtxb  
    for i=ia+1:ib +D[

|Mi  
        for j=ja+1:jb qPh @Bl3  
            k=ka+1; S6k R o^2  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; -^ )0c  
            w=dr-floor(dr); W'aZw9  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); gZa/?[+  
            ex_down(i,j)=ein*(sp2*cp3-cp1*cp2*sp3); Yt]tRqrh;T  
            gN2$;hb?  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; D?+\"lI  
            w=dr-floor(dr); 8nt3Sm  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); |Z]KF>S]  
            ey_down(i,j)=ein*(-cp2*cp3-cp1*sp2*sp3); r5

7&F`{  
            #nK38W#  
            k=ka; $;kFuJF  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; ~0}gRpMW  
            w=dr-floor(dr); "Di27Rq  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); ;[-OMGr]#  
            hx_down(i,j)=hin*(sp2*sp3+cp1*cp2*cp3); c$Vu/dgx  
            1J`<'{*  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; nz{ ;]U1  
            w=dr-floor(dr); l@;UwnI  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); 9YpgzCx Z  
            hy_down(i,j)=hin*(-cp2*sp3+cp1*sp2*cp3); ;kSRv=S  
        end nkI+"$Rz0  
    end bsfYz  
    `Aa}q(}k  
    %%%  left/right connect edge }tq  
    for i=ia+1:ib 9Fo00"q  
        for k=ka+1:kb `pMI@"m  
            j=ja+1; #Tc]L<."
 
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; r[doN{%  
            w=dr-floor(dr); $8/=@E{51  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); BZS%p  
            ex_left(i,k)=ein*(sp2*cp3-cp1*cp2*sp3); *QKxrg  
            oXC|q-(C
 
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; ]><K8N3Z  
            w=dr-floor(dr); n$2IaE;v  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); 8Zj=:;  
            ez_left(i,k)=ein*(sp1*sp3); hv)>HU
&  
            W''%{A/'  
            j=ja; %b
u$t,  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; AcZ{B<  
            w=dr-floor(dr); eAU0 8gM.  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); vY 0EffZ  
            hx_left(i,k)=hin*(sp2*sp3+cp1*cp2*cp3); ew13qpt)<L  
            x)35}mi){L  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; `IP?w&k)  
            w=dr-floor(dr); \a<7DTV  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); jL9g.q4^  
            hz_left(i,k)=hin*(-sp1*cp3); YmrrZ&]q  
        end JD`IPQb~E  
    end L/ L#[  
       TLcev*  
    for i=ia+1:ib &a;{ed1B  
        for k=ka+1:kb R07]{  
            j=jb;  Dno]N  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; NCrNlHIF  
            w=dr-floor(dr); 6?;U[eV  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); BiFU3FlTf  
            ex_right(i,k)=ein*(sp2*cp3-cp1*cp2*sp3); _Y
{8FN(4  
            UL{
+
mp  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; _f5>r(1Q  
            w=dr-floor(dr); OD@k9I[  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); y4\(ynk  
            ez_right(i,k)=ein*(sp1*sp3); Tu!2lHK;  
            OC?a[^hB^)  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; QN4{xf:}S  
            w=dr-floor(dr); *B4?(&0  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); Vy.gr4Cm  
            hx_right(i,k)=hin*(sp2*sp3+cp1*cp2*cp3); \ltbiDP2  
            fL^$G;_?3  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; Z\YCjs%  
            w=dr-floor(dr); -_*ux!  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); {GH0>
1&  
            hz_right(i,k)=hin*(-sp1*cp3); m|') A  
        end ^}w@&Bje  
    end X!aC6gujOH  
    u%t/W0xi  
    %%%  front/back connect edge O]-)?y/  
    for j=ja+1:jb Yvi.l6JL  
        for k=ka+1:kb 'hoEdJ]t5  
            i=ib; tPp9=e2[s  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; ojc m%yd  
            w=dr-floor(dr); l|kGp~  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); OLH[F  
            ey_front(j,k)=ein*(-cp2*cp3-cp1*sp2*sp3); N8[ &1  
            v}cTS@0  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; qjBF]3%t%  
            w=dr-floor(dr); c-jE1y<  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); RJ\'"XQ  
            ez_front(j,k)=ein*(sp1*sp3); Sg&0a$  
            V=i/cI\  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; 0~a9gBG  
            w=dr-floor(dr); ugu|?z*dI  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1);  ff;9P5X  
            hy_front(j,k)=hin*(-cp2*sp3+cp1*sp2*cp3); Ub,5~I+`  
            9QXBz=Fnf  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; >='y+68  
            w=dr-floor(dr); ut*sx9
l  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); MyZ5~jnr\  
            hz_front(j,k)=hin*(-sp1*cp3); e2"<3  
        end <)68ol~<  
    end +$uQ_ve  
    A!^ d8#~.  
    for j=ja+1:jb @u>:(9bp  
        for k=ka+1:kb /f*QxNZ,p  
            i=ia+1; Z|#G+$"QV  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; nF5\iV  
            w=dr-floor(dr); ;aj4V<@  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); 7Dt*++:  
            ey_back(j,k)=ein*(-cp2*cp3-cp1*sp2*sp3); qWdob>u  
            4)Bk:K  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; g[c_rty  
            w=dr-floor(dr); J16t&Ha`  
            ein=(1-w)*e_inc(is+floor(dr))+w*e_inc(is+floor(dr)+1); h5o6G1ur  
            ez_back(j,k)=ein*(sp1*sp3); $etw'c0  
            yI{4h $c  
            i=ia;  *Cj<Vy  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; kEYkd@{  
            w=dr-floor(dr); rq#\x{l  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); ;f!}vo<;  
            hy_back(j,k)=hin*(-cp2*sp3+cp1*sp2*cp3); "C]v  
            iW?z2%#  
            dr=(i-ia-1)*sp1*cp2+(j-ja-1)*sp1*sp2+(k-ka-1)*cp1; On~
w`  
            w=dr-floor(dr); eqY8;/  
            hin=(1-w)*h_inc(is+floor(dr))+w*h_inc(is+floor(dr)+1); gcA,u)z}R  
            hz_back(j,k)=hin*(-sp1*cp3); e"d-$$'e  
        end Nx;Oz  
    end T5_/*`F  
    @e#{Sm  
    %%% Dx I3.cy i  
    for j=2:JE \H4$9lPk  
        for k=2:KE p 2>
\  
            for i=2:ia WvoJ^{\4N*  
                curl_h=hz(i,j,k)-hz(i,j-1,k)-hy(i,j,k)+hy(i,j,k-1); TpGnSD  
                idxl(i,j,k)=idxl(i,j,k)+curl_h;  NZu2D  
                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)); O>@ChQF  
            end Q,LDn%+;B*  
            for i=ia+1:ib HhZ>/5'(  
                curl_h=hz(i,j,k)-hz(i,j-1,k)-hy(i,j,k)+hy(i,j,k-1); =h>jo&=Wad  
                dx(i,j,k)=gj3(j)*gk3(k)*dx(i,j,k)+gj2(j)*gk2(k)*0.5*curl_h; ,%T sfB  
            end D[ v2#2  
            for i=ib+1:IE 'QdDXw5o  
                ixh=i-ib; ysH'X95  
                curl_h=hz(i,j,k)-hz(i,j-1,k)-hy(i,j,k)+hy(i,j,k-1); jmBsPSGIC  
                idxh(ixh,j,k)=idxh(ixh,j,k)+curl_h; :^En\YcU  
                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)); yog(  
            end LOEiV  
        end K'c[r0Ew  
    end ]<BT+6L  
   "Ng%"N
z  
    for i=ia+1:ib-1 xFScj0Y  
        for k=ka+1:kb 5F78)qu6N  
            dx(i,ja+1,k)=dx(i,ja+1,k)-0.5*hz_left(i,k); m4on<5s/  
            dx(i,jb,k)=dx(i,jb,k)+0.5*hz_right(i,k); M:*)l(  
        end iu iVr$E  
    end rqWD#FB=z  
    5v[2R.eT-  
    for i=ia+1:ib-1 w }=LC#le  
        for j=ja+1:jb ts/Ha*h  
            dx(i,j,ka+1)=dx(i,j,ka+1)+0.5*hy_down(i,j); mhgvN-? "h  
            dx(i,j,kb)=dx(i,j,kb)-0.5*hy_up(i,j); M,vCAZ  
        end c6Z"6-}$  
    end 7
PbwCRg  
    {O!B8a   
    %%% Dy *5KDu$'(e  
    for i=2:IE 3rj7]:Vr  
        for k=2:KE B'-n ^';  
            for j=2:ja u#+Is4Vh  
                curl_h=hx(i,j,k)-hx(i,j,k-1)-hz(i,j,k)+hz(i-1,j,k); =/+f3  
                idyl(i,j,k)=idyl(i,j,k)+curl_h; S[
e> 8  
                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)); v05$"Ig  
            end `D? &)Y  
            for j=ja+1:jb j" 5 +"j  
                curl_h=hx(i,j,k)-hx(i,j,k-1)-hz(i,j,k)+hz(i-1,j,k); no,b_0@N  
                dy(i,j,k)=gi3(i)*gk3(k)*dy(i,j,k)+gi2(i)*gk2(k)*0.5*curl_h; ?&JKq^9\I  
            end Tl L,dPM  
            for j=jb+1:JE EX/{W$ &K  
            jyh=j-jb; $EnBigb!  
            curl_h=hx(i,j,k)-hx(i,j,k-1)-hz(i,j,k)+hz(i-1,j,k); -`CE;  
            idyh(i,jyh,k)=idyh(i,jyh,k)+curl_h; BC)1FxsGf  
            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)); b?j\YX[e  
            end G.:QA}FE'  
        end Y~B-dx'V  
    end 2Ah B)8bG  
    Xd@ d$  
    for i=ia+1:ib HE:]zH  
        for j=ja+1:jb-1 7n[0)XR>  
            dy(i,j,ka+1)=dy(i,j,ka+1)-0.5*hx_down(i,j); MmQk@~  
            dy(i,j,kb)=dy(i,j,kb)+0.5*hx_up(i,j); J(5#fo{Q.g  
        end 6Zx)L|B  
    end HP,{/ $i:  
    ^4C djMF-E  
    for j=ja+1:jb-1 GGU>={D)  
        for k=ka+1:kb }E^k*S  
            dy(ia+1,j,k)=dy(ia+1,j,k)+0.5*hz_back(j,k); x3l~kZ(  
            dy(ib,j,k)=dy(ib,j,k)-0.5*hz_front(j,k); Hyb_>n  
        end wW, n~W  
    end Y^QG\6q  
   }BWT21'-Y  
    %%%  Dz C ~Doj  
    for i=2:IE 5i-VnG  
        for j=2:JE .1:B\R((  
            for k=2:ka .naSK`J,`  
                curl_h=hy(i,j,k)-hy(i-1,j,k)-hx(i,j,k)+hx(i,j-1,k); (z:qj/|  
                idzl(i,j,k)=idzl(i,j,k)+curl_h; 8'Iei78Ov  
                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)); Zg3 /,:1  
            end EvptGM  
            for k=ka+1:kb Z4dl'v)9  
                curl_h=hy(i,j,k)-hy(i-1,j,k)-hx(i,j,k)+hx(i,j-1,k); b/d1(B@  
                dz(i,j,k)=gi3(i)*gj3(j)*dz(i,j,k)+gi2(i)*gj2(j)*0.5*curl_h; |~B`[p]5H  
            end "..I$R  
            for k=kb+1:KE moCR64n  
                kzh=k-kb; I&f!>y?,Z  
                curl_h=hy(i,j,k)-hy(i-1,j,k)-hx(i,j,k)+hx(i,j-1,k); 'F^1)Ga$  
                idzh(i,j,kzh)=idzh(i,j,kzh)+curl_h; D$Ao-6QE W  
                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)); g*;zVi  
            end ub]s>aqy   
        end - WQ)rz  
    end ^|(VI0KO  
    ZKJhmk  
    for i=ia+1:ib 0|| 5r#  
        for k=ka+1:kb-1 OCu/w1bc  
            dz(i,ja+1,k)=dz(i,ja+1,k)+0.5*hx_left(i,k); +t8#rT ^B  
            dz(i,jb,k)=dz(i,jb,k)-0.5*hx_right(i,k); Y+DVwz
$  
        end @!*I mN
MI  
    end <q`|,mc  
    6J <.i  
    for j=ja+1:jb s;#,c(  
        for k=ka+1:kb-1 A"6

&  
            dz(ia+1,j,k)=dz(ia+1,j,k)-0.5*hy_back(j,k); xk7VuS*  
            dz(ib,j,k)=dz(ib,j,k)+0.5*hy_front(j,k); M9.FtQhK/  
        end t;}`~B  
    end 0py29>"t  
    ! u9LZ  
    %%%  Calculate E from D ?(Xy 2%v  
    for i=2:IE-1 $U"pdf  
        for j=2:JE-1 :Q}Zb,32  
            for k=2:KE-1 T%yGSk  
                ex(i,j,k)=gax(i,j,k)*(dx(i,j,k)-ix(i,j,k)); ~!P&LZ  
                ix(i,j,k)=ix(i,j,k)+gbx(i,j,k)*ex(i,j,k); 8(]q/g"O  
                -Np}<O`./  
                ey(i,j,k)=gay(i,j,k)*(dy(i,j,k)-iy(i,j,k)); C,B{7s0-  
                iy(i,j,k)=iy(i,j,k)+gby(i,j,k)*ey(i,j,k); EUbyQL  
                deOk>v&U  
                ez(i,j,k)=gaz(i,j,k)*(dz(i,j,k)-iz(i,j,k)); ^@)*voP#G  
                iz(i,j,k)=iz(i,j,k)+gbz(i,j,k)*ez(i,j,k); x$z>.4  
            end i)(-Ad_  
        end ^;[^L=}8$  
    end A&V'WahC@I  
    y;fnC5Q  
    %%% Calculate the Fourier transform of Ex = O|}R
 
    
oK3PA  
    %%% Calculate the Incident field a28`)17z  
    for j=1:END-1 d wku6lCk  
        h_inc(j)=h_inc(j)+0.5*(e_inc(j)-e_inc(j+1)); -0+h&CO  
    end I:M15  
    hX:yn:P~  
    %%% Hx ;P^}2i[q>[  
    for j=1:JE-1 |?v+8QL,;t  
        for k=1:KE-1 n2Y a'YF  
            for i=1:ia }"STc&1  
                curl_e=ey(i,j,k+1)-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k); P"|-)d  
                ihxl(i,j,k)=ihxl(i,j,k)+curl_e; A#gy[.Bb  
                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)); h>A~yDT[  
            end usw(]CnH  
            for i=ia+1:ib bZ}T;!U?I  
                curl_e=ey(i,j,k+1)-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k); \eXuNv_  
                hx(i,j,k)=fj3(j)*fk3(k)*hx(i,j,k)+fj2(j)*fk2(k)*0.5*curl_e; GQYB2{e>  
            end ,WE2MAjhT  
            for i=ib+1:IE cvC 7#i[G  
                ixh=i-ib; 2LS91  
                curl_e=ey(i,j,k+1)-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k); 'VV"$`Fu"  
                ihxh(ixh,j,k)=ihxh(ixh,j,k)+curl_e; v)~!HCG  
                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)); $49;\pBZl  
            end CYdYa|  
        end 6M[OEI5  
    end or(P?Ro  
    %<1_\N7  
    for i=ia+1:ib t\O#5mo  
        for k=ka+1:kb-1 4JL]?75  
            hx(i,ja,k)=hx(i,ja,k)+0.5*ez_left(i,k); UYGO|lkEU  
            hx(i,jb,k)=hx(i,jb,k)-0.5*ez_right(i,k); |`d-;pk!%  
        end eC_i]q&o|  
    end >6|Xvtf  
    HE-ErEtGB  
    for i=ia+1:ib MC1&X
'  
        for j=ja+1:jb-1 q:mqA$n  
            hx(i,j,ka)=hx(i,j,ka)-0.5*ey_down(i,j); {d%hkbN+{  
            hx(i,j,kb)=hx(i,j,kb)+0.5*ey_up(i,j); (d['f]S+&  
        end F#z1 sl'  
    end PJAM_K;  
        +*G<xW :M  
    %%% Hy dvLL~VP  
    for i=1:IE-1 _pkmHj(  
        for k=1:KE-1 ])
#?rRw  
            for j=1:ja QG5c>Q  
                curl_e=ez(i+1,j,k)-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k); =WK's8FB;8  
                ihyl(i,j,k)=ihyl(i,j,k)+curl_e; Y8/&1s_  
                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)); 0y%s\,PsT  
            end d~y]7h|  
            for j=ja+1:jb :H/Rhx=  
                curl_e=ez(i+1,j,k)-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k); <1<0odB  
                hy(i,j,k)=fi3(i)*fk3(k)*hy(i,j,k)+fi2(i)*fk2(k)*0.5*curl_e; d!y_N&z|(  
            end IO"q4(&;P4  
            for j=jb+1:JE OG^#e+  
                jyh=j-jb; wq]vcY9^  
                curl_e=ez(i+1,j,k)-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k); 8-9<r  
                ihyh(i,jyh,k)=ihyh(i,jyh,k)+curl_e; T/tCX[}  
                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)); 'JJ :  
            end juMHc$d17  
        end WL;2&S/{@  
    end < Q6  
    &H%z1Lp  
    for i=ia+1:ib-1 _xM3c&VeG  
        for j=ja+1:jb h5lngw  
            hy(i,j,ka)=hy(i,j,ka)+0.5*ex_down(i,j); V=E5pB`Pr  
            hy(i,j,kb)=hy(i,j,kb)-0.5*ex_up(i,j); o[g]Va*8  
        end ^ R3g7 DG  
    end TlC?
?#  
    w}?,N  
    for j=ja+1:jb |5*:ThC[  
        for k=ka+1:kb-1 &<# ,J4  
            hy(ia,j,k)=hy(ia,j,k)-0.5*ez_back(j,k); D 1(9/;9  
            hy(ib,j,k)=hy(ib,j,k)+0.5*ez_front(j,k); <-HWs@8#  
        end [6%y RQ_  
    end 58]t iP"  
    ySHio;g9  
   %%% Hz 1%7zCM0s  
    for i=1:IE-1 X\5EF7:S  
        for j=1:JE-1 :JK+V2B$H  
            for k=1:ka Q@rlqWgU ~  
                curl_e=ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k); 8< "lEL|  
                ihzl(i,j,k)=ihzl(i,j,k)+curl_e; X3l6b+p  
                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)); ,S1'SCwVdJ  
            end L"}2Y3
 
            for k=ka+1:kb /U="~{*-R  
                curl_e=ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k); FLQ^J3A,I  
                hz(i,j,k)=fi3(i)*fj3(j)*hz(i,j,k)+fi2(i)*fj2(j)*0.5*curl_e; {<Xl57w-Q  
            end %V92q0XW  
            for k=kb+1:KE pEY zB;  
                kzh=k-kb; y27MG  
                curl_e=ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k); _3{8Zg  
                ihzh(i,j,kzh)=ihzh(i,j,kzh)+curl_e; |9mGX9q  
                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)); +XAM2uN5_.  
            end '&Tz8.jp~  
        end v>5TTL~
?  
    end BLb'7`t  
    oP:/%  
    for i=ia+1:ib-1 =g)SZK  
        for k=ka+1:kb ^geY Ay  
            hz(i,ja,k)=hz(i,ja,k)-0.5*ex_left(i,k); 
MpJ]1  
            hz(i,jb,k)=hz(i,jb,k)+0.5*ex_right(i,k); 8< z  
        end /p)y!5e  
    end >~uKkQ_p  
    E#\'$@8j
 
    for j=ja+1:jb-1 W\O.[7JP  
        for k=ka+1:kb dl+c+
w"  
            hz(ia,j,k)=hz(ia,j,k)+0.5*ey_back(j,k); 21hTun"W  
            hz(ib,j,k)=hz(ib,j,k)-0.5*ey_front(j,k); j#9n.i %h  
        end P?zL`czWd  
    end X + B=?|M  
    JR|P
]}  
    %%%  Visualize fields Hm_&``='  
    %timestep=int2str(N); SA~oGgk=P  
    %tview(:,:)=ez(:,:,kc); &zHY0
fxX  
    %pcolor(tview); 3
#idXc  
    %shading interp; <H|]^An!H  
    %caxis([-1,1]); .19_EQ>+  
    %colorbar; 7 ;x to =  
    %axis image; D"P<;@ef  
    %title(['Ez(i,j,k=kc), time step = ',timestep]); UC!"1)~mt`  
    %xlabel('i'); o3.b='HAm  
    %ylabel('j'); DK<}q1
xi  
    %pause(0.05);   BWzo|isv  
        L]=LY  
if N==T-dT 2K2_
-  
    eio=e_inc(is); 2Iq*7n:v0  
    hio=e_inc(is); J:\O .F#Fi  
    for i=iao:ibo 1(/rg  
        for j=jao:jbo            q ;e/gP2  
            exi_oup(i,j)=ex(i,j,kbo); r:,"k:C  
            eyi_oup(i,j)=ey(i,j,kbo); Lp{/  
            `91?^T;\F  
            exi_odown(i,j)=ex(i,j,kao); ]"i^VVw  
            eyi_odown(i,j)=ey(i,j,kao); LF!S`|FF  
            S>R40T=e  
            hxi_oup(i,j)=0.5*(hx(i,j,kbo)+hx(i,j,kbo-1)); Zc=#Y  
            hyi_oup(i,j)=0.5*(hy(i,j,kbo)+hy(i,j,kbo-1)); z*WQ=l2  
            7+m.:~H3}  
            hxi_odown(i,j)=0.5*(hx(i,j,kao)+hx(i,j,kao-1)); $~T|v7Y%  
            hyi_odown(i,j)=0.5*(hy(i,j,kao)+hy(i,j,kao-1)); gO bP  
        end xsg55`  
    end WU6F-{M"?  
    for i=iao:ibo |q 0iX2W  
        for k=kao:kbo a'7RzN ,]  
            exi_oleft(i,k)=ex(i,jao,k); -9(nsaV  
            ezi_oleft(i,k)=ez(i,jao,k); JgB"N/Oz  
            ?IR+OCAA  
            exi_oright(i,k)=ex(i,jbo,k); NZuylQ)0  
            ezi_oright(i,k)=ez(i,jbo,k); / 16 r_l  
            VB~Do?]*k%  
            hxi_oleft(i,k)=0.5*(hx(i,jao,k)+hx(i,jao-1,k)); d4LH`@SUZ-  
            hzi_oleft(i,k)=0.5*(hz(i,jao,k)+hz(i,jao-1,k)); _j]vR  
            Rc~63![O.  
            hxi_oright(i,k)=0.5*(hx(i,jbo,k)+hx(i,jbo-1,k)); ;z9U
_  
            hzi_oright(i,k)=0.5*(hz(i,jbo,k)+hz(i,jbo-1,k)); UL3u2g;d  
        end `%=!_|  
    end R'aA\k-  
    for j=jao:jbo 5-"aK~@+  
        for k=kao:kbo u7 {R; QKw  
            eyi_ofront(j,k)=ey(ibo,j,k); y9K'(/  
            ezi_ofront(j,k)=ez(ibo,j,k); VpB+|%@p  
            D#11 N^-K  
            eyi_oback(j,k)=ey(iao,j,k); m!:sDQn{3  
            ezi_oback(j,k)=ez(iao,j,k); +P=I4-?eX  
            #|6M*;lN|  
            hyi_ofront(j,k)=0.5*(hy(ibo,j,k)+hy(ibo-1,j,k)); m[A$Sp_"-h  
            hzi_ofront(j,k)=0.5*(hz(ibo,j,k)+hz(ibo-1,j,k)); W<H<~wf#  
            @eRv`O"  
            hyi_oback(j,k)=0.5*(hy(iao,j,k)+hy(iao-1,j,k)); I_na^sh*  
            hzi_oback(j,k)=0.5*(hz(iao,j,k)+hz(iao-1,j,k)); Sp@{5  
        end %Pk@`t(3  
    end ~.Cu,>fV  
end Ki%RSW(_`  
if N==T
y&2O)z!B  
    ero=e_inc(is); 5I0j>{U&  
    hro=h_inc(is); <X]dR 6FT  
    for i=iao:ibo I0OfK3!^  
        for j=jao:jbo            N)X51;+  
            exr_oup(i,j)=ex(i,j,kbo); +  $/mh  
            eyr_oup(i,j)=ey(i,j,kbo); zmU>  
            Uy$?B"Z  
            exr_odown(i,j)=ex(i,j,kao); i$z*~SuM#  
            eyr_odown(i,j)=ey(i,j,kao); 0|~3\e/QV  
            34qfP{9!N  
            hxr_oup(i,j)=0.5*(hx(i,j,kbo)+hx(i,j,kbo-1)); ]dnB,  
            hyr_oup(i,j)=0.5*(hy(i,j,kbo)+hy(i,j,kbo-1)); I(+%`{Wv  
            Kk/cI6`W  
            hxr_odown(i,j)=0.5*(hx(i,j,kao)+hx(i,j,kao-1)); 't3nh  
            hyr_odown(i,j)=0.5*(hy(i,j,kao)+hy(i,j,kao-1)); Bmi9U  
        end - s0QEQ  
    end }z_7?dn/  
    for i=iao:ibo \PK}4<x}  
        for k=kao:kbo TL([hR _  
            exr_oleft(i,k)=ex(i,jao,k); n<MreKixE  
            ezr_oleft(i,k)=ez(i,jao,k); L_fiE3G|>  
            qZX\riR  
            exr_oright(i,k)=ex(i,jbo,k); =T|m#*{.L  
            ezr_oright(i,k)=ez(i,jbo,k); uG4Q\,R  
            rq^VOK
|L  
            hxr_oleft(i,k)=0.5*(hx(i,jao,k)+hx(i,jao-1,k)); k E-+#p  
            hzr_oleft(i,k)=0.5*(hz(i,jao,k)+hz(i,jao-1,k)); _Zbgmasb  
            LS4E.Xdn  
            hxr_oright(i,k)=0.5*(hx(i,jbo,k)+hx(i,jbo-1,k)); 5}`e"X  
            hzr_oright(i,k)=0.5*(hz(i,jbo,k)+hz(i,jbo-1,k)); 01o,9_|FL  
        end ;V4f6[<]'z  
    end ;h7
O_|<%  
    for j=jao:jbo DZV U!J  
        for k=kao:kbo {_X&{dZLX  
            eyr_ofront(j,k)=ey(ibo,j,k); +JY]J89  
            ezr_ofront(j,k)=ez(ibo,j,k); 'h;x>r  
            \yY2 mr  
            eyr_oback(j,k)=ey(iao,j,k); <O]B'Wc [  
            ezr_oback(j,k)=ez(iao,j,k); \Gy+y`   
            
I=}R Z9  
            hyr_ofront(j,k)=0.5*(hy(ibo,j,k)+hy(ibo-1,j,k)); T :X A  
            hzr_ofront(j,k)=0.5*(hz(ibo,j,k)+hz(ibo-1,j,k)); jR{Rd}QtQ  
            CM+/.y T  
            hyr_oback(j,k)=0.5*(hy(iao,j,k)+hy(iao-1,j,k)); RTeG\U  
            hzr_oback(j,k)=0.5*(hz(iao,j,k)+hz(iao-1,j,k)); O-]^_LV`  
        end `6y\.6j  
    end ;r3|EA35  
end .0]Odf:@  
    yUu+68Z6  
end ;Y`k-R:E6A  
eo=ero+sqrt(-1)*eio; Eo_;Nc  
eo=abs(eo); "$|ne[b2  
ho=(hro+sqrt(-1)*hio)*exp(-sqrt(-1)*0.5*omiga); d]Mjr2h  
for i=iao:ibo /7Ft1f  
    for j=jao:jbo i UW.$1l  
        ex_oup(i,j)=(exr_oup(i,j)+sqrt(-1)*exi_oup(i,j))/eo; &(rR
)cG  
        ey_oup(i,j)=(eyr_oup(i,j)+sqrt(-1)*eyi_oup(i,j))/eo; ]o*-|[^?  
        <84d Vg  
        ex_odown(i,j)=(exr_odown(i,j)+sqrt(-1)*exi_odown(i,j))/eo; K?acRi  
        ey_odown(i,j)=(eyr_odown(i,j)+sqrt(-1)*eyi_odown(i,j))/eo; HH'5kE0;d  
        *gN
)a%9  
        hx_oup(i,j)=(hxr_oup(i,j)+sqrt(-1)*hxi_oup(i,j))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; _u8d`7$*%  
        hy_oup(i,j)=(hyr_oup(i,j)+sqrt(-1)*hyi_oup(i,j))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; `Qo}4nuRs  
        ,98`tB0  
        hx_odown(i,j)=(hxr_odown(i,j)+sqrt(-1)*hxi_odown(i,j))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; ^I6GH?19>e  
        hy_odown(i,j)=(hyr_odown(i,j)+sqrt(-1)*hyi_odown(i,j))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; &BQ%df<y\  
    end A@JZK+WB}  
end >A1;!kGE#  
for i=iao:ibo R?+:Js
/  
    for k=kao:kbo Dhp|%_>  
        ex_oleft(i,k)=(exr_oleft(i,k)+sqrt(-1)*exi_oleft(i,k))/eo; L@8C t  
        ez_oleft(i,k)=(ezr_oleft(i,k)+sqrt(-1)*ezi_oleft(i,k))/eo; Q#*Pjl  
        FY<77i  
        ex_oright(i,k)=(exr_oright(i,k)+sqrt(-1)*exi_oright(i,k))/eo; :faB7wduW;  
        ez_oright(i,k)=(ezr_oright(i,k)+sqrt(-1)*ezi_oright(i,k))/eo; W&Kjh|[1QZ  
        !!UQ,yU  
        hx_oleft(i,k)=(hxr_oleft(i,k)+sqrt(-1)*hxi_oleft(i,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; IP l]$j>N  
        hz_oleft(i,k)=(hzr_oleft(i,k)+sqrt(-1)*hzi_oleft(i,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; [7gwJiK  
        ^0T[V-PgiD  
        hx_oright(i,k)=(hxr_oright(i,k)+sqrt(-1)*hxi_oright(i,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; gqan]b_  
        hz_oright(i,k)=(hzr_oright(i,k)+sqrt(-1)*hzi_oright(i,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; '@eH)wh@m)  
    end !<j
)D_  
end Zcxj.F(,  
for j=jao:jbo "SN+ ^`  
    for k=kao:kbo 2^=.jML[  
        ey_ofront(j,k)=(eyr_ofront(j,k)+sqrt(-1)*eyi_ofront(j,k))/eo; gD/% l[  
        ez_ofront(j,k)=(ezr_ofront(j,k)+sqrt(-1)*ezi_ofront(j,k))/eo; Fx'E"d  
        B;xGTl@8  
        ey_oback(j,k)=(eyr_oback(j,k)+sqrt(-1)*eyi_oback(j,k))/eo; a1# 'uS9W  
        ez_oback(j,k)=(ezr_oback(j,k)+sqrt(-1)*ezi_oback(j,k))/eo; ~QE-$;  
        )>rHM6-W  
        hy_ofront(j,k)=(hyr_ofront(j,k)+sqrt(-1)*hyi_ofront(j,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; S#^-VZ~U4x  
        hz_ofront(j,k)=(hzr_ofront(j,k)+sqrt(-1)*hzi_ofront(j,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; glP W9q,f  
        FK.Qj P:  
        hy_oback(j,k)=(hyr_oback(j,k)+sqrt(-1)*hyi_oback(j,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; Y!lc/[8  
        hz_oback(j,k)=(hzr_oback(j,k)+sqrt(-1)*hzi_oback(j,k))*exp(-sqrt(-1)*0.5*omiga)/eo/Z; V2_I=]p_  
    end %%f(R7n  
end jd?NN:7  
%%%%%%%%%  Far Field Calculate  %%%%%%% g'1ASMuR  
rk=2*pi/WL; -XARew  
fx=zeros(1,181); -K%~2M<  
fy=zeros(1,181); ~3s\Q%   
fz=zeros(1,181); rgXD>yu(  
fmx=zeros(1,181); z+%74O"c  
fmy=zeros(1,181);
n ^_B0Rkv  
fmz=zeros(1,181); (TFo]c  
rcs=zeros(1,181); {]dH+J7  
est=zeros(1,181); #Vh$u%q3  
efa=zeros(1,181); 9U~sRj=D  
fai_s=0*pi/180; IRl(H_.  
cps2=cos(fai_s); #dxS QmG  
sps2=sin(fai_s); cp<jwcc!  
ss=0; ?<-ins  
for sita_s=0:pi/180:pi dNgjM Q  
    ss=ss+1; 7.tEi}O&_g  
    sps1=sin(sita_s); `d.4L.],  
    cps1=cos(sita_s); !0 -[}vvU  
    rk=2*pi/WL; Yq0=4#_  
    %%% fx、fmx——left/right output boundry K44j-Ypb  
    for i=iao+1:ibo ad[oor/7|  
        for k=kao:kbo X3,+aL`  
            temp1=exp(sqrt(-1)*rk*((i-0.5)*sps1*cps2+jao*sps1*sps2+k*cps1)); ?tFsSU  
            fx(ss)=fx(ss)-(1-0.5*(k==kao)|(k==kbo))*ddx*ddx*hz_oleft(i,k)*temp1; SjB#"A5  
            fmx(ss)=fmx(ss)+(1-0.5*(k==kao)|(k==kbo))*ddx*ddx*ez_oleft(i,k)*temp1; 6pY<,7t0  
            2@#`x"0  
            temp2=exp(sqrt(-1)*rk*((i-0.5)*sps1*cps2+jbo*sps1*sps2+k*cps1)); y]TNjLpo$  
            fx(ss)=fx(ss)+(1-0.5*(k==kao)|(k==kbo))*ddx*ddx*hz_ori .. 7H5t!yk|9  
c-hhA%@Wq  

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

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

好的，谢谢！ v4Ga0]VN$8  
&<pKx!  
根据你的提醒调试，希望能找出错误。 ($nrqAv4