登 录
註 冊
论坛
微波仿真网
注册
登录论坛可查看更多信息
微波仿真论坛
>
时域有限差分法 FDTD
>
大家帮忙看看这个3D_FDTD介质球的RCS有什么问 ..
发帖
回复
1
2
3
4119
阅读
29
回复
大家帮忙看看这个3D_FDTD介质球的RCS有什么问题
离线
cc81372365
UID ：27752
注册：
2009-03-18
登录：
2023-11-09
发帖：
345
等级：
论坛版主
0楼
发表于: 2010-11-29 13:16:05
//****************************************************************//
fNNik7
//
X0+M|8:
//********* 3D_介质球散射（RCS）FDTD程序 **********//
kf >L
//
n&;-rj^qq
//****************************************************************//
^_o9%)RL(
//
;&?l1Vu
//入射波为Y方向的平面波（theta=pi/2 , phi=pi/2), RCS观察角为（theta=pi/2 , phi = 0到2pi）
a]nyZdt`
// 数字0，1代表场值经过傅里叶变换后的复数的实部和虚部
VSQxlAGk@
// 金属方柱置于中心位置（ic,jc）,激励源为高斯脉冲源
n=tg{_9f%
#include <iostream>
QJKVNOo
#include <iomanip>
gW_^GrK pI
#include <fstream>
]CZ&JL
//#include <omp.h>
L; f
#include <time.h>
&bw ``e&c
#include <valarray>
~ @Au <
#include <math.h>
69K*]s
#include <stdio.h>
.>bvI1
#include <stdlib.h>
Ii4lwZnz
3|g]2|~w@h
#include "complex_operation.h"
mbCY\vEl
#define IE 60
@o6^"
#define JE 60
+pMjm&CF
#define KE 60
~n $e
#define ia 15
Bvy(vc=UDW
#define ja 15
j_Z "=
#define ka 15
b*@y/ e\u`
#define length 41
Ow7}&\;^-
using namespace std;
)Og,VXEB
ofstream outFile;
Y/aNrIK7
int main()
F MBzTD
{
+w'{I`QIL0
int i,j,k,ic,jc,kc,ib,jb,kb;
l: X]$2;
int npml,n_pml,T;
d}e/f)(
double radius,v,lamda,Z,f,kk;
jf'#2-
double epsz_r,epsz,muz,muz_r,ddx,dt,pi,sigma;
*q/oS8vavd
double xn,xxn;
@vc vte
double t0,spread,pulse;
3[E) /~-
double ez_inc[JE]={0},hx_inc[JE]={0};
<D!\"C
double ez_low_m1=0,ez_low_m2=0,ez_high_m1=0,ez_high_m2=0;
b6E,u*)"
int ixh,jyh,kzh;
%}%vey
int co=0;
;}=[( eqA
double rcs[361]={0.0};
rBUdHd9
double ine[2] = {0.0};
sG~5O\,E
double ine_fudu = 0.0;
}\oy?_8~
double curl_h,curl_e;
1]l m0bfs
double dist,xdist,ydist,zdist;
Tfba3+V
double dx[IE][JE][KE]={0},dy[IE][JE][KE]={0},dz[IE][JE][KE]={0};
&v# *
double ex[IE][JE][KE]={0},ey[IE][JE][KE]={0},ez[IE][JE][KE]={0};
DMY?'Nts!
double hx[IE][JE][KE]={0},hy[IE][JE][KE]={0},hz[IE][JE][KE]={0};
d;kdw
double ix[IE][JE][KE]={0},iy[IE][JE][KE]={0},iz[IE][JE][KE]={0};
D{!N Tr
double gax[IE][JE][KE]={0},gay[IE][JE][KE]={0},gaz[IE][JE][KE]={0};
WV'FW)%
double gbx[IE][JE][KE]={0},gby[IE][JE][KE]={0},gbz[IE][JE][KE]={0};
ou[_ y
double idxl[ia][JE][KE]={0},idxh[ia][JE][KE]={0};
Zg@N MT
double idyl[IE][ja][KE]={0},idyh[IE][ja][KE]={0};
T; u >]"S
double idzl[IE][JE][ka]={0},idzh[IE][JE][ka]={0};
k!lz_Y
double ihxl[ia][JE][KE]={0},ihxh[ia][JE][KE]={0};
fy`e)?46
double ihyl[IE][ja][KE]={0},ihyh[IE][ja][KE]={0};
f/:XIG
double ihzl[IE][JE][ka]={0},ihzh[IE][JE][ka]={0};
f9Hm2wV
double gi1[IE]={0},gi2[IE]={0},gi3[IE]={0};
XdDy0e4{%<
double gj1[JE]={0},gj2[JE]={0},gj3[JE]={0};
!a[1rQH
double gk1[KE]={0},gk2[KE]={0},gk3[KE]={0};
_= .f+1W
double fi1[IE]={0},fi2[IE]={0},fi3[IE]={0};
2@ <x%T
double fj1[JE]={0},fj2[JE]={0},fj3[JE]={0};
k 8 ;
double fk1[KE]={0},fk2[KE]={0},fk3[KE]={0};
u}7#3JfLn
double Ex_down[length][length][2]={0},Ey_down[length][length][2]={0};
5r:SBt|/
double Hx_down[length][length][2]={0},Hy_down[length][length][2]={0};
aiKZ$KLC
double Ex_up[length][length][2]={0},Ey_up[length][length][2]={0};
n>Rt9
double Hx_up[length][length][2]={0},Hy_up[length][length][2]={0};
6J|Y+Y$
double Ex_left[length][length][2]={0},Ez_left[length][length][2]={0};
#.(6.Li
double Hx_left[length][length][2]={0},Hz_left[length][length][2]={0};
xM/B"SG2
double Ex_right[length][length][2]={0},Ez_right[length][length][2]={0};
ZMHb
double Hx_right[length][length][2]={0},Hz_right[length][length][2]={0};
TuBl9 p'6
double Ey_front[length][length][2]={0},Ez_front[length][length][2]={0};
.cQ<F4)!tu
double Hy_front[length][length][2]={0},Hz_front[length][length][2]={0};
df Dz/sD*
double Ey_back[length][length][2]={0},Ez_back[length][length][2]={0};
}lk_Oe1
double Hy_back[length][length][2]={0},Hz_back[length][length][2]={0};
^/3R/;?
ic = IE/2;
Z5 uetS^
jc = JE/2;
Q`Q%;%t
kc = KE/2;
'y [eH
ib = IE - ia -1;
xsx @aF
jb = JE - ja -1;
As n7 ;x0;
kb = KE - ka -1;
\V|\u= @H
//外推边界位置
%s;#epP$
int dis = 5;
pN) 9 G O5
const int i_front = ib + dis;
TvE M{
const int i_back = ia - dis;
T #&9|
const int j_right = jb + dis;
iR9 $E
const int j_left = ja - dis;
ag-\(i;K]
const int k_up = kb + dis;
LsnM5GU7
const int k_down = ka - dis;
0@yHT-Dy
pi = 3.1415926;
[lqwzW{(UN
epsz = 8.8e-12; // 真空介电常数
4ey m$UWw
epsz_r = 4.0; // 相对介电常数
sS, Swgr
muz = 4*pi*1e-7; // 真空磁导率
Y. ]FVq
muz_r = 1.0; // 相对磁导率
V f JYYR
sigma = 0; // 介质电导率
ngprTMO$&
v = 3.0e8; // 波速
(bxSN@hp2
f = 9.375e9; // 频率 10GHz
M1^?_; B
lamda = v/f; // 波长 0.03m
~n$VCLa
kk = 2*pi/lamda; // 波数
7 u!i)<pn
Z = sqrt(muz/epsz); // 真空波阻抗 377
TK pka]nJ
ddx = lamda/42; // 网格尺寸
E rK5iTSD
dt = ddx/(2.0*v); // 时间步长
bsw0+UY=9
/* Initialize the arrays */
_uuxTNN0x*
for ( i=0; i<IE; i++)
l+'@y (}Q
{
1D 03Nbh|5
for ( j=0; j<JE; j++)
QCFLi n+r
{
[];*9vxW
for ( k=0; k<KE; k++)
`H6-g=C
{
G[1: <Vg8
ex
[j][k] = 0.0;
epsh&)5a*
ey
[j][k] = 0.0;
s l|n]#)
ez
[j][k] = 0.0;
Y]aVa2!Wb
dx
[j][k] = 0.0;
37DyDzW)'
dy
[j][k] = 0.0;
,y}?Z 8?63
dz
[j][k] = 0.0;
B] dvX
hx
[j][k] = 0.0;
=B g
hy
[j][k] = 0.0;
hA.? 19<Z
hz
[j][k] = 0.0;
}>I|\Z0I
ix
[j][k] = 0.0;
!rzbm& @
iy
[j][k] = 0.0;
aA?Uf~ "t
iz
[j][k] = 0.0;
FUic7>
gax
[j][k] = 1.0;
ufE;rcYE
gay
[j][k] = 1.0;
.5*h']iFr1
gaz
[j][k] = 1.0;
*HFRG)[V
gbx
[j][k] = 0.0;
`9BZ))Pg
gby
[j][k] = 0.0;
-&JUg o=
gbz
[j][k] = 0.0;
tLXwszR0 r
}
'M>QA"*48E
}
.}n%gc~ A
}
W 86`R
for ( i=0; i<ia; i++)
4. :2!Q
{
_ KBN
for ( j=0; j<JE; j++)
?.d6!v A
{
0|kkwZVPn
for ( k=0; k<KE; k++)
'thWo wE
{
1zwk0={x-%
idxl
[j][k] = 0.0;
Q>R jv.1
idxh
[j][k] = 0.0;
h+)XL s
ihxl
[j][k] = 0.0;
u`I&&
ihxh
[j][k] = 0.0;
f8yE>qJP
}
RKoM49W
}
b^[Ab:`}[V
}
e&WlJ
for ( i=0; i<IE; i++)
m2 0:{fld
{
@X@?jj&
for ( j=0; j<ja; j++)
1 p&e:v
{
p4QQ5O$;
for ( k=0; k<KE; k++)
-j1?l Y
{
\$GM4:R D
idyl
[j][k] = 0.0;
E9e|+$
idyh
[j][k] = 0.0;
N>kY$ *
ihyl
[j][k] = 0.0;
b&[bfM<
ihyh
[j][k] = 0.0;
a *?b nw?
}
Fk(nf9M%
}
:.8@ xVH
}
2P3,\L
for ( i=0; i<IE; i++)
d:C|laZHn
{
-D*,*L
for ( j=0; j<JE; j++)
{nvF >
{
|>_e& }Y%L
for ( k=0; k<ka; k++)
O# n<`;W
{
Vwxb6,}Z
idzl
[j][k] = 0.0;
A [W3.$s
idzh
[j][k] = 0.0;
^3re*u4b=
ihzl
[j][k] = 0.0;
zh8\ _> +
ihzh
[j][k] = 0.0;
Cngi5._Lb
}
qiEw[3Za]'
}
!1G6ZC:z
}
jSi\ /( E
//PML初始化
t&5N{C:
for ( i=0; i<IE; i++)
{A<pb{<u
{
{}gx;v)
gi1
= 0;
% gBulvg
fi1
= 0;
2F5*C
gi2
= 1;
lICpfcc(+
fi2
= 1;
90g=&O5@O
gi3
= 1;
>\f'Q Q
fi3
= 1;
okwkMd-yW
}
{`RCh]W
for ( j=0; j<JE; j++)
KHnq%#
{
t`F<lOKj
gj1[j] = 0;
6:\0=k5
fj1[j] = 0;
CuT~ Bj
gj2[j] = 1;
iO$Z?Dyg9
fj2[j] = 1;
;og[ q
gj3[j] = 1;
eKpWFP 0
fj3[j] = 1;
?7;_3+T#
}
HNXMM
for ( k=0; k<KE; k++)
'xK ,|U
{
''p7!V?
gk1[k] = 0;
`^d [$IbDW
fk1[k] = 0;
Y \7WCaSgi
gk2[k] = 1;
g>gVO@"b2
fk2[k] = 1;
z5njblUz
gk3[k] = 1;
oItEGJ|
fk3[k] = 1;
"W71 #n+ [
}
le7!:4 /8
radius = 10 ;
+"yt/9AO
n_pml = 7; //PML层厚度
?aMd#.&
npml = 7;
8PRKS J[@K
for (i = 0; i < n_pml; i++)
fx74h{3u
{
}Bk>'
xxn = (double)(npml-i)/(double)npml;
0:C ^-zrx
xn = 0.33*pow(xxn,3);
s?Wkh `b
fi1
= xn;
gdj,e ^
fi1[IE-1-i] = xn;
itU P%
gi2
= 1.0/(1.0 + xn);
XEX-NE"]
gi2[IE-1-i] = 1.0/(1.0 + xn);
`4Db( ~
gi3
= (1.0 - xn)/(1.0 + xn);
xZE% Gf_U
gi3[IE-1-i] = (1.0 - xn)/(1.0 + xn);
=|j~ *6Hd
xxn = (double)(npml-i-.5)/(double)npml;
zn+ 5pn&?
xn = .33*pow(xxn,3);
zWA~0l.2
gi1
= xn;
PI-o)U$Ehv
gi1[IE-2-i] = xn;
0?80V'
fi2
= 1.0/(1.0 + xn);
1yK= Yf%B
fi2[IE-2-i] = 1.0/(1.0 + xn);
9coN >y
fi3
= (1.0 - xn)/(1.0 + xn);
$ca> b X]
fi3[IE-2-i] = (1.0 - xn)/(1.0 + xn);
jhx @6[
}
"e ;wN3/bF
for (j = 0; j < n_pml; j++)
WHk rd8
{
iJaA&z5sr
xxn = (double)(npml-j)/(double)npml;
M[:},?ah0
xn = .33*pow(xxn,3);
_N1UL?
fj1[j] = xn;
ZHN}:W/p
fj1[JE-1-j] = xn;
CqAv^ n7 }
gj2[j] = 1.0/(1.0 + xn);
o0 &pSCK
gj2[JE-1-j] = 1.0/(1.0 + xn);
{G i:W/jJ
gj3[j] = (1.0 - xn)/(1.0 + xn);
8GKqPS+
gj3[JE-1-j] = (1.0 - xn)/(1.0 + xn);
.dw;b~p
xxn = (double)(npml-j-.5)/(double)npml;
JI7 .:k;
xn = .33*pow(xxn,3);
g[D `.
gj1[j] = xn;
X/AA8QV o
gj1[JE-2-j] = xn;
oMV^W^<
fj2[j] = 1.0/(1.0 + xn);
e b*w$|y6"
fj2[JE-2-j] = 1.0/(1.0 + xn);
<&m `)FJ
fj3[j] = (1.0 - xn)/(1.0 + xn);
u.!<)VIJx
fj3[JE-2-j] = (1.0 - xn)/(1.0 + xn);
;p8,= w
}
U#gv ~)\k
0(h'ZV
for (k = 0; k < n_pml; k++)
-e gu5#d>
{
r>kDRIHB
xxn = (double)(npml-k)/(double)npml;
uzO3 _.4Y
xn = .33*pow(xxn,3);
/9k}Ip
fk1[k] = xn;
mA" 82"
fk1[KE-1-k] = xn;
*nYb9.T]i
gk2[k] = 1.0/(1.0 + xn);
'kE^oX_
gk2[KE-1-k] = 1.0/(1.0 + xn);
F7N4qq1
gk3[k] = (1.0 - xn)/(1.0 + xn);
*{%d{x}l
gk3[KE-1-k] = (1.0 - xn)/(1.0 + xn);
j"J[dlm2M
xxn = (double)(npml-k-.5)/(double)npml;
LQ'VhNU
xn = .33*pow(xxn,3);
nep-?7 x
gk1[k] = xn;
Fq` wx
gk1[KE-2-k] = xn;
zK f.jpF^
fk2[k] = 1.0/(1.0 + xn);
\?K>~{)
fk2[KE-2-k] = 1.0/(1.0 + xn);
mEVne.D
fk3[k] = (1.0 - xn)/(1.0 + xn);
-)o0P\cTEt
fk3[KE-2-k] = (1.0 - xn)/(1.0 + xn);
^hIKDc!.m
}
>XomjU[srQ
0o*
/* Calculate gax ,gbx */
O~u@J'4
for ( i = ia; i<=ib; i++)
Ng;F hv+
{
j|c6BdROl
for ( j = ja; j<=jb; j++)
FZJyqqA$_
{
:978D0}{p
for ( k =ka; k<=kb; k++)
P%^\<#Ya7
{
w @Gk#
xdist = (ic-i-.5);
Ag@R 60#
ydist = (jc-j);
n q3B(
zdist = (kc-k);
3^%sz!jK+
dist = sqrt (pow(xdist,2) + pow(ydist,2) + pow(zdist,2));
otSF8[
if ( dist <= radius )
0ofl,mXW
{
2Hj;o
gax
[j][k] = 1./(epsz_r + (sigma*dt/epsz));
o ImW
gbx
[j][k] = sigma*dt/epsz;
u4,b%h.
}
N \Wd 0b
}
}RPeA cbU_
}
K" U!SWv
}
n?z^"vv$i
/* Calculate gay ,gby*/
% m0x ]
for ( i = ia; i<=ib; i++)
= S|^pN
{
[ Cu3D
for ( j = ja; j<=jb; j++)
|&xjuBC
{
{I~[a #^
for ( k =ka; k<=kb; k++)
EHC^ [5
{
J`@# yHL
xdist = (ic-i);
Xy'qgK?
ydist = (jc-j-.5);
$v*0 \O
zdist = (kc-k);
YkqauyV^
dist = sqrt (pow(xdist,2) + pow(ydist,2) + pow(zdist,2));
i<]Y0 _?s
if ( dist <= radius )
JJn+H&[B
{
M_monj}Z
gay
[j][k] = 1./(epsz_r + (sigma*dt/epsz));
J&jNONu?
gby
[j][k] = sigma*dt/epsz;
!YJ^BI
}
G*$a81dAX
}
!&=%#i
}
0Fi&7%
}
ok+-#~VTn
/* Calculate gaz ,gbz */
c(y~,hN&p
for ( i = ia; i<=ib; i++)
X/!37
{
.1 =8c\%
for ( j = ja; j<=jb; j++)
g :Z, ab4
{
NA`EG, 2
for ( k =ka; k<=kb; k++)
d PfD Pb
{
gc6T`O-_;
xdist = (ic-i);
t <Z)D0.
ydist = (jc-j);
w}jH,Ew
zdist = (kc-k-.5);
)x M P
dist = sqrt (pow(xdist,2) + pow(ydist,2) + pow(zdist,2));
~ jqh&u$(
if ( dist <= radius )
D;l)&"| r?
{
l 8_TeO
gaz
[j][k] = 1./(epsz_r +(sigma*dt/epsz));
]CI ZF,
gbz
[j][k] = sigma*dt/epsz;
V^^nJs tV
}
ErDt ~FH
}
$ *u{i4b
}
U2(|/M +
}
[! v| M
//源参数
G?OwhX
t0 = 40.0;
/rqaUC )A
spread = 10.0;
^9Je8 @Yu
2\,vq R
clock_t start, finish;
D(]])4
start = clock();
g}hR q%
/* outFile.open("F:\\ez[40][25][25].txt");*/
s a71Vh{
for ( T = 1; T < 1001; T++ )
jSj (ZU6
{
G)E#wh_S^
/* Calculate the incident buffer */
u/h!i @_w[
for ( j=1; j<JE; j++)
|F'k5Lh
{
e!5nz_J1}
ez_inc[j] = ez_inc[j] + .5*(hx_inc[j-1] - hx_inc[j]);
1Jx|0YmO
}
0*.> >rI
//Fourier Transform of the incident field
5 ^\f[}
ine[0] = ine[0] + dt*cos(2*pi*f*dt*T)*ez_inc[ja-1];
W ,^(FR.
ine[1] = ine[1] - dt*sin(2*pi*f*dt*T)*ez_inc[ja-1];
=lf&mD _/
/* Source */
&dV|~xA6N
//pulse = exp (-.5*(pow(( t0-T )/spread,2)));
-pRyN]YD
pulse = exp (2*pi*f*T*dt);
?=]*r >a3
ez_inc[3] = pulse;
Q.Kr;64G
/* Boundary conditions for the incident buffer */
mYudUn4Wo
ez_inc[0] = ez_low_m2;
g8{?;
ez_low_m2 = ez_low_m1;
"DFj4XKXY9
ez_low_m1 = ez_inc[1];
tN5b rf
ez_inc[JE-1] = ez_high_m2;
cJ%u&2J_
ez_high_m2 = ez_high_m1;
j;O{Hvvz
ez_high_m1 = ez_inc[JE-2];
kd9GHN ;7
/* Calculate the Dx field */
04 2sjt
for ( i=1; i<ia; i++)
l_!.yV{
{
U/B1/96 lJ
for ( j=1; j<JE; j++)
up~l4]b+
{
NAj1ORy4pX
for (k=1; k<KE; k++)
X^#.4:>.
{
,T{(t@
curl_h = ( hz
[j][k] - hz
[j-1][k] - hy
[j][k] + hy
[j][k-1] );
TXi$Q%0W
idxl
[j][k] = idxl
[j][k] + curl_h;
{V!Jj6n
dx
[j][k] = gj3[j]*gk3[k]*dx
[j][k] + gj2[j]*gk2[k]*.5*(curl_h + gi1
*idxl
[j][k]);
eXa a'bTx
}
/ 4K *iq
}
_[SP*" ]H
}
AcN~Q/xU
for ( i=ia; i<=ib; i++)
eo4<RDe<
{
Q r;es,f
for ( j=1;j<JE; j++)
2O|o%` ?
{
N b(f
for ( k=1; k<KE; k++)
v UAYYe
{
$ a6&OH/
curl_h = ( hz
[j][k] - hz
[j-1][k] - hy
[j][k] + hy
[j][k-1] );
C61KY7iyR
dx
[j][k] = gj3[j]*gk3[k]*dx
[j][k] + gj2[j]*gk2[k]*.5*curl_h;
N1UE u,j
}
Y.?|[x0Wh
}
Y 962rZ
}
/FiFtAbb
for ( i=ib+1; i<IE; i++)
J ,fXXi)J
{
FeS6>/
ixh = i-ib-1;
N1Y*IkW"
for ( j=1; j<JE; j++)
R0fZ9_d7}
{
E jB<`yT
for ( k=1; k<KE; k++)
H{x}gBQ
{
&L-y1'i=j
curl_h = ( hz
[j][k] - hz
[j-1][k] - hy
[j][k] + hy
[j][k-1] );
o(~>a
idxh[ixh][j][k] = idxh[ixh][j][k] + curl_h;
:&`Yz
dx
[j][k] = gj3[j]*gk3[k]*dx
[j][k] + gj2[j]*gk2[k]*.5*(curl_h + gi1
*idxh[ixh][j][k]);
q<}5KY
}
5c P]
}
sp[nKo ^
}
gv; =Yhw.c
/* Calculate the Dy field */
fPHv|_XM>
for ( i=1; i<IE; i++)
cW),Y|8
{
:Qu!0tY
for ( j=1; j<ja; j++)
4j)Y>
{
3Fh<%<=
for (k=1; k<KE; k++)
5/mW:G,&
{
)%C482GO-
curl_h = ( hx
[j][k] - hx
[j][k-1] - hz
[j][k] + hz[i-1][j][k] );
C%v@ u$N
idyl
[j][k] = idyl
[j][k] + curl_h;
C["^%0lj
dy
[j][k] = gi3
*gk3[k]*dy
[j][k]+gi2
*gk2[k]*.5*(curl_h+ gj1[j]*idyl
[j][k]);
dH!k {3bL
}
g?(Z+w4A 3
}
4(&00#Yxg2
}
5SX 0g(C
for ( i=1; i<IE; i++)
C4Q ^WU+$j
{
"Q+'lA[}
for ( j=ja;j<=jb; j++)
<P( K,L?r
{
H(]lqvO
for ( k=1; k<KE; k++)
neQ2 +W%oj
{
E]_lYYkA
curl_h = ( hx
[j][k] - hx
[j][k-1] - hz
[j][k] + hz[i-1][j][k] );
&I?1(t~hT
dy
[j][k] = gi3
*gk3[k]*dy
[j][k]+gi2
*gk2[k]*.5*curl_h;
)xP]rOT
}
ZzI^*Nyg
}
@P "q`*
}
p7 !q#o
for ( i=1; i<IE; i++)
g i:;{
{
XFmnZpqXH
for ( j=jb+1; j<JE; j++)
!cnH|ePbI
{
YjS|Ht- >
jyh = j-jb-1;
.F)b9d[?
for ( k=1; k<KE; k++)
Lq LciD
{
Q=gVx S
curl_h = ( hx
[j][k] - hx
[j][k-1] - hz
[j][k] + hz[i-1][j][k] );
[_qBp:_j?s
idyh
[jyh][k] = idyh
[jyh][k] + curl_h;
gSu+]N
dy
[j][k] = gi3
*gk3[k]*dy
[j][k]+gi2
*gk2[k]*.5*(curl_h+ gj1[j]*idyh
[jyh][k]);
%M9^QHyo@
}
D}!U?]la&
}
M.d{:&@`%
}
kOR%<#:J
/* Incident Dy */
^k^%w/fo
for ( i=ia; i<=ib; i++)
K]@^8e$(
{
U*k$pp6\b~
for ( j=ja; j<=jb-1; j++)
"S[VtuxPCU
{
"SyyOD )WA
dy
[j][ka] = dy
[j][ka] - .5*hx_inc[j];
xzMa[D4(
dy
[j][kb+1] = dy
[j][kb+1] + .5*hx_inc[j];
"=| yM~V
}
WLNkO^zb
}
"6pjkE t4
/* Calculate the Dz field */
-.g5|B
for ( i=1; i<IE; i++)
2Pi}<pG~
{
J~<:y Bup}
for ( j=1; j<JE; j++)
4pq >R
{
I|9e4EX{ y
for (k=1; k<ka; k++)
)J (ekfM
{
3Luv$6
curl_h = ( hy
[j][k] - hy[i-1][j][k] - hx
[j][k] + hx
[j-1][k] );
| fTQ\q]W
idzl
[j][k] = idzl
[j][k] + curl_h;
pq-zy6^
dz
[j][k] = gi3
*gj3[j]*dz
[j][k]+gi2
*gj2[j]*.5*(curl_h+ gk1[k]*idzl
[j][k]);
,X\ z# B
}
AGCqJ 8`|T
}
5yxZ 5Ni!
}
EE&~D~yHUL
for ( i=1; i<IE; i++)
wC=IN
{
v_5DeaMF'
for ( j=1;j<JE; j++)
% C 6 H(
{
~v ,LFIT
for ( k=ka; k<=kb; k++)
Rl3KE)<
{
pPh_p @3I
curl_h = ( hy
[j][k] - hy[i-1][j][k] - hx
[j][k] + hx
[j-1][k] );
]pb;q(?^
dz
[j][k] = gi3
*gj3[j]*dz
[j][k]+gi2
*gj2[j]*.5*curl_h;
FNmIXpAn*@
}
d-N<VVcy\
}
N4fuV?E`
}
LC[, K
for ( i=1; i<IE; i++)
giaO7Qh~
{
<ze' o.c
for ( j=1; j<JE; j++)
aQFY Sl
{
6/@ cP /
for ( k=kb+1; k<KE; k++)
$am7 xd
{
4)'5;|pI
kzh = k-kb-1;
Di-"y, [
curl_h = ( hy
[j][k] - hy[i-1][j][k] - hx
[j][k] + hx
[j-1][k] );
.# j)YG
idzh
[j][kzh] = idzh
[j][kzh] + curl_h;
,Z>Rv Ll
dz
[j][k] = gi3
*gj3[j]*dz
[j][k]+gi2
*gj2[j]*.5*(curl_h+ gk1[k]*idzh
[j][kzh]);
Uk;S Y[mU
}
f^%v IB ~[
}
o{UwUMw5`
}
3O#7OL68v
/* Incident Dz */
mvUYp,JECl
for ( i=ia; i<=ib; i++)
R"O9 ~s6N
{
1P2% n [y
for ( k=ka; k<=kb; k++)
G0]q(.sOy
{
/B1< N}
dz
[ja][k] = dz
[ja][k] + .5*hx_inc[ja-1];
p2&KGt X'
dz
[jb][k] = dz
[jb][k] - .5*hx_inc[jb];
h.2!d0j]
}
c}v:X Slh7
}
{b/AOR o
/* Calculate the E from D field */
Z"!C
for ( i=1; i<IE-1; i++)
6Mk@ ,\1
{
?m?e2{]u,
for ( j=1; j<JE-1; j++)
C(,s_Ks
{
'd6hQ4Vw4
for ( k=1; k<KE-1; k++)
yuI 5# VUS
{
Uut,cQ". d
ex
[j][k] = gax
[j][k]*(dx
[j][k]-ix
[j][k]);
TF=S \ Q
ix
[j][k] = ix
[j][k]+gbx
[j][k]*ex
[j][k];
u3k+ Xg:
ey
[j][k] = gay
[j][k]*(dy
[j][k]-iy
[j][k]);
YuPgsJ[m
iy
[j][k] = iy
[j][k]+gby
[j][k]*ey
[j][k];
'a"<uk3DT
ez
[j][k] = gaz
[j][k]*(dz
[j][k]-iz
[j][k]);
Ujj2A^
iz
[j][k] = iz
[j][k]+gbz
[j][k]*ez
[j][k];
;xth#j
}
5>A3;P
}
UdL`. D ,
}
%],.?TS2V
/* Remember: part of the PML is E = 0 at the edges */
jjl4A} *0
for (j = 0;j < JE-1;j++)
-BoN}xE4
{
j*xe ns$)
for (k = 0;k < KE-1;k++)
Iv/yIS
{
] E=JUYf0
ez[0][j][k] = 0.0;
&# @1n
ez[IE-1][j][k] = 0.0;
})!n1k t
}
^x/0*t5};z
}
3tLh{S?uJ
for (i = 0;i < IE-1;i++)
a- QHm;_S
{
B r pin
for (k = 0;k < KE;k++)
bjQfZT (
{
Vx0V6{JX
ez
[0][k] = 0.0;
u: ,B"!
ez
[JE-1][k] = 0.0;
j/B zbjq"
}
*Z/B\n b
}
O _1}LS!
for (i = 0;i < IE-1;i++)
t 8M3V GN
{
X^xu$d6
for (j = 0;j < JE-1;j++)
DK)q Bxc8
{
Q?1 K xD!
ez
[j][0] = 0.0;
w$B7. .r
ez
[j][KE-1] = 0.0;
**s:H'M w_
}
#4c uNX5m%
}
}qlz^s
Yhlk#>I
/*outFile << setw(22) << setprecision(8) << T << setw(22) << setprecision(8) << ez[40][25][25]<< endl;*/
Y_)04dmr@[
h/_z QR-
/* Calculate the incident field */
0IA '8_K
for ( j=0; j<JE-1; j++)
=G=.THRUk
{
mZM5aTQ3
hx_inc[j] = hx_inc[j] + .5*(ez_inc[j]-ez_inc[j+1]);
K8 ThZY%
}
W%h<@@c4,
/* Calculate the Hx field */
KqtI^qC8
for ( i=0; i<ia; i++)
7< 1Y%|x`
{
BAKfs/N
for ( j=0; j<JE-1; j++)
s{iYf :
{
gJF;yW 4
for (k=0; k<KE-1; k++)
}%XNB1/`
{
'QW 0K]il
curl_e = ( ey
[j][k+1] - ey
[j][k] - ez
[j+1][k] + ez
[j][k] );
CWDo_g $
ihxl
[j][k] = ihxl
[j][k] + curl_e;
%5z88-\
hx
[j][k] = fj3[j]*fk3[k]*hx
[j][k]+fj2[j]*fk2[k]*.5*(curl_e+ fi1
*ihxl
[j][k]);
+ZZiZ&y
}
+ubO-A?
}
?xG #4P<C=
}
~RMOEH.o
for ( i=ia; i<=ib; i++)
/? J_7Lg
{
=3J~ F k
for ( j=0;j<JE-1; j++)
u YH{4%
{
xMuy[)b
for ( k=0; k<KE-1; k++)
R/^ rh
{
G&n_vwZ%
curl_e = ( ey
[j][k+1] - ey
[j][k] - ez
[j+1][k] + ez
[j][k] );
UNhD
hx
[j][k] = fj3[j]*fk3[k]*hx
[j][k]+fj2[j]*fk2[k]*.5*curl_e;
_` D_0v(X
}
[yVU p+
}
K1+,y1c
}
xH L{3^
for ( i=ib+1; i<IE; i++)
j5 /H#_ .
{
J J3 vC
ixh = i-ib-1;
i&bttSRNV
for ( j=0; j<JE-1; j++)
XM=`(e o
{
{u5)zVYC,U
for ( k=0; k<KE-1; k++)
P)}:lTe
{
]mqB&{g
curl_e = ( ey
[j][k+1] - ey
[j][k] - ez
[j+1][k] + ez
[j][k] );
^~A r
ihxh[ixh][j][k] = ihxh[ixh][j][k] + curl_e;
t(_XB|AKm
hx
[j][k] = fj3[j]*fk3[k]*hx
[j][k]+fj2[j]*fk2[k]*.5*(curl_e+ fi1
*ihxh[ixh][j][k]);
"thu@~aC
}
1zp ,Suv
}
Ic r'l$PE
}
CRqa[boU*
/* Incident Hx */
` 8 6b
for ( i=ia; i<=ib; i++)
n1;V2k{uV
{
<A6<q&g|E
for ( k=ka; k<=kb; k++)
$ u"K1Q 3
{
Xg%zE
hx
[ja-1][k] = hx
[ja-1][k] + .5*ez_inc[ja];
( h%wO
hx
[jb][k] = hx
[jb][k] - .5*ez_inc[jb];
9B%"7MVn
}
GdP9Uj)n-
}
}3i@5ctQ
.[_&>@bmrP
/* Calculate the Hy field */
dm 2_Fj
for ( i=0; i<IE-1; i++)
i L\\JuY
{
K8RloDjk_A
for ( j=0; j<ja; j++)
j1d=$'a "
{
> voUh;L
for (k=0; k<KE-1; k++)
v~xG *e
{
%%3ugD5i!
curl_e = ( ez[i+1][j][k] - ez
[j][k] - ex
[j][k+1] + ex
[j][k] );
lq9c2xK
ihyl
[j][k] = ihyl
[j][k] + curl_e;
u C}YKT>V7
hy
[j][k] = fi3
*fk3[k]*hy
[j][k]+fi2
*fk2[k]*.5*(curl_e+ fj1[j]*ihyl
[j][k]);
~O|0.)71]
}
x)GoxH~#
}
>:Oo[ {)
}
|LjCtm)@+
for ( i=0; i<IE-1; i++)
My1E@<
{
!nYAyjf
for ( j=ja;j<=jb; j++)
~ \kRW6
{
P] x@h
for ( k=0; k<KE-1; k++)
M,{F/Yu
{
J$P]>By5:
curl_e = ( ez[i+1][j][k] - ez
[j][k] - ex
[j][k+1] + ex
[j][k] );
nv8,O= #s
hy
[j][k] = fi3
*fk3[k]*hy
[j][k]+fi2
*fk2[k]*.5*curl_e;
aY?}4Bx
}
lW2q V R
}
4t+88e
}
>AfJxdd1
for ( i=0; i<IE-1; i++)
1ii.nt1 u
{
']rh0?
for ( j=jb+1; j<JE; j++)
7u}r^+6_o
{
|Y99s)2&N
jyh = j-jb-1;
EfUo<E
for ( k=0; k<KE-1; k++)
]pGr'T~Gj
{
ilHj%h*z
curl_e = ( ez[i+1][j][k] - ez
[j][k] - ex
[j][k+1] + ex
[j][k] );
~/Ry=8
ihyh
[jyh][k] = ihyh
[jyh][k] + curl_e;
o^NQ]BdH8
hy
[j][k] = fi3
*fk3[k]*hy
[j][k]+fi2
*fk2[k]*.5*(curl_e+ fj1[j]*ihyh
[jyh][k]);
Y/hay[ 6
}
/2U.,vw
}
G.N3 R
}
JBg>E 3*N
/* Incident Hy */
ofi']J{R
for ( j=ja; j<=jb; j++)
}@#e D
{
< l6CtK@
for ( k=ka; k<=kb; k++)
C[Fh^
{
npe *A
hy[ia-1][j][k] = hy[ia-1][j][k] - .5*ez_inc[j];
Vw#07P#A
hy[ib][j][k] = hy[ib][j][k] + .5*ez_inc[j];
5 QT9
}
m7@`POI
}
iN)@Cu7
/* Calculate the Hz field */
k |Lm;g
for ( i=0; i<IE-1; i++)
'- ~86Q
{
eQeNlCG
for ( j=0; j<JE-1; j++)
q)vD "{ 0.
{
: L?zk"0C
for (k=0; k<ka; k++)
,L^L uw'7
{
*X>rvAd3
curl_e = ( ex
[j+1][k] - ex
[j][k] - ey[i+1][j][k] + ey
[j][k] );
V{*9fB#4 L
ihzl
[j][k] = ihzl
[j][k] + curl_e;
z/T ZOFaM
hz
[j][k] = fi3
*fj3[j]*hz
[j][k]+fi2
*fj2[j]*.5*(curl_e+ fk1[k]*ihzl
[j][k]);
ht-6_]+ME
}
'nLv0.7*
}
K~E]Fkw!;
}
W"0 #
for ( i=0; i<IE-1; i++)
C 8jZcs#4
{
7oWv'
for ( j=0;j<JE-1; j++)
q\/|nZO4
{
`2Z= Lp
for ( k=ka; k<=kb; k++)
*V\kS
{
ruW6cvsvet
curl_e = ( ex
[j+1][k] - ex
[j][k] - ey[i+1][j][k] + ey
[j][k] );
ULJI` I|m
hz
[j][k] = fi3
*fj3[j]*hz
[j][k]+fi2
*fj2[j]*.5*curl_e;
ywdNwNJ
}
`=>Bo p)
}
0O:TKgb&C.
}
DU tpd|
for ( i=0; i<IE-1; i++)
8[Qw8z5-
{
T]myhNk
for ( j=0; j<JE-1; j++)
C6Mb(&
{
X >**M
for ( k=kb+1; k<KE; k++)
{u1t .+
{
v{jl)?`~w
kzh = k-kb-1;
2<O hO ^
curl_e = ( ex
[j+1][k] - ex
[j][k] - ey[i+1][j][k] + ey
[j][k] );
.xD-eWw3R
ihzh
[j][kzh] = ihzh
[j][kzh] + curl_e;
V&>7i9lEz
hz
[j][k] = fi3
*fj3[j]*hz
[j][k]+fi2
*fj2[j]*.5*(curl_e+ fk1[k]*ihzh
[j][kzh]);
