登 录
註 冊
论坛
微波仿真网
注册
登录论坛可查看更多信息
微波仿真论坛
>
时域有限差分法 FDTD
>
大家帮忙看看这个3D_FDTD介质球的RCS有什么问 ..
发帖
回复
1
2
3
4115
阅读
29
回复
大家帮忙看看这个3D_FDTD介质球的RCS有什么问题
离线
cc81372365
UID ：27752
注册：
2009-03-18
登录：
2023-11-09
发帖：
345
等级：
论坛版主
0楼
发表于: 2010-11-29 13:16:05
//****************************************************************//
`FF8ie 8L
//
8G^<[`.@j
//********* 3D_介质球散射（RCS）FDTD程序 **********//
a$11PBi[9
//
Gs0H@
//****************************************************************//
uXZg1 F)
//
h) (* q+a
//入射波为Y方向的平面波（theta=pi/2 , phi=pi/2), RCS观察角为（theta=pi/2 , phi = 0到2pi）
7 h=QW 5
// 数字0，1代表场值经过傅里叶变换后的复数的实部和虚部
%Q"(/jm?
// 金属方柱置于中心位置（ic,jc）,激励源为高斯脉冲源
(GK pA}~R
#include <iostream>
cD}Sf>
#include <iomanip>
0JyVNuHn
#include <fstream>
o @Z#
//#include <omp.h>
XFw Lz
#include <time.h>
{=:#S+^ER
#include <valarray>
WY
#include <math.h>
K 7TzF&
#include <stdio.h>
f>9s!Hpu_
#include <stdlib.h>
<O~WB
sp9W?IJ 6c
#include "complex_operation.h"
vRQ Os0F;
#define IE 60
VRng=,
#define JE 60
a> mMvc"
#define KE 60
VU.@R,
#define ia 15
7N8H)X
#define ja 15
^V*-1r1
#define ka 15
Vzlh+R>c
#define length 41
c @(&[/q!
using namespace std;
A&5:ATQ/|
ofstream outFile;
fRZ KEIyk
int main()
kQ>^->w
{
JgRYljQi2
int i,j,k,ic,jc,kc,ib,jb,kb;
%8.J =B
int npml,n_pml,T;
77 r(*.O|
double radius,v,lamda,Z,f,kk;
pf" <!O[
double epsz_r,epsz,muz,muz_r,ddx,dt,pi,sigma;
\f WW'
double xn,xxn;
u eb-2[=
double t0,spread,pulse;
{d 3<W N
double ez_inc[JE]={0},hx_inc[JE]={0};
.10y0F L4
double ez_low_m1=0,ez_low_m2=0,ez_high_m1=0,ez_high_m2=0;
0$ .m _0H
int ixh,jyh,kzh;
L5fuM]G`
int co=0;
63WS7s"
double rcs[361]={0.0};
d3Di/Iej
double ine[2] = {0.0};
i&Fiq&V)[
double ine_fudu = 0.0;
|AhF7Mj *
double curl_h,curl_e;
< B g8,;
double dist,xdist,ydist,zdist;
_Q(g(p&
double dx[IE][JE][KE]={0},dy[IE][JE][KE]={0},dz[IE][JE][KE]={0};
V\5 L?}
double ex[IE][JE][KE]={0},ey[IE][JE][KE]={0},ez[IE][JE][KE]={0};
`RRE(SiKU
double hx[IE][JE][KE]={0},hy[IE][JE][KE]={0},hz[IE][JE][KE]={0};
R=j% S!
double ix[IE][JE][KE]={0},iy[IE][JE][KE]={0},iz[IE][JE][KE]={0};
M" lg%j
double gax[IE][JE][KE]={0},gay[IE][JE][KE]={0},gaz[IE][JE][KE]={0};
Z=S>0|`R
double gbx[IE][JE][KE]={0},gby[IE][JE][KE]={0},gbz[IE][JE][KE]={0};
"hz\Z0zg2
double idxl[ia][JE][KE]={0},idxh[ia][JE][KE]={0};
_D7 ]-3uC!
double idyl[IE][ja][KE]={0},idyh[IE][ja][KE]={0};
e( X|3h|
double idzl[IE][JE][ka]={0},idzh[IE][JE][ka]={0};
^]C&tG0 !
double ihxl[ia][JE][KE]={0},ihxh[ia][JE][KE]={0};
6i@\5}m=
double ihyl[IE][ja][KE]={0},ihyh[IE][ja][KE]={0};
2{,n_w?Wy
double ihzl[IE][JE][ka]={0},ihzh[IE][JE][ka]={0};
VaONd0Z I
double gi1[IE]={0},gi2[IE]={0},gi3[IE]={0};
v;1F[?@3Y
double gj1[JE]={0},gj2[JE]={0},gj3[JE]={0};
Q%S9fq,q
double gk1[KE]={0},gk2[KE]={0},gk3[KE]={0};
t ~7V { xk
double fi1[IE]={0},fi2[IE]={0},fi3[IE]={0};
R,2P3lv1v@
double fj1[JE]={0},fj2[JE]={0},fj3[JE]={0};
} YhtUWz].
double fk1[KE]={0},fk2[KE]={0},fk3[KE]={0};
W-~n|PX8+
double Ex_down[length][length][2]={0},Ey_down[length][length][2]={0};
ZAKeEm2A
double Hx_down[length][length][2]={0},Hy_down[length][length][2]={0};
25y6a| `
double Ex_up[length][length][2]={0},Ey_up[length][length][2]={0};
LtgXShp_!
double Hx_up[length][length][2]={0},Hy_up[length][length][2]={0};
rNOES3[~
double Ex_left[length][length][2]={0},Ez_left[length][length][2]={0};
:nY 2O
double Hx_left[length][length][2]={0},Hz_left[length][length][2]={0};
cgu~
double Ex_right[length][length][2]={0},Ez_right[length][length][2]={0};
G bP!9I
double Hx_right[length][length][2]={0},Hz_right[length][length][2]={0};
& BE g
double Ey_front[length][length][2]={0},Ez_front[length][length][2]={0};
0u B'g+MU`
double Hy_front[length][length][2]={0},Hz_front[length][length][2]={0};
e$) 300 o
double Ey_back[length][length][2]={0},Ez_back[length][length][2]={0};
H].y w9
double Hy_back[length][length][2]={0},Hz_back[length][length][2]={0};
vdDludEv
ic = IE/2;
uGN ^!NG-0
jc = JE/2;
* *H&+T/B
kc = KE/2;
,3i,P(?(
ib = IE - ia -1;
E; $+f
jb = JE - ja -1;
\00DqL(Oj`
kb = KE - ka -1;
Z"-L[2E/{!
//外推边界位置
;lX(}2tXW
int dis = 5;
{Lz H&qu
const int i_front = ib + dis;
po9 9 y-
const int i_back = ia - dis;
t(!r8!c u}
const int j_right = jb + dis;
PCrU<J 7
const int j_left = ja - dis;
nz.{P@[Qk
const int k_up = kb + dis;
BP[|n L
const int k_down = ka - dis;
&;T J~r#K
pi = 3.1415926;
^eYqll/U
epsz = 8.8e-12; // 真空介电常数
z&8un% Jt
epsz_r = 4.0; // 相对介电常数
_jU5O;
muz = 4*pi*1e-7; // 真空磁导率
07 g':QU@
muz_r = 1.0; // 相对磁导率
K5t0L!6<+
sigma = 0; // 介质电导率
yF._*9Q3hK
v = 3.0e8; // 波速
IeX^4 rc(
f = 9.375e9; // 频率 10GHz
`Mj}md;O"
lamda = v/f; // 波长 0.03m
bizTd
kk = 2*pi/lamda; // 波数
=DbY? Q<Q
Z = sqrt(muz/epsz); // 真空波阻抗 377
T'-FV
ddx = lamda/42; // 网格尺寸
gR^>3n'
dt = ddx/(2.0*v); // 时间步长
2S3 lsp5!
/* Initialize the arrays */
,|?CU r9Y
for ( i=0; i<IE; i++)
?(6m VyIe
{
hg@}@Wq\)
for ( j=0; j<JE; j++)
4\ c,)U}
{
h0vob_Fdl
for ( k=0; k<KE; k++)
fU3`v\X
{
Ty(@+M~-
ex
[j][k] = 0.0;
QYb33pN|
ey
[j][k] = 0.0;
1SIq[1
ez
[j][k] = 0.0;
S8Fmy1#
dx
[j][k] = 0.0;
RkeltE~u
dy
[j][k] = 0.0;
V D?*h
dz
[j][k] = 0.0;
)D-c]+yt
hx
[j][k] = 0.0;
#zUX yT#X
hy
[j][k] = 0.0;
C Wo1.pV w
hz
[j][k] = 0.0;
NG&_?| OmV
ix
[j][k] = 0.0;
.9[45][FK
iy
[j][k] = 0.0;
M6r^L6$N
iz
[j][k] = 0.0;
H5cV5E0
gax
[j][k] = 1.0;
-D^A: }$
gay
[j][k] = 1.0;
%?C8mA'w
gaz
[j][k] = 1.0;
m(iR|Zx
gbx
[j][k] = 0.0;
g?v\!/~(u
gby
[j][k] = 0.0;
_Us*+ 2(4L
gbz
[j][k] = 0.0;
S/ywA9~3Q
}
X .F^$
}
M]v= -
}
wXjFLg!g?
for ( i=0; i<ia; i++)
gd/H``x|Y
{
?rX]x8iP
for ( j=0; j<JE; j++)
?YM4b5!3T
{
6(|d|Si *c
for ( k=0; k<KE; k++)
nR>r2wMk@
{
6 S i-u
idxl
[j][k] = 0.0;
;^Sr"v6r>u
idxh
[j][k] = 0.0;
~rr 4o k
ihxl
[j][k] = 0.0;
ysIh[1E~%:
ihxh
[j][k] = 0.0;
= VLS/\ A
}
|wE3UWsy
}
x3ERCqTR
}
?Ucu#UO
for ( i=0; i<IE; i++)
cV{%^0? D
{
]j~"mFAP
for ( j=0; j<ja; j++)
}V;+l8
{
^I6^ g
for ( k=0; k<KE; k++)
4F3x@H'
{
Z BUArIC
idyl
[j][k] = 0.0;
wu9=N ^x
idyh
[j][k] = 0.0;
sAJ7R(p
ihyl
[j][k] = 0.0;
f&, {XZ
ihyh
[j][k] = 0.0;
)&{K~i ;:
}
dBV7Te4L
}
5F|8?BkOL^
}
\qvaE +
for ( i=0; i<IE; i++)
x_!0.SU
{
{+<P:jbz;
for ( j=0; j<JE; j++)
6 &Juv
{
fEE /-}d
for ( k=0; k<ka; k++)
88"Sai
{
]{~NO{0@Y
idzl
[j][k] = 0.0;
P |2E2=G
idzh
[j][k] = 0.0;
`cQo0{xK
ihzl
[j][k] = 0.0;
&M/0 g]4p
ihzh
[j][k] = 0.0;
gk BdR +
}
:HxA`@Ok
}
IL /Y c1
}
pC pb;<JG
//PML初始化
+ 1IQYa|
for ( i=0; i<IE; i++)
;;$# )b
{
ZD#9&q'4<
gi1
= 0;
R1:k23 {
fi1
= 0;
8n BL\{'B[
gi2
= 1;
Xb !MaNm)
fi2
= 1;
7?gFy-
gi3
= 1;
tj]9~eJ-
fi3
= 1;
MF3b{|Z
}
Lj\/ Ji_
for ( j=0; j<JE; j++)
-V"22sR]
{
d%I" /8-J
gj1[j] = 0;
Ch()P.n?
fj1[j] = 0;
}0`nvAf
gj2[j] = 1;
Sw`RBN[ yo
fj2[j] = 1;
_sE#)@p
gj3[j] = 1;
0n ~ Zz
fj3[j] = 1;
/WV7gO&L1
}
; R=.iOn
for ( k=0; k<KE; k++)
(C] SH\
{
@pI5lh
gk1[k] = 0;
H xV#WoYKj
fk1[k] = 0;
B<8N96fx
gk2[k] = 1;
`dMqe\o%!
fk2[k] = 1;
:M{Y,~cP
gk3[k] = 1;
m\>a,oZH
fk3[k] = 1;
&bz:K8c
}
GhY1k";
radius = 10 ;
c H-@V<
n_pml = 7; //PML层厚度
%u^ J pC{E
npml = 7;
@=]~\[e\
for (i = 0; i < n_pml; i++)
M aBYk?TR~
{
V~ ~=Qp+.
xxn = (double)(npml-i)/(double)npml;
b?4/#&z]
xn = 0.33*pow(xxn,3);
r dnno
fi1
= xn;
1QZ&Mj^^
fi1[IE-1-i] = xn;
-"Y{$/B
gi2
= 1.0/(1.0 + xn);
o[ENp'r
gi2[IE-1-i] = 1.0/(1.0 + xn);
]u-]'P
gi3
= (1.0 - xn)/(1.0 + xn);
Ed 0I WPx
gi3[IE-1-i] = (1.0 - xn)/(1.0 + xn);
0yx 3OY
xxn = (double)(npml-i-.5)/(double)npml;
Ee1LO#^_6
xn = .33*pow(xxn,3);
M F6 0-VE
gi1
= xn;
BYW^/B Y)
gi1[IE-2-i] = xn;
@ ''GPL@
fi2
= 1.0/(1.0 + xn);
G iq=*D+
fi2[IE-2-i] = 1.0/(1.0 + xn);
b96%")
fi3
= (1.0 - xn)/(1.0 + xn);
Bx"7%[
fi3[IE-2-i] = (1.0 - xn)/(1.0 + xn);
VN0mDh?E
}
to8X=80-3
for (j = 0; j < n_pml; j++)
YI-O{U
{
i`/+, <
xxn = (double)(npml-j)/(double)npml;
2>-S-;i
xn = .33*pow(xxn,3);
=K&\E2kA 4
fj1[j] = xn;
d w~p? [
fj1[JE-1-j] = xn;
dt[k\ !-v
gj2[j] = 1.0/(1.0 + xn);
m|=Ecu
gj2[JE-1-j] = 1.0/(1.0 + xn);
N$Y " c*
gj3[j] = (1.0 - xn)/(1.0 + xn);
IC{e E
gj3[JE-1-j] = (1.0 - xn)/(1.0 + xn);
/MB{Pmk$R
xxn = (double)(npml-j-.5)/(double)npml;
;n=. {[,
xn = .33*pow(xxn,3);
+5 .t. d
gj1[j] = xn;
,< <4*
gj1[JE-2-j] = xn;
PN~@
fj2[j] = 1.0/(1.0 + xn);
; U: {/
fj2[JE-2-j] = 1.0/(1.0 + xn);
h=kQ$`j6
fj3[j] = (1.0 - xn)/(1.0 + xn);
Sp}D ;7
fj3[JE-2-j] = (1.0 - xn)/(1.0 + xn);
6OiSK@<Hk
}
vQF vtwd
Iu~\L0R427
for (k = 0; k < n_pml; k++)
M L7 \BT
{
Sx*oo{Kk%
xxn = (double)(npml-k)/(double)npml;
&bgvy'p
xn = .33*pow(xxn,3);
!|2VWI}
fk1[k] = xn;
v{i 7h|e
fk1[KE-1-k] = xn;
;o/>JHGj
gk2[k] = 1.0/(1.0 + xn);
5!C_X5M
gk2[KE-1-k] = 1.0/(1.0 + xn);
2M)]!lYy
gk3[k] = (1.0 - xn)/(1.0 + xn);
7@.cOB`y@3
gk3[KE-1-k] = (1.0 - xn)/(1.0 + xn);
obK*rdg ,
xxn = (double)(npml-k-.5)/(double)npml;
zJ+8FWy:S
xn = .33*pow(xxn,3);
wdzOFDA
gk1[k] = xn;
'`Bm'Dd
gk1[KE-2-k] = xn;
Bb& ^ {7
fk2[k] = 1.0/(1.0 + xn);
)CI1;
fk2[KE-2-k] = 1.0/(1.0 + xn);
&2 -L. Xb
fk3[k] = (1.0 - xn)/(1.0 + xn);
wiOgyMdx
fk3[KE-2-k] = (1.0 - xn)/(1.0 + xn);
s{ V*1$e~
}
? Y:x[pOe
5dx&Qu'}ZS
/* Calculate gax ,gbx */
4GG>!@|
for ( i = ia; i<=ib; i++)
AeIrr*~]B
{
C^t(^9
for ( j = ja; j<=jb; j++)
PqNFy Qkl
{
RUC PV[{b
for ( k =ka; k<=kb; k++)
sfr(/mp(
{
R j~
xdist = (ic-i-.5);
PCd0 ?c
ydist = (jc-j);
x# ~ x ;)
zdist = (kc-k);
@ P@c.*}s
dist = sqrt (pow(xdist,2) + pow(ydist,2) + pow(zdist,2));
@ZN^1?][
if ( dist <= radius )
~zm 7?_"@]
{
Mh(]3\
gax
[j][k] = 1./(epsz_r + (sigma*dt/epsz));
,1q_pep~?%
gbx
[j][k] = sigma*dt/epsz;
D3Jr3 %>
}
k^$+n _
}
,=(Z00#(
}
FZd.L6q
}
SUWD]k >PH
/* Calculate gay ,gby*/
n[Zz]IO,g
for ( i = ia; i<=ib; i++)
sDg XU@
{
oHM ]
for ( j = ja; j<=jb; j++)
|vte =)%
{
syx\gz
for ( k =ka; k<=kb; k++)
r _,_5 @0e
{
M:E r_,E
xdist = (ic-i);
?^' 7+8C*J
ydist = (jc-j-.5);
K.C> a:J
zdist = (kc-k);
a/~aFmu6b
dist = sqrt (pow(xdist,2) + pow(ydist,2) + pow(zdist,2));
oRN-xn g
if ( dist <= radius )
nqR?l4 DX
{
9[v1h,L
gay
[j][k] = 1./(epsz_r + (sigma*dt/epsz));
lX:| iB
gby
[j][k] = sigma*dt/epsz;
[w -{r+[
}
TUL_TR
}
UH@ a s
}
*KK+X07
}
2c<&eX8 "
/* Calculate gaz ,gbz */
k>}g\a,
for ( i = ia; i<=ib; i++)
jk\ dG16
{
QU{Ech'
for ( j = ja; j<=jb; j++)
UX`]k{Mz
{
<1w/hy&mWN
for ( k =ka; k<=kb; k++)
0U66y6
{
~5ZvOX6L2
xdist = (ic-i);
)oo~m\`
ydist = (jc-j);
d#:3be{|&q
zdist = (kc-k-.5);
H-rWDN#
dist = sqrt (pow(xdist,2) + pow(ydist,2) + pow(zdist,2));
eRbGZYrJ
if ( dist <= radius )
"SN4*
{
q,<[hBri-
gaz
[j][k] = 1./(epsz_r +(sigma*dt/epsz));
@R D+xYm
gbz
[j][k] = sigma*dt/epsz;
GwsY-jf
}
?3z- _8#
}
E< 4l#Z<
}
) T907I|
}
*IlaM' [*
//源参数
JXiZB 8}
t0 = 40.0;
w)E@*h< Z
spread = 10.0;
3>zN/ f
[MX;,%;;
clock_t start, finish;
D0Z\Vvy
start = clock();
{v a aFs
/* outFile.open("F:\\ez[40][25][25].txt");*/
{"*VU3%q
for ( T = 1; T < 1001; T++ )
L+B?~_*
{
j^`X~gE
/* Calculate the incident buffer */
m,3er*t{
for ( j=1; j<JE; j++)
^x*nq3^h\
{
Uu6L~iB
ez_inc[j] = ez_inc[j] + .5*(hx_inc[j-1] - hx_inc[j]);
D/(L
}
?:^mBb) T
//Fourier Transform of the incident field
M_B:{%4
ine[0] = ine[0] + dt*cos(2*pi*f*dt*T)*ez_inc[ja-1];
-7WW[ w
ine[1] = ine[1] - dt*sin(2*pi*f*dt*T)*ez_inc[ja-1];
w&Dv8Wv+Oq
/* Source */
mtic>
//pulse = exp (-.5*(pow(( t0-T )/spread,2)));
@J[6,$UVu
pulse = exp (2*pi*f*T*dt);
"wH) mQnd
ez_inc[3] = pulse;
< (RC|?
/* Boundary conditions for the incident buffer */
SEQ%'E5-'
ez_inc[0] = ez_low_m2;
8$ xPex~2
ez_low_m2 = ez_low_m1;
#L crI
ez_low_m1 = ez_inc[1];
e"-X U@`k1
ez_inc[JE-1] = ez_high_m2;
[\)oo
ez_high_m2 = ez_high_m1;
#Z =tJ
ez_high_m1 = ez_inc[JE-2];
K*K1(_x=
/* Calculate the Dx field */
nGZX7Fx5
for ( i=1; i<ia; i++)
LCXO>MXN
{
k'`m97B
for ( j=1; j<JE; j++)
<Y 4:'L6
{
;S$Ll*f>D
for (k=1; k<KE; k++)
wYeB)1.
{
di9!lS $
curl_h = ( hz
[j][k] - hz
[j-1][k] - hy
[j][k] + hy
[j][k-1] );
MFJE6ei
idxl
[j][k] = idxl
[j][k] + curl_h;
'f8(#n=6qP
dx
[j][k] = gj3[j]*gk3[k]*dx
[j][k] + gj2[j]*gk2[k]*.5*(curl_h + gi1
*idxl
[j][k]);
Y izE5 [*
}
Rg29
}
!=Y;h[J.p
}
n9LGP2#!
for ( i=ia; i<=ib; i++)
fnzy5+9"
{
O|kKwadC
for ( j=1;j<JE; j++)
)j6eE+gF
{
oC*ees g_
for ( k=1; k<KE; k++)
?<X(]I.j
{
0m[dP
curl_h = ( hz
[j][k] - hz
[j-1][k] - hy
[j][k] + hy
[j][k-1] );
Kj=gm .
dx
[j][k] = gj3[j]*gk3[k]*dx
[j][k] + gj2[j]*gk2[k]*.5*curl_h;
fi PIAT}
}
{ PlK@#UN
}
[ {[N( g&d
}
UbJ*'eoX
for ( i=ib+1; i<IE; i++)
:(3|HTz
{
k!jNOqbb
ixh = i-ib-1;
U IJx*
for ( j=1; j<JE; j++)
{hSGv
{
5FvOznK^e
for ( k=1; k<KE; k++)
kg5ev8
{
Mro4`GL
curl_h = ( hz
[j][k] - hz
[j-1][k] - hy
[j][k] + hy
[j][k-1] );
k>4qkigjc
idxh[ixh][j][k] = idxh[ixh][j][k] + curl_h;
dtM[E`PL
dx
[j][k] = gj3[j]*gk3[k]*dx
[j][k] + gj2[j]*gk2[k]*.5*(curl_h + gi1
*idxh[ixh][j][k]);
;SwC&.I
}
Pvo#pY^dXX
}
|wxGpBau
}
{mYP<N BT
/* Calculate the Dy field */
[y) Fc IK}
for ( i=1; i<IE; i++)
78#!Q.##
{
i;\s.wrzH
for ( j=1; j<ja; j++)
=<f-ob8,
{
g]mtFrP
for (k=1; k<KE; k++)
h SqY$P
{
~%w~ - O2
curl_h = ( hx
[j][k] - hx
[j][k-1] - hz
[j][k] + hz[i-1][j][k] );
}pNX@C#De
idyl
[j][k] = idyl
[j][k] + curl_h;
$KV&\Q3\0
dy
[j][k] = gi3
*gk3[k]*dy
[j][k]+gi2
*gk2[k]*.5*(curl_h+ gj1[j]*idyl
[j][k]);
O/A E}]
}
BJjx|VA+
}
"yb WD Wu
}
4YG/` P
for ( i=1; i<IE; i++)
/k\01hc`
{
TDK@)mP
for ( j=ja;j<=jb; j++)
!&kL9A).
{
ZU'!iU|8
for ( k=1; k<KE; k++)
hi{%pi&!T
{
/z."l!u6
curl_h = ( hx
[j][k] - hx
[j][k-1] - hz
[j][k] + hz[i-1][j][k] );
YZ#V#[j'^
dy
[j][k] = gi3
*gk3[k]*dy
[j][k]+gi2
*gk2[k]*.5*curl_h;
qcB){p+UQ
}
S^RUw
}
L6:h.1 U$
}
iMP
for ( i=1; i<IE; i++)
noVa=aU^
{
`ue[q!Qq
for ( j=jb+1; j<JE; j++)
) yee2(S
{
S[M4ukYK
jyh = j-jb-1;
Z%o7f6P0IX
for ( k=1; k<KE; k++)
Lrz3
{
'3tw<k!1{.
curl_h = ( hx
[j][k] - hx
[j][k-1] - hz
[j][k] + hz[i-1][j][k] );
-Q e~)7
idyh
[jyh][k] = idyh
[jyh][k] + curl_h;
d<p 2/ aA
dy
[j][k] = gi3
*gk3[k]*dy
[j][k]+gi2
*gk2[k]*.5*(curl_h+ gj1[j]*idyh
[jyh][k]);
O0l^*nZ46t
}
/4S;QEv
}
^~ =9
}
jdRq6U^
/* Incident Dy */
j28 _Hh T
for ( i=ia; i<=ib; i++)
8Vl!| \x5
{
i`U: gw
for ( j=ja; j<=jb-1; j++)
9 wa,k
{
%(m ])
dy
[j][ka] = dy
[j][ka] - .5*hx_inc[j];
q1Qje%9@t
dy
[j][kb+1] = dy
[j][kb+1] + .5*hx_inc[j];
u/^|XOy
}
}am U[U,
}
V*n==Nb5L
/* Calculate the Dz field */
B l.u=I:Y4
for ( i=1; i<IE; i++)
dZ2 %S''\
{
`{<f rB@
for ( j=1; j<JE; j++)
_]#klL
{
I8m:3fL"
for (k=1; k<ka; k++)
+`en{$%%
{
#mc!Wt 10
curl_h = ( hy
[j][k] - hy[i-1][j][k] - hx
[j][k] + hx
[j-1][k] );
^i%A7pg
idzl
[j][k] = idzl
[j][k] + curl_h;
*DeTqO65
dz
[j][k] = gi3
*gj3[j]*dz
[j][k]+gi2
*gj2[j]*.5*(curl_h+ gk1[k]*idzl
[j][k]);
mL L$|
}
N$aZ== $5
}
ND]S(C"?
}
~&G4)AM
for ( i=1; i<IE; i++)
,4F,:w
{
G"s0GpvQ
for ( j=1;j<JE; j++)
64ox jF)
{
a_T,t' 6
for ( k=ka; k<=kb; k++)
M8W# io
{
7Ur?e p
curl_h = ( hy
[j][k] - hy[i-1][j][k] - hx
[j][k] + hx
[j-1][k] );
^tAO _ ~4
dz
[j][k] = gi3
*gj3[j]*dz
[j][k]+gi2
*gj2[j]*.5*curl_h;
-/{af
}
CDg AGy
}
,TC~~EWq
}
gy _$#e
for ( i=1; i<IE; i++)
t\y-T$\\
{
e$l 6gY
for ( j=1; j<JE; j++)
``4wX-y
{
Rq)BssdF
for ( k=kb+1; k<KE; k++)
_g|acBF
{
a<Uqyilm
kzh = k-kb-1;
h* .w"JO
curl_h = ( hy
[j][k] - hy[i-1][j][k] - hx
[j][k] + hx
[j-1][k] );
DQ6jT@ZDH
idzh
[j][kzh] = idzh
[j][kzh] + curl_h;
;_JH:}j
dz
[j][k] = gi3
*gj3[j]*dz
[j][k]+gi2
*gj2[j]*.5*(curl_h+ gk1[k]*idzh
[j][kzh]);
M]_E
}
s${_K* g6
}
gp
}
[, f)9v)
/* Incident Dz */
)q^ Bj$
for ( i=ia; i<=ib; i++)
!e8i /!}^S
{
] pPz@@xx
for ( k=ka; k<=kb; k++)
90:K#nW;
{
crmUrF#
dz
[ja][k] = dz
[ja][k] + .5*hx_inc[ja-1];
D&{ *AH%Q
dz
[jb][k] = dz
[jb][k] - .5*hx_inc[jb];
vE\lp8 j+
}
<)#kq1b?
}
q^Tis>*u6
/* Calculate the E from D field */
# mT]j""
for ( i=1; i<IE-1; i++)
3{J.xWB@:
{
hY8#b)l~lu
for ( j=1; j<JE-1; j++)
iyM^[/-R6
{
3;/?q
for ( k=1; k<KE-1; k++)
qc8Ta "
{
uDG+SdyN@
ex
[j][k] = gax
[j][k]*(dx
[j][k]-ix
[j][k]);
h8}8Lp(/'
ix
[j][k] = ix
[j][k]+gbx
[j][k]*ex
[j][k];
(Pi-uL<[a
ey
[j][k] = gay
[j][k]*(dy
[j][k]-iy
[j][k]);
X-Ycz 5?
iy
[j][k] = iy
[j][k]+gby
[j][k]*ey
[j][k];
>9 (i)e
ez
[j][k] = gaz
[j][k]*(dz
[j][k]-iz
[j][k]);
Wc'Ehyi;
iz
[j][k] = iz
[j][k]+gbz
[j][k]*ez
[j][k];
?b$3ob"
}
L7q | ^`
}
k$[{n'\@
}
&TT":FPR
/* Remember: part of the PML is E = 0 at the edges */
Aka^e\Y@6*
for (j = 0;j < JE-1;j++)
-CBD|fo[h
{
D1"7s,Hmu
for (k = 0;k < KE-1;k++)
F,}s$ v
{
X , Z eD
ez[0][j][k] = 0.0;
jM U9{Si
ez[IE-1][j][k] = 0.0;
b~Z=:'m8
}
"DckwtG :%
}
W cGXp $M
for (i = 0;i < IE-1;i++)
J/Q|uRpmqr
{
H?:Jq\Ba0
for (k = 0;k < KE;k++)
#ooc)) ,
{
U</+ .$b
ez
[0][k] = 0.0;
!Yh}H<w0
ez
[JE-1][k] = 0.0;
jxY-u+B
}
:m`/Q_y"
}
9 81!2*
for (i = 0;i < IE-1;i++)
h5onRa *7
{
{+=i?
for (j = 0;j < JE-1;j++)
~77 5soN
{
J?jeYW
ez
[j][0] = 0.0;
j[DIz@^
ez
[j][KE-1] = 0.0;
_Y-$}KwY!
}
YFx=b!/ s
}
:C%47qv
9*_uCPR
/*outFile << setw(22) << setprecision(8) << T << setw(22) << setprecision(8) << ez[40][25][25]<< endl;*/
!<Z{@7oH
YNM\pX'
/* Calculate the incident field */
q|{tQJfYg
for ( j=0; j<JE-1; j++)
u%opY<h
{
dF,DiRD
hx_inc[j] = hx_inc[j] + .5*(ez_inc[j]-ez_inc[j+1]);
D@ hmO]5c
}
87 8tI3-
/* Calculate the Hx field */
<xF?~7
for ( i=0; i<ia; i++)
`Cj,HI_/*
{
M {x ie
for ( j=0; j<JE-1; j++)
FmA-OqEpA
{
3 jh|y,
for (k=0; k<KE-1; k++)
q<XcOc5
{
$FCw$ +w
curl_e = ( ey
[j][k+1] - ey
[j][k] - ez
[j+1][k] + ez
[j][k] );
G%OpO.Wf
ihxl
[j][k] = ihxl
[j][k] + curl_e;
. bG{T|
hx
[j][k] = fj3[j]*fk3[k]*hx
[j][k]+fj2[j]*fk2[k]*.5*(curl_e+ fi1
*ihxl
[j][k]);
C4_t_N
}
;D&wh
}
faVS2TN4
}
6y@<?08Q
for ( i=ia; i<=ib; i++)
*|'}v[{v^9
{
oNSz& )LP
for ( j=0;j<JE-1; j++)
d| \#?W&
{
>B~?dT m
for ( k=0; k<KE-1; k++)
OK\] *r
{
dofR)"<p,^
curl_e = ( ey
[j][k+1] - ey
[j][k] - ez
[j+1][k] + ez
[j][k] );
LXxl ?D
hx
[j][k] = fj3[j]*fk3[k]*hx
[j][k]+fj2[j]*fk2[k]*.5*curl_e;
z(UX't (q
}
6D^%'[4t
}
Q-Y@)Mf~?0
}
c~oe, 9
for ( i=ib+1; i<IE; i++)
[+DNM 2A
{
1UyH 0`&
ixh = i-ib-1;
)dg U mN
for ( j=0; j<JE-1; j++)
UCF[oO>v
{
<4NQL*|>
for ( k=0; k<KE-1; k++)
>@Khm"/T
{
A I fk"2
curl_e = ( ey
[j][k+1] - ey
[j][k] - ez
[j+1][k] + ez
[j][k] );
}}s) +d
ihxh[ixh][j][k] = ihxh[ixh][j][k] + curl_e;
'%O\ E{h
hx
[j][k] = fj3[j]*fk3[k]*hx
[j][k]+fj2[j]*fk2[k]*.5*(curl_e+ fi1
*ihxh[ixh][j][k]);
nDn{zea7
}
y Hw!#gWM
}
t^$Div_%G
}
^ K>pT}u
/* Incident Hx */
*]9XDc]{j1
for ( i=ia; i<=ib; i++)
m?vAyi
{
= +Tsknq
for ( k=ka; k<=kb; k++)
f[M"EMy
{
J WS q"N
hx
[ja-1][k] = hx
[ja-1][k] + .5*ez_inc[ja];
2o}8W7y
hx
[jb][k] = hx
[jb][k] - .5*ez_inc[jb];
Vx(;|/:
}
R7t bxC
}
UJs?9]x>
Bcm=G""
/* Calculate the Hy field */
dh,7 iQ s
for ( i=0; i<IE-1; i++)
<Am^z~[
{
.Bijc G
for ( j=0; j<ja; j++)
/2'c>
{
<^'{ G
for (k=0; k<KE-1; k++)
_ ^3@PM>
{
B#4 S/d{/
curl_e = ( ez[i+1][j][k] - ez
[j][k] - ex
[j][k+1] + ex
[j][k] );
]eL~L_[G\
ihyl
[j][k] = ihyl
[j][k] + curl_e;
XAlD ww
hy
[j][k] = fi3
*fk3[k]*hy
[j][k]+fi2
*fk2[k]*.5*(curl_e+ fj1[j]*ihyl
[j][k]);
<M>#qd@c
}
[r~~=b7*[
}
A>y#}^l]
}
<on)"{W13
for ( i=0; i<IE-1; i++)
iy9VruT< x
{
c$.T<r)Z
for ( j=ja;j<=jb; j++)
/K& wr6
{
?(M\:`G'
for ( k=0; k<KE-1; k++)
Es,0'\m&
{
L`bo#,eg6
curl_e = ( ez[i+1][j][k] - ez
[j][k] - ex
[j][k+1] + ex
[j][k] );
+{bh
hy
[j][k] = fi3
*fk3[k]*hy
[j][k]+fi2
*fk2[k]*.5*curl_e;
u"joCZ7`kG
}
<v -92?
}
.=aMjrME
}
8Nv-/VQ/b
for ( i=0; i<IE-1; i++)
#|D:f~"d3
{
h@/>?Va
for ( j=jb+1; j<JE; j++)
<sF!]R&4
{
^)'| |Ly
jyh = j-jb-1;
.P x,=56$X
for ( k=0; k<KE-1; k++)
n]Ebwznt-
{
B C&^]M
curl_e = ( ez[i+1][j][k] - ez
[j][k] - ex
[j][k+1] + ex
[j][k] );
vKV{ $|
ihyh
[jyh][k] = ihyh
[jyh][k] + curl_e;
k dUc&
hy
[j][k] = fi3
*fk3[k]*hy
[j][k]+fi2
*fk2[k]*.5*(curl_e+ fj1[j]*ihyh
[jyh][k]);
1 pYsjo~
}
s-o~@(r6
}
~M(pCSJ[
}
)%1&/uN)
/* Incident Hy */
$o: :PDQ?
for ( j=ja; j<=jb; j++)
pe|\'<>i
{
N}1- 2
for ( k=ka; k<=kb; k++)
.;}pU!S~R
{
gGD]t;<u
hy[ia-1][j][k] = hy[ia-1][j][k] - .5*ez_inc[j];
6q!7i%fK?
hy[ib][j][k] = hy[ib][j][k] + .5*ez_inc[j];
@~Rk^/0
}
r6/<&1[
}
UvoG< ;
/* Calculate the Hz field */
PK9Qm'W b
for ( i=0; i<IE-1; i++)
QT JrJD
{
;+! xZOmm
for ( j=0; j<JE-1; j++)
%hrsE5k^,
{
%k ~C -+
for (k=0; k<ka; k++)
@V*au:
{
ec,z6v^9
curl_e = ( ex
[j+1][k] - ex
[j][k] - ey[i+1][j][k] + ey
[j][k] );
X]>[Qz)K^
ihzl
[j][k] = ihzl
[j][k] + curl_e;
( rMTW+,
hz
[j][k] = fi3
*fj3[j]*hz
[j][k]+fi2
*fj2[j]*.5*(curl_e+ fk1[k]*ihzl
[j][k]);
Oym]&SrbS
}
F\l!A'Q+t
}
ROcY'-
}
j= p|'`
for ( i=0; i<IE-1; i++)
9 ?[4i'
{
$::51# ^Wg
for ( j=0;j<JE-1; j++)
047*gn.b
{
"otr+.{`*
for ( k=ka; k<=kb; k++)
il<gjlyR]L
{
[ )B@
curl_e = ( ex
[j+1][k] - ex
[j][k] - ey[i+1][j][k] + ey
[j][k] );
I%C]>ZZh
hz
[j][k] = fi3
*fj3[j]*hz
[j][k]+fi2
*fj2[j]*.5*curl_e;
qrOB_Nz
}
$4& 8U ~Zs
}
N>h]mX6
}
6NP`P j R
for ( i=0; i<IE-1; i++)
!G@V<'F
{
~p0 e=u
for ( j=0; j<JE-1; j++)
_y.mpX&
{
nT(!HDH
for ( k=kb+1; k<KE; k++)
:^C#-O
{
30:HRF( :
kzh = k-kb-1;
%YsRm%q
curl_e = ( ex
[j+1][k] - ex
[j][k] - ey[i+1][j][k] + ey
[j][k] );
U5T^S
ihzh
[j][kzh] = ihzh
[j][kzh] + curl_e;
oKZ [0(4<
hz
[j][k] = fi3
*fj3[j]*hz
[j][k]+fi2
*fj2[j]*.5*(curl_e+ fk1[k]*ihzh
[j][kzh]);
