登 录
註 冊
论坛
微波仿真网
注册
登录论坛可查看更多信息
微波仿真论坛
>
时域有限差分法 FDTD
>
大家帮忙看看这个3D_FDTD介质球的RCS有什么问 ..
发帖
回复
1
2
3
4122
阅读
29
回复
大家帮忙看看这个3D_FDTD介质球的RCS有什么问题
离线
cc81372365
UID ：27752
注册：
2009-03-18
登录：
2023-11-09
发帖：
345
等级：
论坛版主
0楼
发表于: 2010-11-29 13:16:05
//****************************************************************//
Wk}D]o0^@
//
E`tQe5K
//********* 3D_介质球散射（RCS）FDTD程序 **********//
FZpsL-yx^N
//
PRwu
//****************************************************************//
*qG $19b
//
"c,!v c4
//入射波为Y方向的平面波（theta=pi/2 , phi=pi/2), RCS观察角为（theta=pi/2 , phi = 0到2pi）
qC..\{z
// 数字0，1代表场值经过傅里叶变换后的复数的实部和虚部
WO@H*
// 金属方柱置于中心位置（ic,jc）,激励源为高斯脉冲源
~5ubh 2{
#include <iostream>
ywEDy|Wn$~
#include <iomanip>
|YRY!V_w
#include <fstream>
l DnMjK\M
//#include <omp.h>
D M}s0O$ 0
#include <time.h>
7 W{~f? Sh
#include <valarray>
o!utZmk$
#include <math.h>
SF^x=[ir
#include <stdio.h>
8)Zk24:])_
#include <stdlib.h>
Y9tV%
n$YE !D'
#include "complex_operation.h"
|Y tg
#define IE 60
H*rx{ F?
#define JE 60
7b46 t2W<
#define KE 60
[fxuUmU
#define ia 15
b_l3+'#ofM
#define ja 15
<5I1 DF[
#define ka 15
Ft) lp>3gv
#define length 41
jN6b*-2
using namespace std;
HlPG3LD!
ofstream outFile;
P x#$uU
int main()
<m0{'xw
{
YDF CGA
int i,j,k,ic,jc,kc,ib,jb,kb;
,gmH2.
int npml,n_pml,T;
/[iG5~G
double radius,v,lamda,Z,f,kk;
zP&q7 t;>
double epsz_r,epsz,muz,muz_r,ddx,dt,pi,sigma;
TP{Gt.e
double xn,xxn;
88g47>{X
double t0,spread,pulse;
/l6\^Xf{
double ez_inc[JE]={0},hx_inc[JE]={0};
rrcwtLNb u
double ez_low_m1=0,ez_low_m2=0,ez_high_m1=0,ez_high_m2=0;
.H2qs{N!
int ixh,jyh,kzh;
?+Q$#pb
int co=0;
74_xR
double rcs[361]={0.0};
_ 88QgThb
double ine[2] = {0.0};
Gqt-_gga
double ine_fudu = 0.0;
.LObOR 5J7
double curl_h,curl_e;
\?&A u
double dist,xdist,ydist,zdist;
, 1 P[
double dx[IE][JE][KE]={0},dy[IE][JE][KE]={0},dz[IE][JE][KE]={0};
bDWe U}
double ex[IE][JE][KE]={0},ey[IE][JE][KE]={0},ez[IE][JE][KE]={0};
U(~U!O}
double hx[IE][JE][KE]={0},hy[IE][JE][KE]={0},hz[IE][JE][KE]={0};
/$:U$JVb?l
double ix[IE][JE][KE]={0},iy[IE][JE][KE]={0},iz[IE][JE][KE]={0};
_q)!B,y-/N
double gax[IE][JE][KE]={0},gay[IE][JE][KE]={0},gaz[IE][JE][KE]={0};
"yW&<7 u1
double gbx[IE][JE][KE]={0},gby[IE][JE][KE]={0},gbz[IE][JE][KE]={0};
l{5O5%\,
double idxl[ia][JE][KE]={0},idxh[ia][JE][KE]={0};
HIGNRm
double idyl[IE][ja][KE]={0},idyh[IE][ja][KE]={0};
vbp-`M(
double idzl[IE][JE][ka]={0},idzh[IE][JE][ka]={0};
9 mPIykAj8
double ihxl[ia][JE][KE]={0},ihxh[ia][JE][KE]={0};
E/mw* c^
double ihyl[IE][ja][KE]={0},ihyh[IE][ja][KE]={0};
Mlj# b8
double ihzl[IE][JE][ka]={0},ihzh[IE][JE][ka]={0};
%b =p< h'(
double gi1[IE]={0},gi2[IE]={0},gi3[IE]={0};
5V @&o`!=h
double gj1[JE]={0},gj2[JE]={0},gj3[JE]={0};
VFSz-<L
double gk1[KE]={0},gk2[KE]={0},gk3[KE]={0};
g~ !$i`_b
double fi1[IE]={0},fi2[IE]={0},fi3[IE]={0};
*,lh:
double fj1[JE]={0},fj2[JE]={0},fj3[JE]={0};
ue6d~8&
double fk1[KE]={0},fk2[KE]={0},fk3[KE]={0};
W2 eAhz&
double Ex_down[length][length][2]={0},Ey_down[length][length][2]={0};
]'k[u
double Hx_down[length][length][2]={0},Hy_down[length][length][2]={0};
7~b!4x|Z
double Ex_up[length][length][2]={0},Ey_up[length][length][2]={0};
XsUUJu CG
double Hx_up[length][length][2]={0},Hy_up[length][length][2]={0};
Yj|]Uff8O
double Ex_left[length][length][2]={0},Ez_left[length][length][2]={0};
/C$ xH@bb
double Hx_left[length][length][2]={0},Hz_left[length][length][2]={0};
Z>2]Xx% \
double Ex_right[length][length][2]={0},Ez_right[length][length][2]={0};
:Q-oV8t{
double Hx_right[length][length][2]={0},Hz_right[length][length][2]={0};
MCU9O
double Ey_front[length][length][2]={0},Ez_front[length][length][2]={0};
G o <'
double Hy_front[length][length][2]={0},Hz_front[length][length][2]={0};
8bOT*^b$H
double Ey_back[length][length][2]={0},Ez_back[length][length][2]={0};
2]2H++
double Hy_back[length][length][2]={0},Hz_back[length][length][2]={0};
6>,# 6{?jl
ic = IE/2;
rl?7W];
jc = JE/2;
a4&:@ ` =
kc = KE/2;
nm @']
ib = IE - ia -1;
%!y89x=E
jb = JE - ja -1;
VE]6wwV2
kb = KE - ka -1;
.G CR!V
//外推边界位置
?4G(N=/&
int dis = 5;
JMlV@t7y<
const int i_front = ib + dis;
w!`Umll2
const int i_back = ia - dis;
RtC 'v";6
const int j_right = jb + dis;
[M:S`{SbY
const int j_left = ja - dis;
#mLuU
const int k_up = kb + dis;
ia4k :\
const int k_down = ka - dis;
9&K/GaG
pi = 3.1415926;
.N"~zOV<#
epsz = 8.8e-12; // 真空介电常数
(A{NF(
epsz_r = 4.0; // 相对介电常数
r5 yO 5W
muz = 4*pi*1e-7; // 真空磁导率
iF9_b
muz_r = 1.0; // 相对磁导率
mW4%2fD [
sigma = 0; // 介质电导率
m<: IFx#
v = 3.0e8; // 波速
}H y ~i
f = 9.375e9; // 频率 10GHz
PZ,z15PG]
lamda = v/f; // 波长 0.03m
;S JF%@x
kk = 2*pi/lamda; // 波数
Hm9< fQuM
Z = sqrt(muz/epsz); // 真空波阻抗 377
DMZ` Sx
ddx = lamda/42; // 网格尺寸
JE!("]&
dt = ddx/(2.0*v); // 时间步长
mp\%M 1<
/* Initialize the arrays */
B4yh3 cf
for ( i=0; i<IE; i++)
PAZ$_eSK6
{
7pet Hi
for ( j=0; j<JE; j++)
XmWlv{T+
{
i:Y \`J
for ( k=0; k<KE; k++)
</s,pe79B
{
}U[-44r:
ex
[j][k] = 0.0;
z[9UQU~x?
ey
[j][k] = 0.0;
jG D%r~lN
ez
[j][k] = 0.0;
*)> do L
dx
[j][k] = 0.0;
E>pVn2|
dy
[j][k] = 0.0;
kus}W J
dz
[j][k] = 0.0;
4 :1)~z
hx
[j][k] = 0.0;
2dbRE:v5
hy
[j][k] = 0.0;
.Yx_:h=u
hz
[j][k] = 0.0;
f=aIXhiYU
ix
[j][k] = 0.0;
#?&0D>E ?k
iy
[j][k] = 0.0;
'RT tE
iz
[j][k] = 0.0;
(oaYF+T
gax
[j][k] = 1.0;
^%#grX#
gay
[j][k] = 1.0;
;&f1vi4
gaz
[j][k] = 1.0;
M{G xjmdx
gbx
[j][k] = 0.0;
r]LCv sVa
gby
[j][k] = 0.0;
HZZD v+
gbz
[j][k] = 0.0;
2P9J' L
}
q`Q }yE> 9
}
}1 QF+C f
}
7[KCWJ
for ( i=0; i<ia; i++)
6RK\}@^=K
{
3z[ $4L'.
for ( j=0; j<JE; j++)
Q $a
{
saW!9HQj
for ( k=0; k<KE; k++)
[B9 ;?G
{
c tI{^f:
idxl
[j][k] = 0.0;
kuKa8c
idxh
[j][k] = 0.0;
B8V,)rn
ihxl
[j][k] = 0.0;
G!Zyl^
ihxh
[j][k] = 0.0;
qLjT.7 .x
}
<KQ(c`KW7
}
WRp yr
}
&[j]Bp?
for ( i=0; i<IE; i++)
`ZGcgO< c\
{
!wh&>3~
for ( j=0; j<ja; j++)
x(~< t X~
{
j*zD0I]
for ( k=0; k<KE; k++)
8&E}n(XE
{
OW`STp!
idyl
[j][k] = 0.0;
V'StvU
idyh
[j][k] = 0.0;
'M/ ([|@
ihyl
[j][k] = 0.0;
z~N a-N
ihyh
[j][k] = 0.0;
<DZcra
}
T~ k)uQ
}
nvyB/
}
9DE)S)e8
for ( i=0; i<IE; i++)
-]XP 2}#d
{
7SS07$B
for ( j=0; j<JE; j++)
X9XI;c;b-
{
$j*Qo/x d
for ( k=0; k<ka; k++)
'*!L!VJ
{
!buz<h
idzl
[j][k] = 0.0;
3M<!?%v\A
idzh
[j][k] = 0.0;
`# ^0cW
ihzl
[j][k] = 0.0;
Zdn! qyR`
ihzh
[j][k] = 0.0;
n&!+wcJ;Yt
}
SSmHEy*r)
}
JP'= UZ '
}
^o eJKjJ
//PML初始化
%Q4i%: Qi
for ( i=0; i<IE; i++)
ycX{NDGs
{
d`%M g&
gi1
= 0;
g7V_ [R(6
fi1
= 0;
1;wb(DN*c
gi2
= 1;
;n*J$B
fi2
= 1;
=2 jhII
gi3
= 1;
kGHC]Fb)
fi3
= 1;
|_zO_F rtp
}
q#[`KOPV
for ( j=0; j<JE; j++)
MR;X&Up6!
{
g>0vm2 |
gj1[j] = 0;
Lwd V3 vb#
fj1[j] = 0;
'Q|c@t
gj2[j] = 1;
=r6qX
fj2[j] = 1;
2 ZG@!Y|
gj3[j] = 1;
qnIew ?-*
fj3[j] = 1;
Fp VV4D
}
_T*A C.
for ( k=0; k<KE; k++)
{#hVD4$b
{
( ?Q|s,
gk1[k] = 0;
)bD nbO$s_
fk1[k] = 0;
3,6O x45
gk2[k] = 1;
I3) Zr+
fk2[k] = 1;
K=\&+at1
gk3[k] = 1;
(:sZ b?*
fk3[k] = 1;
+AI`R`Tm
}
U[|| ~FW'
radius = 10 ;
/[pqI0sf<A
n_pml = 7; //PML层厚度
FvXqggfGv
npml = 7;
<avQR 9'&
for (i = 0; i < n_pml; i++)
pP.'wSj
{
LW mB, Zf/
xxn = (double)(npml-i)/(double)npml;
hh"-w3+
xn = 0.33*pow(xxn,3);
wE Qi0!
fi1
= xn;
F ?=9eISLJ
fi1[IE-1-i] = xn;
Ai kf|)D[
gi2
= 1.0/(1.0 + xn);
Y25uU%6t_
gi2[IE-1-i] = 1.0/(1.0 + xn);
2\@ Z5m3B
gi3
= (1.0 - xn)/(1.0 + xn);
d2jr8U
gi3[IE-1-i] = (1.0 - xn)/(1.0 + xn);
m5-9yQ=.
xxn = (double)(npml-i-.5)/(double)npml;
(tCib 4
xn = .33*pow(xxn,3);
@ 4UxRp6+
gi1
= xn;
JN[0L :
gi1[IE-2-i] = xn;
|wp ,f%WK
fi2
= 1.0/(1.0 + xn);
srmKaa|
fi2[IE-2-i] = 1.0/(1.0 + xn);
4hAJ !7[A.
fi3
= (1.0 - xn)/(1.0 + xn);
PK:2xN:=
fi3[IE-2-i] = (1.0 - xn)/(1.0 + xn);
w^;DG
}
oXgdLtsu
for (j = 0; j < n_pml; j++)
5PiOH"!19
{
E zT`,#b
xxn = (double)(npml-j)/(double)npml;
sb.J bE8
xn = .33*pow(xxn,3);
OO-k|\{ |
fj1[j] = xn;
o 8S"&O ?
fj1[JE-1-j] = xn;
n=!uN u 7
gj2[j] = 1.0/(1.0 + xn);
2@Lb fo A
gj2[JE-1-j] = 1.0/(1.0 + xn);
' [p)N,
gj3[j] = (1.0 - xn)/(1.0 + xn);
+5X DF
gj3[JE-1-j] = (1.0 - xn)/(1.0 + xn);
.C!vr@@]
xxn = (double)(npml-j-.5)/(double)npml;
, 8VU&?`<}
xn = .33*pow(xxn,3);
k<Sl1v K
gj1[j] = xn;
Ge \["`;i
gj1[JE-2-j] = xn;
Oh; Jw
fj2[j] = 1.0/(1.0 + xn);
3W&S.$l
fj2[JE-2-j] = 1.0/(1.0 + xn);
8<dOMp;}r
fj3[j] = (1.0 - xn)/(1.0 + xn);
FFcB54ALTf
fj3[JE-2-j] = (1.0 - xn)/(1.0 + xn);
.SS<MDcqIt
}
HEY4$Lf(I
Yl({)qK{
for (k = 0; k < n_pml; k++)
D>#l -{d
{
;y1Q6eN
xxn = (double)(npml-k)/(double)npml;
A lwtmDa
xn = .33*pow(xxn,3);
%SA!p;
fk1[k] = xn;
#R5U
fk1[KE-1-k] = xn;
Ypm Yxd^
gk2[k] = 1.0/(1.0 + xn);
jRQ+2@n{E
gk2[KE-1-k] = 1.0/(1.0 + xn);
kiUk4&1
gk3[k] = (1.0 - xn)/(1.0 + xn);
K* vU5S
gk3[KE-1-k] = (1.0 - xn)/(1.0 + xn);
HW[L [&/
xxn = (double)(npml-k-.5)/(double)npml;
%oof}=MxCL
xn = .33*pow(xxn,3);
!Q %P%P<$
gk1[k] = xn;
(! KG)!
gk1[KE-2-k] = xn;
p3]Q^KFS
fk2[k] = 1.0/(1.0 + xn);
jRj=Awy
fk2[KE-2-k] = 1.0/(1.0 + xn);
}[!92WS/ee
fk3[k] = (1.0 - xn)/(1.0 + xn);
ls|LCQPx
fk3[KE-2-k] = (1.0 - xn)/(1.0 + xn);
pv# 2]v
}
I.V:q!4*
PHr a+NY#A
/* Calculate gax ,gbx */
:tj-gDa\Y
for ( i = ia; i<=ib; i++)
>k(MU mhX
{
WUoOGbA `
for ( j = ja; j<=jb; j++)
i[nF .I5*f
{
,YLF+^w-
for ( k =ka; k<=kb; k++)
T *>` ,}J
{
= r4!V>
xdist = (ic-i-.5);
7 y<1 LQ;}
ydist = (jc-j);
Uems\I0
zdist = (kc-k);
: 1 aL ?
dist = sqrt (pow(xdist,2) + pow(ydist,2) + pow(zdist,2));
@;^Y7po6u
if ( dist <= radius )
|:[9O`U)s
{
Mr3-q
gax
[j][k] = 1./(epsz_r + (sigma*dt/epsz));
#&Is GyU
gbx
[j][k] = sigma*dt/epsz;
? Rr2/W#F
}
~[W#/kd1n
}
@,OT/egF4:
}
9F##F-%x
}
SW 8x]B
/* Calculate gay ,gby*/
P3o @g kXP
for ( i = ia; i<=ib; i++)
4IsG=7
{
WeJ@x L
for ( j = ja; j<=jb; j++)
;hDa@3|]34
{
?6N3tk-2
for ( k =ka; k<=kb; k++)
hYg'2OG
{
!a- b6A a
xdist = (ic-i);
MDO$m g
ydist = (jc-j-.5);
M\2"gT-LV
zdist = (kc-k);
4na8
dist = sqrt (pow(xdist,2) + pow(ydist,2) + pow(zdist,2));
WlVl[/qt
if ( dist <= radius )
+t!S'|C
{
S2^>6/[xM
gay
[j][k] = 1./(epsz_r + (sigma*dt/epsz));
*3+-W
gby
[j][k] = sigma*dt/epsz;
8OFj0S1r`
}
eV!L^>>>
}
w# y2_
}
f ~h~5
}
q{' ~+Nq
/* Calculate gaz ,gbz */
o( G"k
for ( i = ia; i<=ib; i++)
-n))*.V
{
D$r Uid
for ( j = ja; j<=jb; j++)
?dq#e9
{
EbK0 j?
for ( k =ka; k<=kb; k++)
-j`LhS~|
{
3Qn!y\#
xdist = (ic-i);
VLvS$0(}Z
ydist = (jc-j);
gPXa>C
zdist = (kc-k-.5);
`!4,jd
dist = sqrt (pow(xdist,2) + pow(ydist,2) + pow(zdist,2));
F4C!CUI
if ( dist <= radius )
2}WDw>V
{
^twivNB
gaz
[j][k] = 1./(epsz_r +(sigma*dt/epsz));
+wfVL|.Wq
gbz
[j][k] = sigma*dt/epsz;
T3,1m=S
}
+ ,4" u
}
e@]-D FG
}
ff2d @P,!
}
_eQ- `?
//源参数
uh3) 0.nR
t0 = 40.0;
5J .0&Dda
spread = 10.0;
h #.N3o
*I*i> ==Z
clock_t start, finish;
Paf% rv2
start = clock();
:nJgwp()@
/* outFile.open("F:\\ez[40][25][25].txt");*/
{7"0,2 Hb?
for ( T = 1; T < 1001; T++ )
x_<qzlQt
{
NC"yDWnO'
/* Calculate the incident buffer */
CboLH0Fa
for ( j=1; j<JE; j++)
i! H!;z#
{
5LW}h^N
ez_inc[j] = ez_inc[j] + .5*(hx_inc[j-1] - hx_inc[j]);
L'A)6^d@S
}
\=e8%.#@J
//Fourier Transform of the incident field
p9[6^rj x8
ine[0] = ine[0] + dt*cos(2*pi*f*dt*T)*ez_inc[ja-1];
5>_5]t {
ine[1] = ine[1] - dt*sin(2*pi*f*dt*T)*ez_inc[ja-1];
L36Yx7gT<
/* Source */
^ 4>k%d
//pulse = exp (-.5*(pow(( t0-T )/spread,2)));
4 [2_,9}
pulse = exp (2*pi*f*T*dt);
yi 6 N-7
ez_inc[3] = pulse;
O<*5$,K9
/* Boundary conditions for the incident buffer */
h0 |}TV^UJ
ez_inc[0] = ez_low_m2;
Ao .\
ez_low_m2 = ez_low_m1;
FY'dJY3O
ez_low_m1 = ez_inc[1];
`N8 7 h"
ez_inc[JE-1] = ez_high_m2;
B(5c9DI`
ez_high_m2 = ez_high_m1;
VG@};dwbz*
ez_high_m1 = ez_inc[JE-2];
x?,9_va]
/* Calculate the Dx field */
DtxE@,
for ( i=1; i<ia; i++)
O1y|v[-BW
{
\=2m7v#E
for ( j=1; j<JE; j++)
1%~ZRmd e
{
onei4c>@
for (k=1; k<KE; k++)
c4mh EE-
{
9U_ks[Qa
curl_h = ( hz
[j][k] - hz
[j-1][k] - hy
[j][k] + hy
[j][k-1] );
_^S]g mE
idxl
[j][k] = idxl
[j][k] + curl_h;
zc+@lJy
dx
[j][k] = gj3[j]*gk3[k]*dx
[j][k] + gj2[j]*gk2[k]*.5*(curl_h + gi1
*idxl
[j][k]);
iz^qR={bW
}
Qu\@Y[eia5
}
X&\ d)/Y
}
rNN j0zw>
for ( i=ia; i<=ib; i++)
l|`^*%W@u6
{
x~ DLW 1I
for ( j=1;j<JE; j++)
QfpuZEUK
{
=?Fkn4t
for ( k=1; k<KE; k++)
lU4}B`#"v
{
` }gbc 69
curl_h = ( hz
[j][k] - hz
[j-1][k] - hy
[j][k] + hy
[j][k-1] );
6z,Dyy]tl
dx
[j][k] = gj3[j]*gk3[k]*dx
[j][k] + gj2[j]*gk2[k]*.5*curl_h;
1T|$BK@)
}
GHc/Zc"iX
}
?A*Kg;IU
}
d j5 hv~
for ( i=ib+1; i<IE; i++)
hsNWqk qys
{
Pcd *">v
ixh = i-ib-1;
q0xjA
for ( j=1; j<JE; j++)
^rAa"p 9
{
I*|P@ 0
for ( k=1; k<KE; k++)
&l2TeC@;
{
WC q /c6 D
curl_h = ( hz
[j][k] - hz
[j-1][k] - hy
[j][k] + hy
[j][k-1] );
+ %*&.@z_
idxh[ixh][j][k] = idxh[ixh][j][k] + curl_h;
:z"!k zdJ
dx
[j][k] = gj3[j]*gk3[k]*dx
[j][k] + gj2[j]*gk2[k]*.5*(curl_h + gi1
*idxh[ixh][j][k]);
<, @%*G1-
}
mLbN/M
}
cu 9Qwm
}
3z =^(Y
/* Calculate the Dy field */
7L(e h7
for ( i=1; i<IE; i++)
vDj;>VE2b
{
!X^Ce)1K
for ( j=1; j<ja; j++)
Sb&lhgW]c
{
J >0b1
for (k=1; k<KE; k++)
s95F#>dr
{
,XK Cz ]8V
curl_h = ( hx
[j][k] - hx
[j][k-1] - hz
[j][k] + hz[i-1][j][k] );
1} m3 ;
idyl
[j][k] = idyl
[j][k] + curl_h;
oy.[+EI`|
dy
[j][k] = gi3
*gk3[k]*dy
[j][k]+gi2
*gk2[k]*.5*(curl_h+ gj1[j]*idyl
[j][k]);
2b@tj 5
}
s|yVAt|=
}
b 'p4wE>
}
8 ;g Xg
for ( i=1; i<IE; i++)
m\u26 ` M
{
47B y`Jh71
for ( j=ja;j<=jb; j++)
o+ O}Te
{
@E9" Zv-$
for ( k=1; k<KE; k++)
Yc Q=vt{
{
Tbbz'b;{
curl_h = ( hx
[j][k] - hx
[j][k-1] - hz
[j][k] + hz[i-1][j][k] );
B|=|.qp$)
dy
[j][k] = gi3
*gk3[k]*dy
[j][k]+gi2
*gk2[k]*.5*curl_h;
K >tf,
}
(U1]:tZ<.
}
wFS2P+e;X
}
-6~' cm
for ( i=1; i<IE; i++)
A$/\1282
{
q3e8# R)l
for ( j=jb+1; j<JE; j++)
$hrIO+
{
p^)B0[P9
jyh = j-jb-1;
k-LEI}h
for ( k=1; k<KE; k++)
lHfe<j]
{
f'En#-?O
curl_h = ( hx
[j][k] - hx
[j][k-1] - hz
[j][k] + hz[i-1][j][k] );
C}}/)BYi
idyh
[jyh][k] = idyh
[jyh][k] + curl_h;
Y g|lq9gD
dy
[j][k] = gi3
*gk3[k]*dy
[j][k]+gi2
*gk2[k]*.5*(curl_h+ gj1[j]*idyh
[jyh][k]);
vv0Q$ O->
}
Bik*b)9y2
}
2B1xUj ]
}
GC@+V|u
/* Incident Dy */
VU.@R,
for ( i=ia; i<=ib; i++)
@J 'YV{]
{
r4}*l 7Q
for ( j=ja; j<=jb-1; j++)
0?Q_@Y
{
u0s8yPA
dy
[j][ka] = dy
[j][ka] - .5*hx_inc[j];
qi[Z,&
dy
[j][kb+1] = dy
[j][kb+1] + .5*hx_inc[j];
5N7H{vT_
}
^ -)txC5{T
}
AC%JC+
/* Calculate the Dz field */
k;y w#Af8
for ( i=1; i<IE; i++)
pV[ ''
{
n`7f"'/:
for ( j=1; j<JE; j++)
AG6K da J
{
O'*@ Ytn
for (k=1; k<ka; k++)
CON0E~"
{
W m&
curl_h = ( hy
[j][k] - hy[i-1][j][k] - hx
[j][k] + hx
[j-1][k] );
NaU r!s
idzl
[j][k] = idzl
[j][k] + curl_h;
3)Ac"nuyqH
dz
[j][k] = gi3
*gj3[j]*dz
[j][k]+gi2
*gj2[j]*.5*(curl_h+ gk1[k]*idzl
[j][k]);
[]kN16F
}
Of}|ib^t
}
d=* x#In
}
N5U)*U'-u
for ( i=1; i<IE; i++)
dR^"X3$
{
-~f51 1<
for ( j=1;j<JE; j++)
j+4H}XyE
{
.3qu9eP
for ( k=ka; k<=kb; k++)
cW8 \d
{
BHFY%6J!
curl_h = ( hy
[j][k] - hy[i-1][j][k] - hx
[j][k] + hx
[j-1][k] );
`_X; .U.Mv
dz
[j][k] = gi3
*gj3[j]*dz
[j][k]+gi2
*gj2[j]*.5*curl_h;
yl&UM qI(
}
TX8<J>x
}
\Gp*x\<^Z
}
+/{L#e>
for ( i=1; i<IE; i++)
hcCp,b
{
RD,5AShP
for ( j=1; j<JE; j++)
"B7`'jz
{
Qeq5 gN]
for ( k=kb+1; k<KE; k++)
+_l^ #?o,
{
n'FwM\
kzh = k-kb-1;
;QCrHqRT`
curl_h = ( hy
[j][k] - hy[i-1][j][k] - hx
[j][k] + hx
[j-1][k] );
z;\d L
idzh
[j][kzh] = idzh
[j][kzh] + curl_h;
nR;D#"p%
dz
[j][k] = gi3
*gj3[j]*dz
[j][k]+gi2
*gj2[j]*.5*(curl_h+ gk1[k]*idzh
[j][kzh]);
DPn=n9n2
}
C#pZw [
}
cu!W4Ub<
}
RI n9(r
/* Incident Dz */
_Xcn N:Rt
for ( i=ia; i<=ib; i++)
Y k7-`
{
N 5zlT
for ( k=ka; k<=kb; k++)
4-GXmC
{
GwU?wII j^
dz
[ja][k] = dz
[ja][k] + .5*hx_inc[ja-1];
'])2k@o@
dz
[jb][k] = dz
[jb][k] - .5*hx_inc[jb];
arK_oh0B
}
lK7m=[ j
}
2.e vx
/* Calculate the E from D field */
266oTER]v:
for ( i=1; i<IE-1; i++)
& [mZD,
{
0IkM
for ( j=1; j<JE-1; j++)
o|pT;1a"
{
%qf V+^
for ( k=1; k<KE-1; k++)
AV8T
{
~p0c3*
ex
[j][k] = gax
[j][k]*(dx
[j][k]-ix
[j][k]);
0U/,aHvhP
ix
[j][k] = ix
[j][k]+gbx
[j][k]*ex
[j][k];
bvBHYf:^
ey
[j][k] = gay
[j][k]*(dy
[j][k]-iy
[j][k]);
sM[I4 .A3
iy
[j][k] = iy
[j][k]+gby
[j][k]*ey
[j][k];
YGvUwj'2a
ez
[j][k] = gaz
[j][k]*(dz
[j][k]-iz
[j][k]);
{svn=H /
iz
[j][k] = iz
[j][k]+gbz
[j][k]*ez
[j][k];
58xnB !h\}
}
$ <8~k^
}
n.oUVr=nX
}
F^$led1/F
/* Remember: part of the PML is E = 0 at the edges */
UO Ug 4
for (j = 0;j < JE-1;j++)
K5t0L!6<+
{
yF._*9Q3hK
for (k = 0;k < KE-1;k++)
IeX^4 rc(
{
`Mj}md;O"
ez[0][j][k] = 0.0;
bizTd
ez[IE-1][j][k] = 0.0;
/t<@"BoV
}
?L) !pP]
}
i+3fhV
for (i = 0;i < IE-1;i++)
1nknSw#
{
G2{O 9
for (k = 0;k < KE;k++)
:G w~7v_
{
LjFqZrH
ez
[0][k] = 0.0;
d5 Edu44
ez
[JE-1][k] = 0.0;
`Hx~UH)
}
5+Mdh`
}
7xo4-fIuT
for (i = 0;i < IE-1;i++)
zLw{ {|
{
7}O.wUKw%
for (j = 0;j < JE-1;j++)
[Qt?W gPj
{
|! SO G
ez
[j][0] = 0.0;
YV4#%I!<
ez
[j][KE-1] = 0.0;
~ I}9;XT
}
Q?t^ @
}
3oZ=k]\
'|>9C^E9X
/*outFile << setw(22) << setprecision(8) << T << setw(22) << setprecision(8) << ez[40][25][25]<< endl;*/
2Se?J)MN
[k$*4 u >
/* Calculate the incident field */
<+#o BN
for ( j=0; j<JE-1; j++)
9i5,2~
{
b#)U UGmI
hx_inc[j] = hx_inc[j] + .5*(ez_inc[j]-ez_inc[j+1]);
J<gJc *Q
}
Q:C$&-$
/* Calculate the Hx field */
Lw7=+h)
for ( i=0; i<ia; i++)
o\otgyoh
{
gy: %l
for ( j=0; j<JE-1; j++)
'' OfS D_g
{
qB]i6*
for (k=0; k<KE-1; k++)
s pLZ2]A
{
FbnO/! $8
curl_e = ( ey
[j][k+1] - ey
[j][k] - ez
[j+1][k] + ez
[j][k] );
|%a4` w
ihxl
[j][k] = ihxl
[j][k] + curl_e;
DHC+C4
hx
[j][k] = fj3[j]*fk3[k]*hx
[j][k]+fj2[j]*fk2[k]*.5*(curl_e+ fi1
*ihxl
[j][k]);
RPnRVJ&"Z
}
d'6|: z9c
}
DciwQcG
}
ysIh[1E~%:
for ( i=ia; i<=ib; i++)
= VLS/\ A
{
vR6B n
for ( j=0;j<JE-1; j++)
hJz):d>Im
{
Y<"7x#AB!
for ( k=0; k<KE-1; k++)
\'CN
{
iTh:N2/-vc
curl_e = ( ey
[j][k+1] - ey
[j][k] - ez
[j+1][k] + ez
[j][k] );
J/!cGr( B~
hx
[j][k] = fj3[j]*fk3[k]*hx
[j][k]+fj2[j]*fk2[k]*.5*curl_e;
^\:8w0Y^
}
^I mP`*X
}
22P$ ~ch
}
wu9=N ^x
for ( i=ib+1; i<IE; i++)
jDM w2#<
{
f&, {XZ
ixh = i-ib-1;
$2D uB
for ( j=0; j<JE-1; j++)
;!< Znw
{
lOwS&4UT
for ( k=0; k<KE-1; k++)
}}<Z,/O
{
Fk;o E'"D
curl_e = ( ey
[j][k+1] - ey
[j][k] - ez
[j+1][k] + ez
[j][k] );
i|)Su4Dw
ihxh[ixh][j][k] = ihxh[ixh][j][k] + curl_e;
/%aiEhL
hx
[j][k] = fj3[j]*fk3[k]*hx
[j][k]+fj2[j]*fk2[k]*.5*(curl_e+ fi1
*ihxh[ixh][j][k]);
Si 9Z>MR
}
+y2[msBs
}
}{ 9&:!uA
}
X;F?:Iw \
/* Incident Hx */
1=7jz ] t
for ( i=ia; i<=ib; i++)
g42)7
{
%Pqk63QF
for ( k=ka; k<=kb; k++)
oY~ D g
{
$eV$2p3H
hx
[ja-1][k] = hx
[ja-1][k] + .5*ez_inc[ja];
l=N2lHU
hx
[jb][k] = hx
[jb][k] - .5*ez_inc[jb];
4n1; Bh$
}
|]Z:&[D]i
}
Rv,JU6>i
IPSF]"}~
/* Calculate the Hy field */
w=h1pwY
for ( i=0; i<IE-1; i++)
ajRSMcKb7i
{
=$8nUX`
for ( j=0; j<ja; j++)
mV73 \P6K
{
8mLP5s!7
for (k=0; k<KE-1; k++)
{K{EOB_ u
{
9z;HsU v
curl_e = ( ez[i+1][j][k] - ez
[j][k] - ex
[j][k+1] + ex
[j][k] );
e^YHJ>@
ihyl
[j][k] = ihyl
[j][k] + curl_e;
SQ,?N XZ
hy
[j][k] = fi3
*fk3[k]*hy
[j][k]+fi2
*fk2[k]*.5*(curl_e+ fj1[j]*ihyl
[j][k]);
K ]OK:hY4
}
C9DJO:f.2y
}
t%zpNd2lk
}
wfvU0]wk}
for ( i=0; i<IE-1; i++)
$uui:wU%Q
{
@;xMs8@
for ( j=ja;j<=jb; j++)
K-<^ $VWh
{
)Dp/('Z2
for ( k=0; k<KE-1; k++)
+(D$9{y
{
LW sP ya
curl_e = ( ez[i+1][j][k] - ez
[j][k] - ex
[j][k+1] + ex
[j][k] );
H xV#WoYKj
hy
[j][k] = fi3
*fk3[k]*hy
[j][k]+fi2
*fk2[k]*.5*curl_e;
B<8N96fx
}
j6~nE'sQ
}
zOEY6lAwI
}
Q(d9n 8
for ( i=0; i<IE-1; i++)
&bz:K8c
{
*{undZ?(>
for ( j=jb+1; j<JE; j++)
{HC@u{K -
{
H+}"q$
jyh = j-jb-1;
49Df?sx
for ( k=0; k<KE-1; k++)
Uq_j\A;c
{
bb42v7?
curl_e = ( ez[i+1][j][k] - ez
[j][k] - ex
[j][k+1] + ex
[j][k] );
;VO.!5W@eg
ihyh
[jyh][k] = ihyh
[jyh][k] + curl_e;
m]_FQWfet
hy
[j][k] = fi3
*fk3[k]*hy
[j][k]+fi2
*fk2[k]*.5*(curl_e+ fj1[j]*ihyh
[jyh][k]);
C.^Ven
}
-"Y{$/B
}
D9mz9
}
HBys
/* Incident Hy */
I]Tsz'T!9
for ( j=ja; j<=jb; j++)
/<CSVJ_r
{
N!Qg; (
for ( k=ka; k<=kb; k++)
+#b:d=v!
{
cWGDee(
hy[ia-1][j][k] = hy[ia-1][j][k] - .5*ez_inc[j];
._ wkj
hy[ib][j][k] = hy[ib][j][k] + .5*ez_inc[j];
xCzebG["
}
=&0wr6
}
5F t5@UF~
/* Calculate the Hz field */
cr?7O;,
for ( i=0; i<IE-1; i++)
5 G0 $
{
LN l#h
for ( j=0; j<JE-1; j++)
OX%MP!#KU
{
TvNY:m6.%
for (k=0; k<ka; k++)
rV({4cIe9R
{
UY2X
curl_e = ( ex
[j+1][k] - ex
[j][k] - ey[i+1][j][k] + ey
[j][k] );
rPc7(,o*
ihzl
[j][k] = ihzl
[j][k] + curl_e;
L{l6Dd43q
hz
[j][k] = fi3
*fj3[j]*hz
[j][k]+fi2
*fj2[j]*.5*(curl_e+ fk1[k]*ihzl
[j][k]);
@cRZk`|1n
}
. *$OQA
}
=z {JgD/
}
; UiwH
for ( i=0; i<IE-1; i++)
MRr</o
{
p5O",3,A4
for ( j=0;j<JE-1; j++)
S.B<pj gt
{
2,vB'C AI
for ( k=ka; k<=kb; k++)
1iL 'V-y
{
bi ozZ
curl_e = ( ex
[j+1][k] - ex
[j][k] - ey[i+1][j][k] + ey
[j][k] );
]J9cVp
hz
[j][k] = fi3
*fj3[j]*hz
[j][k]+fi2
*fj2[j]*.5*curl_e;
4`N t{
}
zL5d0_E9
}
FLoNE>q
}
&bgvy'p
for ( i=0; i<IE-1; i++)
!|2VWI}
{
v{i 7h|e
for ( j=0; j<JE-1; j++)
;o/>JHGj
{
5!C_X5M
for ( k=kb+1; k<KE; k++)
2M)]!lYy
{
x 5dWBGH
kzh = k-kb-1;
obK*rdg ,
curl_e = ( ex
[j+1][k] - ex
[j][k] - ey[i+1][j][k] + ey
[j][k] );
*'"T$ib
ihzh
[j][kzh] = ihzh
[j][kzh] + curl_e;
Nf3.\eR
hz
[j][k] = fi3
*fj3[j]*hz
[j][k]+fi2
*fj2[j]*.5*(curl_e+ fk1[k]*ihzh
[j][kzh]);
