登 录
註 冊
论坛
微波仿真网
注册
登录论坛可查看更多信息
微波仿真论坛
>
时域有限差分法 FDTD
>
大家帮忙看看这个3D_FDTD介质球的RCS有什么问 ..
发帖
回复
1
2
3
4123
阅读
29
回复
大家帮忙看看这个3D_FDTD介质球的RCS有什么问题
离线
cc81372365
UID ：27752
注册：
2009-03-18
登录：
2023-11-09
发帖：
345
等级：
论坛版主
0楼
发表于: 2010-11-29 13:16:05
//****************************************************************//
.R8 HZ}3
//
[kzd (u
//********* 3D_介质球散射（RCS）FDTD程序 **********//
,\d6VB P&
//
2Nm>5l
//****************************************************************//
|R@~-Ht
//
I3:[= ,5
//入射波为Y方向的平面波（theta=pi/2 , phi=pi/2), RCS观察角为（theta=pi/2 , phi = 0到2pi）
he -Ji
// 数字0，1代表场值经过傅里叶变换后的复数的实部和虚部
p,<&zHb >K
// 金属方柱置于中心位置（ic,jc）,激励源为高斯脉冲源
`)h6j)xiQ
#include <iostream>
q;D+ai
#include <iomanip>
7HJS .047
#include <fstream>
<uk1 ?Q g
//#include <omp.h>
{(#%N5%
#include <time.h>
G"59cv8z4R
#include <valarray>
Af5D>/
#include <math.h>
CBKkBuKuk
#include <stdio.h>
(ihP `k-.
#include <stdlib.h>
lvpc*d|K
q )lnS )
#include "complex_operation.h"
70Yjv 1 i
#define IE 60
5P h X"7
#define JE 60
!.<T"8BUpv
#define KE 60
h~,JdDV8l*
#define ia 15
J3b4cxm
#define ja 15
-/dEsgO
#define ka 15
b7\ cxgRq
#define length 41
#4h+ j%y[H
using namespace std;
;t> Z +O%
ofstream outFile;
omf Rs
int main()
$]&(7@'qo
{
4YMX|1wd)
int i,j,k,ic,jc,kc,ib,jb,kb;
Tv]<SI<B[
int npml,n_pml,T;
sz):oea@f@
double radius,v,lamda,Z,f,kk;
>x@P|\
double epsz_r,epsz,muz,muz_r,ddx,dt,pi,sigma;
#W 2[
double xn,xxn;
m&0"<V!H/B
double t0,spread,pulse;
y3;q_4.
double ez_inc[JE]={0},hx_inc[JE]={0};
Koln9'tB
double ez_low_m1=0,ez_low_m2=0,ez_high_m1=0,ez_high_m2=0;
o1OBwPj
int ixh,jyh,kzh;
Q2%QLM:.,
int co=0;
zT* .jv
double rcs[361]={0.0};
,kp\(X[J
double ine[2] = {0.0};
w_4]xgS:
double ine_fudu = 0.0;
RF!1 oZ
double curl_h,curl_e;
^, i>'T
double dist,xdist,ydist,zdist;
rf 9 _eP
double dx[IE][JE][KE]={0},dy[IE][JE][KE]={0},dz[IE][JE][KE]={0};
G>Em! 4h
double ex[IE][JE][KE]={0},ey[IE][JE][KE]={0},ez[IE][JE][KE]={0};
w$5A|%Y+V}
double hx[IE][JE][KE]={0},hy[IE][JE][KE]={0},hz[IE][JE][KE]={0};
Dli^2hD
double ix[IE][JE][KE]={0},iy[IE][JE][KE]={0},iz[IE][JE][KE]={0};
z ggB$ 5
double gax[IE][JE][KE]={0},gay[IE][JE][KE]={0},gaz[IE][JE][KE]={0};
QIn/,Yd
double gbx[IE][JE][KE]={0},gby[IE][JE][KE]={0},gbz[IE][JE][KE]={0};
i=32KI(%
double idxl[ia][JE][KE]={0},idxh[ia][JE][KE]={0};
MZSxQ 8
double idyl[IE][ja][KE]={0},idyh[IE][ja][KE]={0};
x2#qg>`l
double idzl[IE][JE][ka]={0},idzh[IE][JE][ka]={0};
W~B5>;y
double ihxl[ia][JE][KE]={0},ihxh[ia][JE][KE]={0};
Lx U= {Y0
double ihyl[IE][ja][KE]={0},ihyh[IE][ja][KE]={0};
Drvt H+e
double ihzl[IE][JE][ka]={0},ihzh[IE][JE][ka]={0};
q 'a
double gi1[IE]={0},gi2[IE]={0},gi3[IE]={0};
>(tn "2
double gj1[JE]={0},gj2[JE]={0},gj3[JE]={0};
sK=}E=
double gk1[KE]={0},gk2[KE]={0},gk3[KE]={0};
OAZ#|U
double fi1[IE]={0},fi2[IE]={0},fi3[IE]={0};
'7B"(dA&C
double fj1[JE]={0},fj2[JE]={0},fj3[JE]={0};
7W9d6i)
double fk1[KE]={0},fk2[KE]={0},fk3[KE]={0};
z N_:nY>
double Ex_down[length][length][2]={0},Ey_down[length][length][2]={0};
! sA_?2$
double Hx_down[length][length][2]={0},Hy_down[length][length][2]={0};
$O:w(U
double Ex_up[length][length][2]={0},Ey_up[length][length][2]={0};
i~{ _eQV
double Hx_up[length][length][2]={0},Hy_up[length][length][2]={0};
M}"r#Plq
double Ex_left[length][length][2]={0},Ez_left[length][length][2]={0};
0gF!!m
double Hx_left[length][length][2]={0},Hz_left[length][length][2]={0};
%im#ww L%
double Ex_right[length][length][2]={0},Ez_right[length][length][2]={0};
Q1&P@Io$
double Hx_right[length][length][2]={0},Hz_right[length][length][2]={0};
.`Zf} [5[
double Ey_front[length][length][2]={0},Ez_front[length][length][2]={0};
) I@L+
double Hy_front[length][length][2]={0},Hz_front[length][length][2]={0};
2o[IHO]
double Ey_back[length][length][2]={0},Ez_back[length][length][2]={0};
nx ap\Lf
double Hy_back[length][length][2]={0},Hz_back[length][length][2]={0};
td -3 h,\\
ic = IE/2;
qkP/Nl. u
jc = JE/2;
bv dR"G
kc = KE/2;
D 0]a \,aZ
ib = IE - ia -1;
;}.Kb
jb = JE - ja -1;
baoD(0d
kb = KE - ka -1;
>1irSUj"~
//外推边界位置
i(wgB\9i4
int dis = 5;
<B&R6<]T
const int i_front = ib + dis;
:H[\;Z1_
const int i_back = ia - dis;
qJT0Y/l:(
const int j_right = jb + dis;
9f}X Rz
const int j_left = ja - dis;
qq0?e0H
const int k_up = kb + dis;
Q I!i
const int k_down = ka - dis;
9<]a!:!^
pi = 3.1415926;
E)t
epsz = 8.8e-12; // 真空介电常数
N\1/ JW+
epsz_r = 4.0; // 相对介电常数
jB^OP1
muz = 4*pi*1e-7; // 真空磁导率
4`2$_T$ F
muz_r = 1.0; // 相对磁导率
=j /hl
sigma = 0; // 介质电导率
!m{2WW-
v = 3.0e8; // 波速
vV`|!5x
f = 9.375e9; // 频率 10GHz
88a<{5 :z
lamda = v/f; // 波长 0.03m
.Nx W=79t
kk = 2*pi/lamda; // 波数
zy N (4
Z = sqrt(muz/epsz); // 真空波阻抗 377
[ij,RE7,T
ddx = lamda/42; // 网格尺寸
-05U%l1e
dt = ddx/(2.0*v); // 时间步长
fRg=!<#%
/* Initialize the arrays */
re,.@${H
for ( i=0; i<IE; i++)
[~k]{[NJ
{
gw3NS8 A+
for ( j=0; j<JE; j++)
PG )_L.7rJ
{
2+9 2Q_ +
for ( k=0; k<KE; k++)
; a/cty0Ch
{
DcV<y-`'1
ex
[j][k] = 0.0;
F}rPY:
ey
[j][k] = 0.0;
NyI ;v =
ez
[j][k] = 0.0;
qGPb
dx
[j][k] = 0.0;
K{ }4zuZ
dy
[j][k] = 0.0;
2 6:evid
dz
[j][k] = 0.0;
02]x Jo
hx
[j][k] = 0.0;
/k7wwZiY@
hy
[j][k] = 0.0;
f'dK73Xof
hz
[j][k] = 0.0;
?9u4a_x
ix
[j][k] = 0.0;
s "0b%0?A
iy
[j][k] = 0.0;
JAn 1{<Ky
iz
[j][k] = 0.0;
! Zw f 397
gax
[j][k] = 1.0;
0v"&G<J
gay
[j][k] = 1.0;
LE c8NQs
gaz
[j][k] = 1.0;
1}`LTPW9
gbx
[j][k] = 0.0;
{B yn{?w
gby
[j][k] = 0.0;
{.#zHL ;
gbz
[j][k] = 0.0;
3BMS_,P
}
j/<??v4F4
}
:?r*p>0$
}
`h;}3r#R{
for ( i=0; i<ia; i++)
.:;fAJPf
{
H$- $2?5
for ( j=0; j<JE; j++)
u5gZxO1J5
{
^w6eW zI
for ( k=0; k<KE; k++)
d0C8*ifFO
{
sdF3cX
idxl
[j][k] = 0.0;
co B 6 rW
idxh
[j][k] = 0.0;
u2`xC4>c
ihxl
[j][k] = 0.0;
wt7.oKbW
ihxh
[j][k] = 0.0;
Z[,`"}}hv=
}
d7bjbJwu
}
.9\Cy4_qSd
}
m(3);)d
for ( i=0; i<IE; i++)
Ww87
{
'rV2Bt,
for ( j=0; j<ja; j++)
~440# kj<
{
.w FU:y4r
for ( k=0; k<KE; k++)
Tj,Nmb>Q7'
{
*=^[VV!
idyl
[j][k] = 0.0;
8SD}nFQ
idyh
[j][k] = 0.0;
K M]Wl_z
ihyl
[j][k] = 0.0;
uU+s!C9r
ihyh
[j][k] = 0.0;
: 2q ?>\
}
^L~ [+|
}
eECj_eH-
}
8)Tj H'
for ( i=0; i<IE; i++)
wd`R4CKhP]
{
< J# R3{
for ( j=0; j<JE; j++)
`QCD$=
{
#DaP=k"XV
for ( k=0; k<ka; k++)
f0F#Yi{fw
{
c=t*I0-OVS
idzl
[j][k] = 0.0;
"<dN9l>
idzh
[j][k] = 0.0;
rZ866\0
ihzl
[j][k] = 0.0;
s}b*5@8|tA
ihzh
[j][k] = 0.0;
4 RO Wz
}
FYeEG
}
t +}uIp42<
}
px&=((Z7>
//PML初始化
H*qD: N
for ( i=0; i<IE; i++)
gO{W#%
{
r|8V @ .@i
gi1
= 0;
w:tGPort
fi1
= 0;
$wXih#7
gi2
= 1;
3GWrn ,f
fi2
= 1;
QBj Y&(vY
gi3
= 1;
p[P[#IeL
fi3
= 1;
U5wT Gv4S|
}
2<.Vv\ =
for ( j=0; j<JE; j++)
'o8\`\'H!
{
0O ['w<_
gj1[j] = 0;
bf^ly6ml
fj1[j] = 0;
CI1m5g [P
gj2[j] = 1;
F9D"kG;Dk
fj2[j] = 1;
#HcI4j:s!
gj3[j] = 1;
(j`l5r#X#/
fj3[j] = 1;
xDe47&qKM
}
5fxbA2\
for ( k=0; k<KE; k++)
5k `e^ARf
{
H5 ^Y->
gk1[k] = 0;
'?L%F{g/9
fk1[k] = 0;
2vXGO|W
gk2[k] = 1;
w2<*$~C]
fk2[k] = 1;
bG1 of sU
gk3[k] = 1;
6(5c 7R#
fk3[k] = 1;
xf U hS t
}
Y=WR6!{
radius = 10 ;
= P8 ~n2V
n_pml = 7; //PML层厚度
0- Yeu5A
npml = 7;
W,'3D~g8
for (i = 0; i < n_pml; i++)
bf c.rZ
{
L (Rorf~V
xxn = (double)(npml-i)/(double)npml;
\1kh yF'
xn = 0.33*pow(xxn,3);
{ ;' :h
fi1
= xn;
+wjlAqMQ
fi1[IE-1-i] = xn;
p d%LL?O
gi2
= 1.0/(1.0 + xn);
iDvpXn
gi2[IE-1-i] = 1.0/(1.0 + xn);
7#/|VQX<A
gi3
= (1.0 - xn)/(1.0 + xn);
_9qEZV
gi3[IE-1-i] = (1.0 - xn)/(1.0 + xn);
$:HLRl{2E
xxn = (double)(npml-i-.5)/(double)npml;
I9s$bRbT
xn = .33*pow(xxn,3);
;W4:#/~14
gi1
= xn;
%[+/>e/m
gi1[IE-2-i] = xn;
'XME?H:q a
fi2
= 1.0/(1.0 + xn);
qSCTFJ0
fi2[IE-2-i] = 1.0/(1.0 + xn);
4}PeP ^pj
fi3
= (1.0 - xn)/(1.0 + xn);
4 jD\]Q="1
fi3[IE-2-i] = (1.0 - xn)/(1.0 + xn);
Uc>LFX& -B
}
$u- lo|
for (j = 0; j < n_pml; j++)
"$%{}{ #W0
{
{C, #rj
xxn = (double)(npml-j)/(double)npml;
z+2u- jG
xn = .33*pow(xxn,3);
IM|Se4;x
fj1[j] = xn;
OUKj@~T
fj1[JE-1-j] = xn;
IM2/(N.%
gj2[j] = 1.0/(1.0 + xn);
|3W3+Rn!
gj2[JE-1-j] = 1.0/(1.0 + xn);
s2%0#6c'c
gj3[j] = (1.0 - xn)/(1.0 + xn);
zVn* !c
gj3[JE-1-j] = (1.0 - xn)/(1.0 + xn);
> fhSaeN
xxn = (double)(npml-j-.5)/(double)npml;
2w;G4
xn = .33*pow(xxn,3);
6;[1Jz]?i
gj1[j] = xn;
pIrv$^
gj1[JE-2-j] = xn;
{K6Kx36
fj2[j] = 1.0/(1.0 + xn);
(#lm#?<)
fj2[JE-2-j] = 1.0/(1.0 + xn);
)<tzm'Rc
fj3[j] = (1.0 - xn)/(1.0 + xn);
+&zb^C`J
fj3[JE-2-j] = (1.0 - xn)/(1.0 + xn);
c[6 zX#{`
}
Y~}QJ +`?
.M`LUb"!
for (k = 0; k < n_pml; k++)
U0ns3LirP
{
@b>YkJDk
xxn = (double)(npml-k)/(double)npml;
\_)02ZT:
xn = .33*pow(xxn,3);
?\l!]vu*
fk1[k] = xn;
5! 2J;.&
fk1[KE-1-k] = xn;
V =Ww>
gk2[k] = 1.0/(1.0 + xn);
)OS>9 kFH
gk2[KE-1-k] = 1.0/(1.0 + xn);
f a/P%9d b
gk3[k] = (1.0 - xn)/(1.0 + xn);
($,iAb
gk3[KE-1-k] = (1.0 - xn)/(1.0 + xn);
</2,2AV4q*
xxn = (double)(npml-k-.5)/(double)npml;
~m3V]v(q7
xn = .33*pow(xxn,3);
c@)p Ki#W
gk1[k] = xn;
``/y=k/a u
gk1[KE-2-k] = xn;
fjF!>Dy
fk2[k] = 1.0/(1.0 + xn);
B1up ^ (?
fk2[KE-2-k] = 1.0/(1.0 + xn);
WCWSLEAz a
fk3[k] = (1.0 - xn)/(1.0 + xn);
ejDCmD
fk3[KE-2-k] = (1.0 - xn)/(1.0 + xn);
h7)VJY
}
kz3?j<
X~` . }
/* Calculate gax ,gbx */
.) ?2)Fl
for ( i = ia; i<=ib; i++)
8dY k 3 sk
{
/,-h%gj
for ( j = ja; j<=jb; j++)
oT$(<$&<
{
MJpP!a^Q
for ( k =ka; k<=kb; k++)
_-YL!oP
{
I}I}K~se*
xdist = (ic-i-.5);
F?!};~$=Z
ydist = (jc-j);
$&c<T4 $d
zdist = (kc-k);
I>(;bNgN E
dist = sqrt (pow(xdist,2) + pow(ydist,2) + pow(zdist,2));
P((S2"D<4
if ( dist <= radius )
hG< a
{
A;b=E[i v
gax
[j][k] = 1./(epsz_r + (sigma*dt/epsz));
GC,vQ\
gbx
[j][k] = sigma*dt/epsz;
`,hW;p> -
}
5 >0\e_V
}
AD0ptHUBa
}
wG Z>iLe:
}
ZJ)3GF}4
/* Calculate gay ,gby*/
xop-f#U*
for ( i = ia; i<=ib; i++)
cS. 7\0$
{
-%7Jj;yA
for ( j = ja; j<=jb; j++)
~t1O]aO(
{
{-:4O\/
for ( k =ka; k<=kb; k++)
_m) gO/02A
{
RcKQER
xdist = (ic-i);
#3A Yz82w
ydist = (jc-j-.5);
n$}R/*
zdist = (kc-k);
& bp#1KR)
dist = sqrt (pow(xdist,2) + pow(ydist,2) + pow(zdist,2));
3W%f#d$`
if ( dist <= radius )
sk i1f
{
OcyiL)tv 5
gay
[j][k] = 1./(epsz_r + (sigma*dt/epsz));
A8CIP:Z
gby
[j][k] = sigma*dt/epsz;
qBf wN 1
}
#r78Ym'aI
}
.eZPp~[lAN
}
ym- lT|>Z
}
Ew)n~!s
/* Calculate gaz ,gbz */
j=!(F`/
for ( i = ia; i<=ib; i++)
YMd&To 0s
{
g]oc(RM
for ( j = ja; j<=jb; j++)
HMl!?%%
{
^\Ue7,H-
for ( k =ka; k<=kb; k++)
OtrXYiKB
{
:e5:\|5*5
xdist = (ic-i);
*B)Jv9
ydist = (jc-j);
w%%6[<3%
zdist = (kc-k-.5);
FkB6*dm-
dist = sqrt (pow(xdist,2) + pow(ydist,2) + pow(zdist,2));
EN5G:hD
if ( dist <= radius )
{Zd)U "
{
;C7BoHB9
gaz
[j][k] = 1./(epsz_r +(sigma*dt/epsz));
UeutFNp
gbz
[j][k] = sigma*dt/epsz;
\#IJ=+z
}
n0>5'm%ES
}
;p?42rCIcl
}
`[g# Mxw
}
(25^ r
//源参数
1"~O"m sb
t0 = 40.0;
ZEXj |w C
spread = 10.0;
/*mFP.en
tk] _QX %
clock_t start, finish;
%M4XbSN|
start = clock();
zz+M1n-;o
/* outFile.open("F:\\ez[40][25][25].txt");*/
8iII) +
for ( T = 1; T < 1001; T++ )
}R?v"6aBS
{
QiQ2XW\E
/* Calculate the incident buffer */
T<9 dW?'|
for ( j=1; j<JE; j++)
wz|Q%.%?[
{
DkF@XK0c3
ez_inc[j] = ez_inc[j] + .5*(hx_inc[j-1] - hx_inc[j]);
Wme1Uid
}
)- Wn'C' Z
//Fourier Transform of the incident field
>Rz#g*@E
ine[0] = ine[0] + dt*cos(2*pi*f*dt*T)*ez_inc[ja-1];
Wfi:wCqZG
ine[1] = ine[1] - dt*sin(2*pi*f*dt*T)*ez_inc[ja-1];
71}L# nQ
/* Source */
\]~kyy
//pulse = exp (-.5*(pow(( t0-T )/spread,2)));
[TpA26#TTO
pulse = exp (2*pi*f*T*dt);
@[[ C s*-
ez_inc[3] = pulse;
|zRoXO`]-*
/* Boundary conditions for the incident buffer */
xC=3| ,U
ez_inc[0] = ez_low_m2;
TV[6+i*#
ez_low_m2 = ez_low_m1;
6cg pg+-a
ez_low_m1 = ez_inc[1];
5KA FUR0
ez_inc[JE-1] = ez_high_m2;
Rd;~'gbG
ez_high_m2 = ez_high_m1;
y5Z<uwXc
ez_high_m1 = ez_inc[JE-2];
*h5ld P
/* Calculate the Dx field */
3f7t%
for ( i=1; i<ia; i++)
*1 J #Mdd
{
Hz;jJ&S
for ( j=1; j<JE; j++)
0|w KR|z W
{
nEZ-h7lzl(
for (k=1; k<KE; k++)
m,"cbJ /
{
\M3Nas Z
curl_h = ( hz
[j][k] - hz
[j-1][k] - hy
[j][k] + hy
[j][k-1] );
%i]uW\~U
idxl
[j][k] = idxl
[j][k] + curl_h;
_ z"ci$[
dx
[j][k] = gj3[j]*gk3[k]*dx
[j][k] + gj2[j]*gk2[k]*.5*(curl_h + gi1
*idxl
[j][k]);
y:^>(l #;
}
lN= m$ J
}
GakmROZ@9
}
"\R@l Ux.Y
for ( i=ia; i<=ib; i++)
f6aT[Nw<
{
H*:r>Lm=
for ( j=1;j<JE; j++)
<(6-9(zHa
{
z wniS6R1
for ( k=1; k<KE; k++)
L`VQ{|&3V
{
y [ rB"
curl_h = ( hz
[j][k] - hz
[j-1][k] - hy
[j][k] + hy
[j][k-1] );
B [U.CAUn
dx
[j][k] = gj3[j]*gk3[k]*dx
[j][k] + gj2[j]*gk2[k]*.5*curl_h;
@ po MK:
}
sWpRX2 {5,
}
)sz 2 9
}
;)bF# @Q
for ( i=ib+1; i<IE; i++)
!m/Dd0
{
6LF^[b/u
ixh = i-ib-1;
L2V $%*6
for ( j=1; j<JE; j++)
ys"mP* wD
{
%+j]vP
for ( k=1; k<KE; k++)
x:&L?eOT
{
4.J aw+
curl_h = ( hz
[j][k] - hz
[j-1][k] - hy
[j][k] + hy
[j][k-1] );
c+G :@%
idxh[ixh][j][k] = idxh[ixh][j][k] + curl_h;
(VF4FC
dx
[j][k] = gj3[j]*gk3[k]*dx
[j][k] + gj2[j]*gk2[k]*.5*(curl_h + gi1
*idxh[ixh][j][k]);
V+"*A
}
T/spUlWu
}
Vg C9' "|
}
EL)/5-=S
/* Calculate the Dy field */
}IalgQ(i
for ( i=1; i<IE; i++)
K:lT- *+S
{
>sl1 cC
for ( j=1; j<ja; j++)
U K]{ ]-
{
^F{)&# 4
for (k=1; k<KE; k++)
vSHIl"h
{
8d*<Aki?;
curl_h = ( hx
[j][k] - hx
[j][k-1] - hz
[j][k] + hz[i-1][j][k] );
f >, Qhl
idyl
[j][k] = idyl
[j][k] + curl_h;
VdN+~+A:
dy
[j][k] = gi3
*gk3[k]*dy
[j][k]+gi2
*gk2[k]*.5*(curl_h+ gj1[j]*idyl
[j][k]);
TckR_0LNV
}
6jy n,GU
}
4- ?`#
}
+ke42Jwt
for ( i=1; i<IE; i++)
'zD ;:wT
{
lDX&v$
for ( j=ja;j<=jb; j++)
=QxE -)v
{
lLwQridFXh
for ( k=1; k<KE; k++)
0Ts_"p
{
<=GzK:4L
curl_h = ( hx
[j][k] - hx
[j][k-1] - hz
[j][k] + hz[i-1][j][k] );
f'?6D+Yw~
dy
[j][k] = gi3
*gk3[k]*dy
[j][k]+gi2
*gk2[k]*.5*curl_h;
'%|20 j
}
q0KXu MK
}
tRrY)eElS
}
K:mL%o2J
for ( i=1; i<IE; i++)
l4B O@
{
Xy(SzJ %
for ( j=jb+1; j<JE; j++)
Xt a>
{
%n`iA7j$W
jyh = j-jb-1;
aK=3`q
for ( k=1; k<KE; k++)
7<C~D,x6
{
qKb- aP-
curl_h = ( hx
[j][k] - hx
[j][k-1] - hz
[j][k] + hz[i-1][j][k] );
pl^"1 Z=*
idyh
[jyh][k] = idyh
[jyh][k] + curl_h;
] F) -}
dy
[j][k] = gi3
*gk3[k]*dy
[j][k]+gi2
*gk2[k]*.5*(curl_h+ gj1[j]*idyh
[jyh][k]);
@x>$_:]
}
/>j+7ts
}
C ,e$g
}
KA*l6 `(
/* Incident Dy */
YC,.Y{oY{
for ( i=ia; i<=ib; i++)
i!+3uHWu`)
{
PTc\I
for ( j=ja; j<=jb-1; j++)
(A<sFw?
{
kBQenMm
dy
[j][ka] = dy
[j][ka] - .5*hx_inc[j];
_F$t#.o
dy
[j][kb+1] = dy
[j][kb+1] + .5*hx_inc[j];
&.bR1wX
}
r :Ma AT<
}
}W>[OY0^A
/* Calculate the Dz field */
`] dx%
for ( i=1; i<IE; i++)
F8r455_W"
{
:$Di.|l@7
for ( j=1; j<JE; j++)
%dWFg<< |
{
a.*j8T
for (k=1; k<ka; k++)
tH|Q4C
{
OD!CnK
curl_h = ( hy
[j][k] - hy[i-1][j][k] - hx
[j][k] + hx
[j-1][k] );
Gy Xs{*
idzl
[j][k] = idzl
[j][k] + curl_h;
=K< I)2
dz
[j][k] = gi3
*gj3[j]*dz
[j][k]+gi2
*gj2[j]*.5*(curl_h+ gk1[k]*idzl
[j][k]);
z%gt V'
}
Rb>RjHo S
}
hq[ gj?P
}
nJ0eZBgB]
for ( i=1; i<IE; i++)
'; T5[l,
{
~esEql=Q3'
for ( j=1;j<JE; j++)
< }'B-k9
{
p 'c<v)ia
for ( k=ka; k<=kb; k++)
r D!.N
{
\ f Fy$
curl_h = ( hy
[j][k] - hy[i-1][j][k] - hx
[j][k] + hx
[j-1][k] );
3Os3=Ix
dz
[j][k] = gi3
*gj3[j]*dz
[j][k]+gi2
*gj2[j]*.5*curl_h;
`h{mj|~
}
8&[<pbN)
}
rZCAj
}
"Ohpb!J9
for ( i=1; i<IE; i++)
ydFhw}1>
{
#1hz=~YO
for ( j=1; j<JE; j++)
pj-HLuZR
{
R~ c vml
for ( k=kb+1; k<KE; k++)
H5MAN,`
{
TOF62,
kzh = k-kb-1;
<XcMc<h~
curl_h = ( hy
[j][k] - hy[i-1][j][k] - hx
[j][k] + hx
[j-1][k] );
TdOWdPvYj
idzh
[j][kzh] = idzh
[j][kzh] + curl_h;
b 0x0CMf
dz
[j][k] = gi3
*gj3[j]*dz
[j][k]+gi2
*gj2[j]*.5*(curl_h+ gk1[k]*idzh
[j][kzh]);
< |.! Px86
}
NF.6(PG|
}
)tQ6rd'
}
Bst>9V&R
/* Incident Dz */
*SG2k .$
for ( i=ia; i<=ib; i++)
KGLhl;a
{
}Z$G=;3#
for ( k=ka; k<=kb; k++)
qy (/
{
5i-;bLm
dz
[ja][k] = dz
[ja][k] + .5*hx_inc[ja-1];
8!`.%)- 4
dz
[jb][k] = dz
[jb][k] - .5*hx_inc[jb];
*RE-K36m|u
}
]f @LhC1 x
}
~I ^[rP~
/* Calculate the E from D field */
1,!\7@<CT
for ( i=1; i<IE-1; i++)
Zg f||,
{
+=04X F:
for ( j=1; j<JE-1; j++)
iwx0V
{
W>s9Mp
for ( k=1; k<KE-1; k++)
34M.xB
{
?- & D'
ex
[j][k] = gax
[j][k]*(dx
[j][k]-ix
[j][k]);
?D 9#dGK
ix
[j][k] = ix
[j][k]+gbx
[j][k]*ex
[j][k];
Z/UVKJm>:
ey
[j][k] = gay
[j][k]*(dy
[j][k]-iy
[j][k]);
6uE1&-:L
iy
[j][k] = iy
[j][k]+gby
[j][k]*ey
[j][k];
3EX&.OL!
ez
[j][k] = gaz
[j][k]*(dz
[j][k]-iz
[j][k]);
^* v{ t?u
iz
[j][k] = iz
[j][k]+gbz
[j][k]*ez
[j][k];
0ao Hv
}
N}<U[nh'
}
4 J2F>m40
}
_<}5[(qu
/* Remember: part of the PML is E = 0 at the edges */
gbv[*R{<%
for (j = 0;j < JE-1;j++)
w*kFtNBfU
{
t SLl'XeN
for (k = 0;k < KE-1;k++)
!w\;Q8irN
{
U$J_:~
ez[0][j][k] = 0.0;
W$&Ets8z o
ez[IE-1][j][k] = 0.0;
:q[n1 O[Ch
}
L*oL KigT
}
8VGXw;(Y,d
for (i = 0;i < IE-1;i++)
]=VI"v<X
{
tpzdYokh >
for (k = 0;k < KE;k++)
/ H/Ne )r
{
T.N7`
ez
[0][k] = 0.0;
k3h53QTmC
ez
[JE-1][k] = 0.0;
\@" . GM%
}
q x }fn/:
}
/1 %0A
for (i = 0;i < IE-1;i++)
h#;K9#x6
{
YH tI%
for (j = 0;j < JE-1;j++)
o5@P>\ u>
{
`%I{l
ez
[j][0] = 0.0;
vX24W*7
ez
[j][KE-1] = 0.0;
rX d2[pp
}
/4BXF4ksi,
}
6]49kHgMhe
Z`KXXlJ^i
/*outFile << setw(22) << setprecision(8) << T << setw(22) << setprecision(8) << ez[40][25][25]<< endl;*/
!0+Ex F
g' U^fN
/* Calculate the incident field */
b<cM[GaV~
for ( j=0; j<JE-1; j++)
Rk0 rHC6[
{
vXZz=E AH
hx_inc[j] = hx_inc[j] + .5*(ez_inc[j]-ez_inc[j+1]);
HCy} '}d
}
,qqV11P]
/* Calculate the Hx field */
5Cka."bQ
for ( i=0; i<ia; i++)
U4pvQE.m<
{
3YF]o9
for ( j=0; j<JE-1; j++)
Y b g`Z
{
zQfxw? ~A
for (k=0; k<KE-1; k++)
#e|kA&+8M
{
l:/V%{sx
curl_e = ( ey
[j][k+1] - ey
[j][k] - ez
[j+1][k] + ez
[j][k] );
\uI C<#o"N
ihxl
[j][k] = ihxl
[j][k] + curl_e;
j !n> d
hx
[j][k] = fj3[j]*fk3[k]*hx
[j][k]+fj2[j]*fk2[k]*.5*(curl_e+ fi1
*ihxl
[j][k]);
=W^L8!BE'
}
+~]g&Mf6o
}
{oeQK
}
3<E$m *
for ( i=ia; i<=ib; i++)
<4b o7XH
{
xF31%b`z:
for ( j=0;j<JE-1; j++)
k~F/Ho+R&
{
7B :aJfxM
for ( k=0; k<KE-1; k++)
D 4-U[l+K>
{
oob0^}^
curl_e = ( ey
[j][k+1] - ey
[j][k] - ez
[j+1][k] + ez
[j][k] );
lY?d*qED
hx
[j][k] = fj3[j]*fk3[k]*hx
[j][k]+fj2[j]*fk2[k]*.5*curl_e;
DQQjx>CK
}
~[,TLg 6
}
y_r6T XnGL
}
z:7F5!Z
for ( i=ib+1; i<IE; i++)
dAt[i \S
{
/YR $#&N2
ixh = i-ib-1;
{XW>:EU'N
for ( j=0; j<JE-1; j++)
55KL^+-~
{
j"=jK^
for ( k=0; k<KE-1; k++)
e-t`\5b;
{
9xp ;$ 14
curl_e = ( ey
[j][k+1] - ey
[j][k] - ez
[j+1][k] + ez
[j][k] );
W5g!`f
ihxh[ixh][j][k] = ihxh[ixh][j][k] + curl_e;
\Nyxi7
hx
[j][k] = fj3[j]*fk3[k]*hx
[j][k]+fj2[j]*fk2[k]*.5*(curl_e+ fi1
*ihxh[ixh][j][k]);
_ 9 O'
}
mmK_xu~f28
}
|{"7/~*[
}
]gk1 h=Y~h
/* Incident Hx */
}Pu O$ L
for ( i=ia; i<=ib; i++)
T ua @w+
{
:M`BVZ1t
for ( k=ka; k<=kb; k++)
)j*qGsOg
{
8r ,%! 70
hx
[ja-1][k] = hx
[ja-1][k] + .5*ez_inc[ja];
}H.vH
hx
[jb][k] = hx
[jb][k] - .5*ez_inc[jb];
EHjhe z
}
j(2T,WM
}
"G(/MT^C
e:R[
/* Calculate the Hy field */
AV! cCQ
for ( i=0; i<IE-1; i++)
c:T P7"vG
{
31Du@h8YX
for ( j=0; j<ja; j++)
(k45k/PAP
{
4(IP
for (k=0; k<KE-1; k++)
x{5*%}lX8
{
-lEh}r
curl_e = ( ez[i+1][j][k] - ez
[j][k] - ex
[j][k+1] + ex
[j][k] );
g H.^NO5\'
ihyl
[j][k] = ihyl
[j][k] + curl_e;
SQx):L)P6
hy
[j][k] = fi3
*fk3[k]*hy
[j][k]+fi2
*fk2[k]*.5*(curl_e+ fj1[j]*ihyl
[j][k]);
[i _x 1
}
z<*]h^ !3
}
J n/=v\K@
}
u9(AT>HxT
for ( i=0; i<IE-1; i++)
np(<Ap r
{
2uEu, Y C
for ( j=ja;j<=jb; j++)
G8'3.;"W5
{
mYX) =B{
for ( k=0; k<KE-1; k++)
Dh<e9s:
{
X9wi:
curl_e = ( ez[i+1][j][k] - ez
[j][k] - ex
[j][k+1] + ex
[j][k] );
e)7r
hy
[j][k] = fi3
*fk3[k]*hy
[j][k]+fi2
*fk2[k]*.5*curl_e;
v_ h{_b 8
}
%9M49 s
}
( PlNaasV
}
,Fiiw
for ( i=0; i<IE-1; i++)
+(0eOO'\M
{
2M p;/b!
for ( j=jb+1; j<JE; j++)
{D< ?.'
{
IgJC> ;]u
jyh = j-jb-1;
]Hq%Q~cE
for ( k=0; k<KE-1; k++)
6 [E"
{
y0]O 6.{
curl_e = ( ez[i+1][j][k] - ez
[j][k] - ex
[j][k+1] + ex
[j][k] );
@ RW %EXKt
ihyh
[jyh][k] = ihyh
[jyh][k] + curl_e;
&|eQLY #l
hy
[j][k] = fi3
*fk3[k]*hy
[j][k]+fi2
*fk2[k]*.5*(curl_e+ fj1[j]*ihyh
[jyh][k]);
SF7 Scd
}
HS9U.G >
}
RPwSo.c4
}
[j39A`t7 o
/* Incident Hy */
*^() el,d
for ( j=ja; j<=jb; j++)
_L>n!"E/
{
o~p^`5#
for ( k=ka; k<=kb; k++)
|n8^Xsx4w
{
R qR X
hy[ia-1][j][k] = hy[ia-1][j][k] - .5*ez_inc[j];
L zC~> Uj
hy[ib][j][k] = hy[ib][j][k] + .5*ez_inc[j];
l:HuG!
}
uyIA]OtyN
}
oef(i}8O@
/* Calculate the Hz field */
V- 0Y~T
for ( i=0; i<IE-1; i++)
/8T{bJ5
{
<D}k@M Z
for ( j=0; j<JE-1; j++)
db|$7 ]!w
{
j*}xe'#
for (k=0; k<ka; k++)
3&f{lsLAC
{
uWE@7e4'I
curl_e = ( ex
[j+1][k] - ex
[j][k] - ey[i+1][j][k] + ey
[j][k] );
7dZ!GX?\y
ihzl
[j][k] = ihzl
[j][k] + curl_e;
g;T`~
hz
[j][k] = fi3
*fj3[j]*hz
[j][k]+fi2
*fj2[j]*.5*(curl_e+ fk1[k]*ihzl
[j][k]);
T8LwDqio
}
>{Djx
}
|.5d ^z
}
7 pV 3#fQ
for ( i=0; i<IE-1; i++)
we3t,?`rk7
{
aL}_j#m{
for ( j=0;j<JE-1; j++)
*d 4D9(
{
u QCS%|8C
for ( k=ka; k<=kb; k++)
_nUuiB>
{
4<|]k?@
curl_e = ( ex
[j+1][k] - ex
[j][k] - ey[i+1][j][k] + ey
[j][k] );
fsoS!6h0k
hz
[j][k] = fi3
*fj3[j]*hz
[j][k]+fi2
*fj2[j]*.5*curl_e;
d|3[MnU[a
}
X +R_TC
}
A\ >qoR!Y
}
zF 7T5 Ge
for ( i=0; i<IE-1; i++)
gO%3~f!vY#
{
l]gf T&
for ( j=0; j<JE-1; j++)
%VCHM GP=
{
t)h3G M
for ( k=kb+1; k<KE; k++)
?fG Y,<c
{
zuw6YY8kQ
kzh = k-kb-1;
NG ~sE&,7
curl_e = ( ex
[j+1][k] - ex
[j][k] - ey[i+1][j][k] + ey
[j][k] );
w'C(? ?mH
ihzh
[j][kzh] = ihzh
[j][kzh] + curl_e;
KMa?2cJH#
hz
[j][k] = fi3
*fj3[j]*hz
[j][k]+fi2
*fj2[j]*.5*(curl_e+ fk1[k]*ihzh
[j][kzh]);
