登 录
註 冊
论坛
微波仿真网
注册
登录论坛可查看更多信息
微波仿真论坛
>
时域有限差分法 FDTD
>
大家帮忙看看这个3D_FDTD介质球的RCS有什么问 ..
发帖
回复
1
2
3
4116
阅读
29
回复
大家帮忙看看这个3D_FDTD介质球的RCS有什么问题
离线
cc81372365
UID ：27752
注册：
2009-03-18
登录：
2023-11-09
发帖：
345
等级：
论坛版主
0楼
发表于: 2010-11-29 13:16:05
//****************************************************************//
-z>Z0viA
//
zR3 lX}g
//********* 3D_介质球散射（RCS）FDTD程序 **********//
Z>>gXh<e[
//
@1J51< x
//****************************************************************//
Fc A)RsMI*
//
p[af[!
//入射波为Y方向的平面波（theta=pi/2 , phi=pi/2), RCS观察角为（theta=pi/2 , phi = 0到2pi）
mK7^:(<.LO
// 数字0，1代表场值经过傅里叶变换后的复数的实部和虚部
!%Z)eO~Z
// 金属方柱置于中心位置（ic,jc）,激励源为高斯脉冲源
EyKkjEXx_
#include <iostream>
Et @=Ic^E
#include <iomanip>
hB2s$QS
#include <fstream>
@&i# S}%/
//#include <omp.h>
Xa@ _^oL
#include <time.h>
zezofW]a
#include <valarray>
2m)kyQ
#include <math.h>
" kE:T.,
#include <stdio.h>
m<,y-bQ*(
#include <stdlib.h>
lzr>WbM{{p
0,"n-5Im
#include "complex_operation.h"
V rnK)za*H
#define IE 60
)$9 C` d[
#define JE 60
4acP*LkkQ
#define KE 60
"FLD%3 l
#define ia 15
d @^o/w8
#define ja 15
)$lSG}WD
#define ka 15
I9Edw]
#define length 41
8a":[Q[
using namespace std;
[+:mt</HN
ofstream outFile;
2heWE
int main()
G]m[ S -
{
?C#F?N0
int i,j,k,ic,jc,kc,ib,jb,kb;
`/gEKrhL-
int npml,n_pml,T;
F=*BvI "+
double radius,v,lamda,Z,f,kk;
K 8yyxJ
double epsz_r,epsz,muz,muz_r,ddx,dt,pi,sigma;
pZ@W6}
double xn,xxn;
||*&g2Y
double t0,spread,pulse;
YP .%CD(K
double ez_inc[JE]={0},hx_inc[JE]={0};
EE{#S
double ez_low_m1=0,ez_low_m2=0,ez_high_m1=0,ez_high_m2=0;
:Z*02J wK
int ixh,jyh,kzh;
YN5OuKMUd'
int co=0;
:) Fp B"
double rcs[361]={0.0};
O_ s9
double ine[2] = {0.0};
=@ L5
double ine_fudu = 0.0;
Us@ {w`T
double curl_h,curl_e;
,,wyydG
double dist,xdist,ydist,zdist;
SS45<!i y
double dx[IE][JE][KE]={0},dy[IE][JE][KE]={0},dz[IE][JE][KE]={0};
lo>-}xd
double ex[IE][JE][KE]={0},ey[IE][JE][KE]={0},ez[IE][JE][KE]={0};
sa0^1$(<
double hx[IE][JE][KE]={0},hy[IE][JE][KE]={0},hz[IE][JE][KE]={0};
n{qw ]/
double ix[IE][JE][KE]={0},iy[IE][JE][KE]={0},iz[IE][JE][KE]={0};
=JySY@?9
double gax[IE][JE][KE]={0},gay[IE][JE][KE]={0},gaz[IE][JE][KE]={0};
"$,}|T?Y`
double gbx[IE][JE][KE]={0},gby[IE][JE][KE]={0},gbz[IE][JE][KE]={0};
!Wdt:MUI8
double idxl[ia][JE][KE]={0},idxh[ia][JE][KE]={0};
0N3tsIm>
double idyl[IE][ja][KE]={0},idyh[IE][ja][KE]={0};
@6I[{{>X
double idzl[IE][JE][ka]={0},idzh[IE][JE][ka]={0};
SCI-jf3 WN
double ihxl[ia][JE][KE]={0},ihxh[ia][JE][KE]={0};
)z*$`?)k
double ihyl[IE][ja][KE]={0},ihyh[IE][ja][KE]={0};
B 4*X0x
double ihzl[IE][JE][ka]={0},ihzh[IE][JE][ka]={0};
sGDrMAQt
double gi1[IE]={0},gi2[IE]={0},gi3[IE]={0};
J7 *G/F
double gj1[JE]={0},gj2[JE]={0},gj3[JE]={0};
a dz;N;rIY
double gk1[KE]={0},gk2[KE]={0},gk3[KE]={0};
]' "^M
double fi1[IE]={0},fi2[IE]={0},fi3[IE]={0};
"z(fBnv
double fj1[JE]={0},fj2[JE]={0},fj3[JE]={0};
2Xj-A\Oh~
double fk1[KE]={0},fk2[KE]={0},fk3[KE]={0};
!w;A=
double Ex_down[length][length][2]={0},Ey_down[length][length][2]={0};
-f |+
double Hx_down[length][length][2]={0},Hy_down[length][length][2]={0};
S#0|#Z5qD
double Ex_up[length][length][2]={0},Ey_up[length][length][2]={0};
i-"h"nF"
double Hx_up[length][length][2]={0},Hy_up[length][length][2]={0};
[V #&sAe
double Ex_left[length][length][2]={0},Ez_left[length][length][2]={0};
v%l|S{>(
double Hx_left[length][length][2]={0},Hz_left[length][length][2]={0};
v &CO#vK5.
double Ex_right[length][length][2]={0},Ez_right[length][length][2]={0};
t UAY]BJ*s
double Hx_right[length][length][2]={0},Hz_right[length][length][2]={0};
e G*s1uQl
double Ey_front[length][length][2]={0},Ez_front[length][length][2]={0};
P7BJ?x
double Hy_front[length][length][2]={0},Hz_front[length][length][2]={0};
ru6H nLhL
double Ey_back[length][length][2]={0},Ez_back[length][length][2]={0};
z*k 3q`=>
double Hy_back[length][length][2]={0},Hz_back[length][length][2]={0};
r /s&e e
ic = IE/2;
No92Y^~/
jc = JE/2;
%;B(_ht<-w
kc = KE/2;
d)17r\*>I
ib = IE - ia -1;
Ygk d~ g
jb = JE - ja -1;
}v( H E%~}
kb = KE - ka -1;
:;hBq4h
//外推边界位置
J G+g88
int dis = 5;
5@%=LPV
const int i_front = ib + dis;
g8iB;%6
const int i_back = ia - dis;
<g|nmu)o$
const int j_right = jb + dis;
\A ~I>x
const int j_left = ja - dis;
v_1JH<GJ-
const int k_up = kb + dis;
7PQ j7&m
const int k_down = ka - dis;
te)g',#lT
pi = 3.1415926;
ETH#IM8J
epsz = 8.8e-12; // 真空介电常数
9T0wdK]
epsz_r = 4.0; // 相对介电常数
[(e`b
muz = 4*pi*1e-7; // 真空磁导率
JY050FL
muz_r = 1.0; // 相对磁导率
8L9S^ '
sigma = 0; // 介质电导率
>`,#%MH#
v = 3.0e8; // 波速
dPdHY&#`
f = 9.375e9; // 频率 10GHz
u):Rw
lamda = v/f; // 波长 0.03m
m(7_ZiL=
kk = 2*pi/lamda; // 波数
r^o} Y
Z = sqrt(muz/epsz); // 真空波阻抗 377
Q&U= jX
ddx = lamda/42; // 网格尺寸
1r-,V X7
dt = ddx/(2.0*v); // 时间步长
CsST-qxg
/* Initialize the arrays */
$Bwvw)(%
for ( i=0; i<IE; i++)
/+JP ~ K
{
r<f-v_bxF
for ( j=0; j<JE; j++)
j|$ y)FBX
{
"6 Hj ji@A
for ( k=0; k<KE; k++)
?H7p6m u
{
"*>QxA%c4
ex
[j][k] = 0.0;
~BYEeUo;%v
ey
[j][k] = 0.0;
@l6 d J
ez
[j][k] = 0.0;
YWBP'Mo
dx
[j][k] = 0.0;
N(({ 2'Rr
dy
[j][k] = 0.0;
hwol7B>
dz
[j][k] = 0.0;
lj}1'K@M
hx
[j][k] = 0.0;
G3?a~n^b
hy
[j][k] = 0.0;
)mo|.L0
hz
[j][k] = 0.0;
cX =b q_
ix
[j][k] = 0.0;
?k7/`g U
iy
[j][k] = 0.0;
?$8 ,j +&I
iz
[j][k] = 0.0;
d7-F&!sQ
gax
[j][k] = 1.0;
}"; hz*a
gay
[j][k] = 1.0;
gu .))3D9
gaz
[j][k] = 1.0;
]I*#R9
gbx
[j][k] = 0.0;
@ :
gby
[j][k] = 0.0;
jQwg)E+o;
gbz
[j][k] = 0.0;
iNrmh iql
}
~zOU/8n ,F
}
4#Nd;gM2
}
Rx"VscB6z
for ( i=0; i<ia; i++)
UNH}*]u4`
{
Sm I8&c
for ( j=0; j<JE; j++)
pcxl2I
{
}y=n#%|i.
for ( k=0; k<KE; k++)
>5@ 0lYhH
{
qP .VK?jF|
idxl
[j][k] = 0.0;
E6&uZr
idxh
[j][k] = 0.0;
zZ"')+7q&%
ihxl
[j][k] = 0.0;
[?r\b
ihxh
[j][k] = 0.0;
&WWO13\qd
}
0#[Nfe*
}
`9"jHw`D
}
FVLA^$5c
for ( i=0; i<IE; i++)
Z(`K6`K M
{
FU v)<rK
for ( j=0; j<ja; j++)
1nM?> j%k
{
* jNu?$
for ( k=0; k<KE; k++)
"@ '9+$i6
{
Fb4S /_ V
idyl
[j][k] = 0.0;
GH)+y D[o
idyh
[j][k] = 0.0;
AxqTPx7`|
ihyl
[j][k] = 0.0;
gdn,nL`dP
ihyh
[j][k] = 0.0;
58gt*yVu
}
?B31 t9
}
|c+N)F B
}
nv|y@! (
for ( i=0; i<IE; i++)
HBkQ`T
{
JC?V].) y5
for ( j=0; j<JE; j++)
C g,w6<7
{
b:d.Lf{y7
for ( k=0; k<ka; k++)
O!cO/]<
{
&7y1KwfXn
idzl
[j][k] = 0.0;
WRyv >Y
idzh
[j][k] = 0.0;
T=KrT7
ihzl
[j][k] = 0.0;
C&K(({5O
ihzh
[j][k] = 0.0;
M(f*hOG{Y
}
K>p:?w
}
9!OCilG
}
Z" 8cGN'
//PML初始化
(YH{%8 Z0
for ( i=0; i<IE; i++)
Tf]VcE F
{
{}gk4 xr
gi1
= 0;
B!C32~[
fi1
= 0;
mp$II?hZ*
gi2
= 1;
v?' k)B
fi2
= 1;
!{=%l+^.
gi3
= 1;
sy#j+gZ
fi3
= 1;
?=o]Wx0(9
}
(HgdmN%
for ( j=0; j<JE; j++)
si4=C
{
*} 4;1OVT
gj1[j] = 0;
'$ nGtB5
fj1[j] = 0;
"Q!{8 9Y
gj2[j] = 1;
3mn-dKe((
fj2[j] = 1;
CM f~Yv
gj3[j] = 1;
j W]c9 u
fj3[j] = 1;
J x[e{o)o
}
^\t">NJ^
for ( k=0; k<KE; k++)
7.'j~hJL
{
C2%3+
gk1[k] = 0;
Z"&ODVP
fk1[k] = 0;
dufHd
gk2[k] = 1;
~bQ:gArk
fk2[k] = 1;
UVRV7^eTe
gk3[k] = 1;
8pKPbi;(2
fk3[k] = 1;
5N }|VGN
}
9jFDBy+
radius = 10 ;
L.&Vi"M <@
n_pml = 7; //PML层厚度
.:Xe* Q
npml = 7;
N9hWx()v
for (i = 0; i < n_pml; i++)
jtC ob'n8
{
<}>-ip?
xxn = (double)(npml-i)/(double)npml;
k5< 0M'
xn = 0.33*pow(xxn,3);
{.'g!{SHp
fi1
= xn;
[9H m][|Ph
fi1[IE-1-i] = xn;
"tbBbEj?d
gi2
= 1.0/(1.0 + xn);
B ~bU7.Cd
gi2[IE-1-i] = 1.0/(1.0 + xn);
UE](`|4H
gi3
= (1.0 - xn)/(1.0 + xn);
*VAi! 3Rx;
gi3[IE-1-i] = (1.0 - xn)/(1.0 + xn);
ag4^y&
xxn = (double)(npml-i-.5)/(double)npml;
d`*vJ#$> 2
xn = .33*pow(xxn,3);
3]82gZG G
gi1
= xn;
KUV{]?'
gi1[IE-2-i] = xn;
( I~Xw P&
fi2
= 1.0/(1.0 + xn);
J\},o|WI
fi2[IE-2-i] = 1.0/(1.0 + xn);
V^H47O;VC
fi3
= (1.0 - xn)/(1.0 + xn);
m/c~2?-;
fi3[IE-2-i] = (1.0 - xn)/(1.0 + xn);
}{PtQc6RL!
}
_E-GHj>k z
for (j = 0; j < n_pml; j++)
o[*ih\d
{
C /#/F#C
xxn = (double)(npml-j)/(double)npml;
Hd}t=6
xn = .33*pow(xxn,3);
: :t !W7W
fj1[j] = xn;
+n]Knfi
fj1[JE-1-j] = xn;
@ y %qQe/g
gj2[j] = 1.0/(1.0 + xn);
AYtcN4\/
gj2[JE-1-j] = 1.0/(1.0 + xn);
{_D'\i(Y_
gj3[j] = (1.0 - xn)/(1.0 + xn);
[_GR'x '0x
gj3[JE-1-j] = (1.0 - xn)/(1.0 + xn);
hIHO a
xxn = (double)(npml-j-.5)/(double)npml;
{i?K~| h
xn = .33*pow(xxn,3);
hRiGW_t
gj1[j] = xn;
Z.19v>-c
gj1[JE-2-j] = xn;
g4932_t C
fj2[j] = 1.0/(1.0 + xn);
4ov~y1Da)
fj2[JE-2-j] = 1.0/(1.0 + xn);
:~(^b;yhZ
fj3[j] = (1.0 - xn)/(1.0 + xn);
N1+]3kt ~
fj3[JE-2-j] = (1.0 - xn)/(1.0 + xn);
V3nv5/6
}
z>q_]U0
r+obm)Qtp
for (k = 0; k < n_pml; k++)
jF0BWPL
{
Q 7?4GxMj
xxn = (double)(npml-k)/(double)npml;
^E8Hv
xn = .33*pow(xxn,3);
+DDvM;31w
fk1[k] = xn;
UzRF'<TWf
fk1[KE-1-k] = xn;
!T1i_
gk2[k] = 1.0/(1.0 + xn);
,*#M%Pv1t
gk2[KE-1-k] = 1.0/(1.0 + xn);
dv>zK#!
gk3[k] = (1.0 - xn)/(1.0 + xn);
U; Wm x
gk3[KE-1-k] = (1.0 - xn)/(1.0 + xn);
x@x@0k`A2
xxn = (double)(npml-k-.5)/(double)npml;
0 N^V&k
xn = .33*pow(xxn,3);
lK SI5d
gk1[k] = xn;
4iPg_+
gk1[KE-2-k] = xn;
.5zJ bZ9
fk2[k] = 1.0/(1.0 + xn);
~~ w4854
fk2[KE-2-k] = 1.0/(1.0 + xn);
Kqjeqr@)
fk3[k] = (1.0 - xn)/(1.0 + xn);
mQ)l`w Gh
fk3[KE-2-k] = (1.0 - xn)/(1.0 + xn);
e1 a %Rj~
}
?q6eV~P
;#Nci%<J\
/* Calculate gax ,gbx */
uSbg*OA
for ( i = ia; i<=ib; i++)
kF7Al]IgT
{
Np)!23 "
for ( j = ja; j<=jb; j++)
L 32ki}2
{
>l[N]CQ
for ( k =ka; k<=kb; k++)
Gj0NN:
{
YRX e j
xdist = (ic-i-.5);
2Ki/K(
ydist = (jc-j);
'f!Jh<i
zdist = (kc-k);
q3:' 69
dist = sqrt (pow(xdist,2) + pow(ydist,2) + pow(zdist,2));
mcracj[ B
if ( dist <= radius )
r-uI FhV^
{
m]vr|:{6/
gax
[j][k] = 1./(epsz_r + (sigma*dt/epsz));
ZaV@}=Rd8
gbx
[j][k] = sigma*dt/epsz;
(/Bkwb JyE
}
3 s Nq3I
}
S]{Z_|h*j
}
YDL)F<Y
}
IVKE dwA
/* Calculate gay ,gby*/
]YKWa"
for ( i = ia; i<=ib; i++)
8t"DQ Y-R
{
y%AJ>@/;
for ( j = ja; j<=jb; j++)
7u(i4O& k
{
U3QnWPt}>
for ( k =ka; k<=kb; k++)
p*,mwKN:
{
F3'G9Xf8Q=
xdist = (ic-i);
%;gW l1&5
ydist = (jc-j-.5);
Fx\Re]~n
zdist = (kc-k);
YEj U3^@
dist = sqrt (pow(xdist,2) + pow(ydist,2) + pow(zdist,2));
n40 Z
if ( dist <= radius )
> skS`/6
{
DKH9 O
gay
[j][k] = 1./(epsz_r + (sigma*dt/epsz));
)<&QcO_
gby
[j][k] = sigma*dt/epsz;
h)8_sC
}
K a jyQ"j
}
"+60B0>sc
}
HBu>BSv:
}
C, rZ }-
/* Calculate gaz ,gbz */
bvKi0-
for ( i = ia; i<=ib; i++)
'r0gqtB
{
'%4,!
for ( j = ja; j<=jb; j++)
`V bG%y&I
{
z@ j&v W
for ( k =ka; k<=kb; k++)
_!D$Aj
{
.j.=|5nVo4
xdist = (ic-i);
N"{o3QmA
ydist = (jc-j);
$P~ a
zdist = (kc-k-.5);
VN;M ;fMs
dist = sqrt (pow(xdist,2) + pow(ydist,2) + pow(zdist,2));
AJ}Q,E
if ( dist <= radius )
F'!}$oT"
{
Wov_jVdN\
gaz
[j][k] = 1./(epsz_r +(sigma*dt/epsz));
.?TVBbc%5
gbz
[j][k] = sigma*dt/epsz;
VP&lWPA}\$
}
ShP V!$0
}
`.XU|J*z,
}
Ab)7h CUW
}
Z5K,y19/~
//源参数
cPSpPx
t0 = 40.0;
M`F L&Ac
spread = 10.0;
bx > D
"M]]H^ r5
clock_t start, finish;
d+;~x*
start = clock();
d 0d2QRX
/* outFile.open("F:\\ez[40][25][25].txt");*/
j3U 8@tuG
for ( T = 1; T < 1001; T++ )
D1a4+AyI
{
,\b5M`<c
/* Calculate the incident buffer */
aDESO 5
for ( j=1; j<JE; j++)
i\b2P2 `B
{
sj~'.Zs%
ez_inc[j] = ez_inc[j] + .5*(hx_inc[j-1] - hx_inc[j]);
:gn!3P}p?
}
%,BJkNV
//Fourier Transform of the incident field
xOH@V4 z:
ine[0] = ine[0] + dt*cos(2*pi*f*dt*T)*ez_inc[ja-1];
<T>f@Dn,
ine[1] = ine[1] - dt*sin(2*pi*f*dt*T)*ez_inc[ja-1];
.u&g2Y
/* Source */
5Ml}m
//pulse = exp (-.5*(pow(( t0-T )/spread,2)));
jw)t"S/E
pulse = exp (2*pi*f*T*dt);
$ hapSrS
ez_inc[3] = pulse;
Bd i~ B")
/* Boundary conditions for the incident buffer */
#kA+Yqy \)
ez_inc[0] = ez_low_m2;
c2QC`h(Wb
ez_low_m2 = ez_low_m1;
YDJc@*D
ez_low_m1 = ez_inc[1];
`7`iCYiTy
ez_inc[JE-1] = ez_high_m2;
n@f@-d$m\<
ez_high_m2 = ez_high_m1;
~}fpe>M:
ez_high_m1 = ez_inc[JE-2];
uU0'y4=
/* Calculate the Dx field */
TaH9Nu
for ( i=1; i<ia; i++)
W=/B[@3'
{
r=;k[*;{
for ( j=1; j<JE; j++)
:1Ay_ b_J
{
9,IGZ55C
for (k=1; k<KE; k++)
Bb6_[ 'y
{
3$Tp I5A
curl_h = ( hz
[j][k] - hz
[j-1][k] - hy
[j][k] + hy
[j][k-1] );
;S/fe(C
idxl
[j][k] = idxl
[j][k] + curl_h;
0FB ifK
dx
[j][k] = gj3[j]*gk3[k]*dx
[j][k] + gj2[j]*gk2[k]*.5*(curl_h + gi1
*idxl
[j][k]);
d>f5T l\E
}
W=drp>Uj
}
uWSG+
}
Z3Y%VHB_F(
for ( i=ia; i<=ib; i++)
a,.9 eHf
{
Q6 E80>
for ( j=1;j<JE; j++)
B s*s8}6
{
vnE , }(M
for ( k=1; k<KE; k++)
O \?ei+(H7
{
:HC{6W`$
curl_h = ( hz
[j][k] - hz
[j-1][k] - hy
[j][k] + hy
[j][k-1] );
sE% n=Ww
dx
[j][k] = gj3[j]*gk3[k]*dx
[j][k] + gj2[j]*gk2[k]*.5*curl_h;
i*3'O :Gq
}
e5dw q
}
!#QD;,SE+
}
b1 +N m
for ( i=ib+1; i<IE; i++)
/@K?W=w4
{
JkTL+obu
ixh = i-ib-1;
@U,cj>K
for ( j=1; j<JE; j++)
TSeAC[%pL
{
gyIPG2d
for ( k=1; k<KE; k++)
!Qy3fs
{
'` n\YO.N
curl_h = ( hz
[j][k] - hz
[j-1][k] - hy
[j][k] + hy
[j][k-1] );
iPao54Z
idxh[ixh][j][k] = idxh[ixh][j][k] + curl_h;
h4U .wk
dx
[j][k] = gj3[j]*gk3[k]*dx
[j][k] + gj2[j]*gk2[k]*.5*(curl_h + gi1
*idxh[ixh][j][k]);
pLRHwL.
}
(Xx @_
}
Z=ho7i
}
nZ]d[
/* Calculate the Dy field */
Rh39x-`Z
for ( i=1; i<IE; i++)
oPi)#|jcb
{
a,X3=+_ K
for ( j=1; j<ja; j++)
<5 okwcJ^
{
1mf_1 spB
for (k=1; k<KE; k++)
"qF8'58
{
&A>J> b
curl_h = ( hx
[j][k] - hx
[j][k-1] - hz
[j][k] + hz[i-1][j][k] );
-1[ri8t;nV
idyl
[j][k] = idyl
[j][k] + curl_h;
vjW S35i
dy
[j][k] = gi3
*gk3[k]*dy
[j][k]+gi2
*gk2[k]*.5*(curl_h+ gj1[j]*idyl
[j][k]);
4&tY5m >
}
I>:'5V
}
|e!Sm{#!
}
2HUoT\M
for ( i=1; i<IE; i++)
Xa36O5$4]9
{
) s M}BY
for ( j=ja;j<=jb; j++)
,%u\2M
{
nG2RBeJV
for ( k=1; k<KE; k++)
4JucNGv
{
u >x2
curl_h = ( hx
[j][k] - hx
[j][k-1] - hz
[j][k] + hz[i-1][j][k] );
H4UnF5G
dy
[j][k] = gi3
*gk3[k]*dy
[j][k]+gi2
*gk2[k]*.5*curl_h;
v+Hu =RZE
}
]>!]X*\9
}
7kO 1d{u6b
}
LGK}oL'
for ( i=1; i<IE; i++)
w1.~N`g$
{
-k'=s{iy
for ( j=jb+1; j<JE; j++)
P;h/)-q8
{
7l?-2I'c
jyh = j-jb-1;
;kdJxxUox
for ( k=1; k<KE; k++)
'8Yx
{
8<xJmcTEwO
curl_h = ( hx
[j][k] - hx
[j][k-1] - hz
[j][k] + hz[i-1][j][k] );
Upf1*$p
idyh
[jyh][k] = idyh
[jyh][k] + curl_h;
9K)2OX;$w
dy
[j][k] = gi3
*gk3[k]*dy
[j][k]+gi2
*gk2[k]*.5*(curl_h+ gj1[j]*idyh
[jyh][k]);
]RV6( |U4_
}
mIOx)`$
}
p[c C%3
}
I;=}@]9
/* Incident Dy */
0&<{o!>k
for ( i=ia; i<=ib; i++)
.S[5CO^
{
qk+:p]2
for ( j=ja; j<=jb-1; j++)
u4%-e )$X
{
P )`-cfg
dy
[j][ka] = dy
[j][ka] - .5*hx_inc[j];
qr :[y
dy
[j][kb+1] = dy
[j][kb+1] + .5*hx_inc[j];
lgU7jn
}
V XC _Y
}
*<J**FhcMu
/* Calculate the Dz field */
y2)~ljR
for ( i=1; i<IE; i++)
} pA0mW9
{
$L(,q!DvH
for ( j=1; j<JE; j++)
o|*,<5t
{
bBGLf)fsTG
for (k=1; k<ka; k++)
[mX\Q`)QP
{
8L0#<"'0
curl_h = ( hy
[j][k] - hy[i-1][j][k] - hx
[j][k] + hx
[j-1][k] );
<[(xGrEZV
idzl
[j][k] = idzl
[j][k] + curl_h;
)U5AnL
dz
[j][k] = gi3
*gj3[j]*dz
[j][k]+gi2
*gj2[j]*.5*(curl_h+ gk1[k]*idzl
[j][k]);
rN OwB2e
}
iYFM@ta
}
V F;%Z
}
(ue;O~
for ( i=1; i<IE; i++)
\SyfEcSf2v
{
jQc$>M<"o
for ( j=1;j<JE; j++)
%ePInpb
{
ug OcK Gf
for ( k=ka; k<=kb; k++)
,w c|YI)E
{
w q JH
curl_h = ( hy
[j][k] - hy[i-1][j][k] - hx
[j][k] + hx
[j-1][k] );
VsFRG;:\U
dz
[j][k] = gi3
*gj3[j]*dz
[j][k]+gi2
*gj2[j]*.5*curl_h;
d^4!=^HN
}
^,,|ED\M{m
}
" @.hz@>
}
l]GLkE
for ( i=1; i<IE; i++)
; R} :2
{
&r jMGk"&
for ( j=1; j<JE; j++)
B \BP:;"
{
1_ %3cN.
for ( k=kb+1; k<KE; k++)
RxcX \:
{
s(-$|f+s
kzh = k-kb-1;
NU81 V0:jG
curl_h = ( hy
[j][k] - hy[i-1][j][k] - hx
[j][k] + hx
[j-1][k] );
L ej3? k
idzh
[j][kzh] = idzh
[j][kzh] + curl_h;
ho 4~-xmN
dz
[j][k] = gi3
*gj3[j]*dz
[j][k]+gi2
*gj2[j]*.5*(curl_h+ gk1[k]*idzh
[j][kzh]);
Y'58.8hl
}
fi`*r\
}
g+U6E6}1
}
~r+;i ,,X
/* Incident Dz */
(|<}q-wO
for ( i=ia; i<=ib; i++)
T2GJ oJ!
{
ge*f<#|0U-
for ( k=ka; k<=kb; k++)
GkVV%0;&J1
{
7N 7W0Ky
dz
[ja][k] = dz
[ja][k] + .5*hx_inc[ja-1];
L -<!,CASW
dz
[jb][k] = dz
[jb][k] - .5*hx_inc[jb];
Lk4g js, V
}
pFhznH{0
}
;=aj)lemCr
/* Calculate the E from D field */
1ih|b8)Dn
for ( i=1; i<IE-1; i++)
2e#hJ-/`-
{
Z+JPxe#7
for ( j=1; j<JE-1; j++)
bvi Y.G3
{
eWE7 >kwh
for ( k=1; k<KE-1; k++)
|} =xA%)
{
*$Bx#0J8
ex
[j][k] = gax
[j][k]*(dx
[j][k]-ix
[j][k]);
IOomBy:
ix
[j][k] = ix
[j][k]+gbx
[j][k]*ex
[j][k];
y @e/G 3
ey
[j][k] = gay
[j][k]*(dy
[j][k]-iy
[j][k]);
ZE#A?5lb
iy
[j][k] = iy
[j][k]+gby
[j][k]*ey
[j][k];
rdSkGb
ez
[j][k] = gaz
[j][k]*(dz
[j][k]-iz
[j][k]);
b@ "#A8M
iz
[j][k] = iz
[j][k]+gbz
[j][k]*ez
[j][k];
>E6w,Ab
}
Dy mf
}
p{NVJ^! +
}
WQK ~;GV-
/* Remember: part of the PML is E = 0 at the edges */
m >DBO|`
for (j = 0;j < JE-1;j++)
V& _
{
lh{U@,/
for (k = 0;k < KE-1;k++)
)$XcO]
{
AKu]c-
ez[0][j][k] = 0.0;
\9jEpE^Ju(
ez[IE-1][j][k] = 0.0;
l*^c?lp)
}
H+6+ I53
}
3P N<J
for (i = 0;i < IE-1;i++)
w`x4i fZ0q
{
/{|JQ'gqX
for (k = 0;k < KE;k++)
67?O} ~jbG
{
tP^2NTs%]
ez
[0][k] = 0.0;
?t46TV'G
ez
[JE-1][k] = 0.0;
D.su^m_1
}
0M-=3 T
}
rUmaKh?v|X
for (i = 0;i < IE-1;i++)
16Cd0[h?
{
K +p7yZJ
for (j = 0;j < JE-1;j++)
@vs+)aRa
{
\ 2#7B8
ez
[j][0] = 0.0;
dv1Y2 [
ez
[j][KE-1] = 0.0;
v4OroG =^
}
k> SPtiAs
}
}fU"s"
rv ouE:
/*outFile << setw(22) << setprecision(8) << T << setw(22) << setprecision(8) << ez[40][25][25]<< endl;*/
b9X "p*'p
8`Ih> D c
/* Calculate the incident field */
&?"E"GH
for ( j=0; j<JE-1; j++)
EIug)S~
{
epXvk &
hx_inc[j] = hx_inc[j] + .5*(ez_inc[j]-ez_inc[j+1]);
clyZD`*
}
{el[W,CT#
/* Calculate the Hx field */
zsOOx% +
for ( i=0; i<ia; i++)
Z4KY VHD,
{
K+|G9
for ( j=0; j<JE-1; j++)
w kc)2z
{
.Yg7V'R1
for (k=0; k<KE-1; k++)
%m )vQ\Vtx
{
(.^ KuXd
curl_e = ( ey
[j][k+1] - ey
[j][k] - ez
[j+1][k] + ez
[j][k] );
Qa>t$`o`
ihxl
[j][k] = ihxl
[j][k] + curl_e;
)jm!bR`
hx
[j][k] = fj3[j]*fk3[k]*hx
[j][k]+fj2[j]*fk2[k]*.5*(curl_e+ fi1
*ihxl
[j][k]);
@kBy|5
}
2D;2QdO
}
'A2^K5`3
}
bS7%%8C
for ( i=ia; i<=ib; i++)
2K wr=t
{
1)Z4 (_
for ( j=0;j<JE-1; j++)
E?v:7p<
{
B :oE&Ahh{
for ( k=0; k<KE-1; k++)
^Q,-4\ec
{
_D.4=2@|l8
curl_e = ( ey
[j][k+1] - ey
[j][k] - ez
[j+1][k] + ez
[j][k] );
C)hS^ D:
hx
[j][k] = fj3[j]*fk3[k]*hx
[j][k]+fj2[j]*fk2[k]*.5*curl_e;
q |h#J}\
}
ON [F
}
t[}&*2"$/
}
r"YOA@
for ( i=ib+1; i<IE; i++)
1# }}:
{
D6N 32q@
ixh = i-ib-1;
[+; >u|
for ( j=0; j<JE-1; j++)
JP{Y Q:N F
{
@KJ~M3d0l
for ( k=0; k<KE-1; k++)
S3 Dmc\f
{
o&-D[|E|
curl_e = ( ey
[j][k+1] - ey
[j][k] - ez
[j+1][k] + ez
[j][k] );
0Sx$6:-~
ihxh[ixh][j][k] = ihxh[ixh][j][k] + curl_e;
p d^"MG
hx
[j][k] = fj3[j]*fk3[k]*hx
[j][k]+fj2[j]*fk2[k]*.5*(curl_e+ fi1
*ihxh[ixh][j][k]);
%/. yGAPkx
}
\l-JU
}
zPm|$d
}
@}\wec_
/* Incident Hx */
yVfF * nG
for ( i=ia; i<=ib; i++)
]Ub"NLYV
{
f 0M5^
for ( k=ka; k<=kb; k++)
&+|bAn9AJ
{
YUjKOPN
hx
[ja-1][k] = hx
[ja-1][k] + .5*ez_inc[ja];
5LaF'>1yY
hx
[jb][k] = hx
[jb][k] - .5*ez_inc[jb];
Tji* \<?
}
Y+)qb);
}
++\s0A(e
Ky[s& >02
/* Calculate the Hy field */
&a8%j+j
for ( i=0; i<IE-1; i++)
_D"V^4^yqu
{
$&Lw 2 c0
for ( j=0; j<ja; j++)
;suY
{
)kIZm Q|f1
for (k=0; k<KE-1; k++)
OjWg>v\ v
{
gk hm Qd
curl_e = ( ez[i+1][j][k] - ez
[j][k] - ex
[j][k+1] + ex
[j][k] );
.VXadgM
ihyl
[j][k] = ihyl
[j][k] + curl_e;
4LXC;gZ
hy
[j][k] = fi3
*fk3[k]*hy
[j][k]+fi2
*fk2[k]*.5*(curl_e+ fj1[j]*ihyl
[j][k]);
s+0n0C
}
,,8'29yEq
}
aE(DNeG-H
}
R+r;V ]-/
for ( i=0; i<IE-1; i++)
U2G[uDa;
{
2=IZD `{!
for ( j=ja;j<=jb; j++)
Yg V" *~
{
{ iXQUj
for ( k=0; k<KE-1; k++)
}y>/#]X
{
=C#22xqQ.
curl_e = ( ez[i+1][j][k] - ez
[j][k] - ex
[j][k+1] + ex
[j][k] );
?+_Y!*J2b
hy
[j][k] = fi3
*fk3[k]*hy
[j][k]+fi2
*fk2[k]*.5*curl_e;
XcFu:B
}
GahIR9_2
}
>zX `qv&>
}
k 5g vo
for ( i=0; i<IE-1; i++)
f D< 0V
{
ROi_k4Fj
for ( j=jb+1; j<JE; j++)
|\FJ
{
Qa#Em1co
jyh = j-jb-1;
#sm@|' Q%
for ( k=0; k<KE-1; k++)
q~. .Z Y`7
{
x~?|bnM#3
curl_e = ( ez[i+1][j][k] - ez
[j][k] - ex
[j][k+1] + ex
[j][k] );
blx"WVqo
ihyh
[jyh][k] = ihyh
[jyh][k] + curl_e;
IE:;`e:\D
hy
[j][k] = fi3
*fk3[k]*hy
[j][k]+fi2
*fk2[k]*.5*(curl_e+ fj1[j]*ihyh
[jyh][k]);
ak _n
}
:d1Kq _\K
}
lk4U/:
}
QC Jf
/* Incident Hy */
^V0I!&7lx
for ( j=ja; j<=jb; j++)
TUN6`/"
{
tS3!cO\
for ( k=ka; k<=kb; k++)
@HQqHO&N
{
!ZS5}/ZU
hy[ia-1][j][k] = hy[ia-1][j][k] - .5*ez_inc[j];
wd*V,ZN7
hy[ib][j][k] = hy[ib][j][k] + .5*ez_inc[j];
tq=7HM
}
L_5o7~` 0
}
&8HJ4Vj2
/* Calculate the Hz field */
|U`A So
for ( i=0; i<IE-1; i++)
o*5b]XWw
{
bQ`2ll*(
for ( j=0; j<JE-1; j++)
@{@DGc
{
3[UaK`/1C
for (k=0; k<ka; k++)
4Q (w D
{
o0kK f+[
curl_e = ( ex
[j+1][k] - ex
[j][k] - ey[i+1][j][k] + ey
[j][k] );
Lvb'qZ6n
ihzl
[j][k] = ihzl
[j][k] + curl_e;
a|-B# S
hz
[j][k] = fi3
*fj3[j]*hz
[j][k]+fi2
*fj2[j]*.5*(curl_e+ fk1[k]*ihzl
[j][k]);
8a;;MJ)
}
Hy1pIU sx
}
7V%b!R}
}
U!b~vrr^
for ( i=0; i<IE-1; i++)
0ZJrK\K;
{
|Hv8GT
for ( j=0;j<JE-1; j++)
,a3M*}Y ~3
{
R^6]v`j;
for ( k=ka; k<=kb; k++)
Wgm{ ]9Q
{
uoq|l
curl_e = ( ex
[j+1][k] - ex
[j][k] - ey[i+1][j][k] + ey
[j][k] );
jwp?eL!7
hz
[j][k] = fi3
*fj3[j]*hz
[j][k]+fi2
*fj2[j]*.5*curl_e;
iw^"?:'%
}
Y;=GM:*H
}
nNhb,J
}
!-Uq# Ea0/
for ( i=0; i<IE-1; i++)
vhrURY.
{
,8.zbr
for ( j=0; j<JE-1; j++)
uQl=?0 85
{
ts &sr
for ( k=kb+1; k<KE; k++)
| MXRNA~
{
T5eJIc3a"
kzh = k-kb-1;
Wb#ON|.2
curl_e = ( ex
[j+1][k] - ex
[j][k] - ey[i+1][j][k] + ey
[j][k] );
&UQP9wS4v
ihzh
[j][kzh] = ihzh
[j][kzh] + curl_e;
? @- t.N
hz
[j][k] = fi3
*fj3[j]*hz
[j][k]+fi2
*fj2[j]*.5*(curl_e+ fk1[k]*ihzh
[j][kzh]);
