リスト3-6-1が遠方界を計算する関数です。
予め境界面の電磁界(複素数)を計算しておきます。複素数の演算は関数を通して行います。
ソースコードの説明
20-39行:r,θ,φ方向の単位ベクトルを計算します。
51行:境界面のデータは構造体配列surfaceに入っています。
53-64行:境界面の電界と磁界から電流と磁流を計算します。
70-78行:ポテンシャルN,LのX,Y,Z成分を計算します。
81-87行:座標変換によりポテンシャルN,Lのθ,φ成分を計算します。
リスト3-6-1 遠方界の計算(farfield.c)
1 /*
2 farfield.c
3
4 far field
5 */
6
7 #include "ofd.h"
8
9 typedef struct {
10 d_complex_t *ex, *ey, *ez;
11 d_complex_t *hx, *hy, *hz;
12 double nx, ny, nz;
13 double px, py, pz;
14 double ds;
15 } surface_t;
16 surface_t *surface;
17
18 void calc_farfield(double theta, double phi, int ifreq, double kwave, d_complex_t *etheta, d_complex_t *ephi)
19 {
20 // unit vector in r, theta, phi
21
22 double r1[3], t1[3], p1[3];
23
24 double cos_t = cos(theta * DTOR);
25 double sin_t = sin(theta * DTOR);
26 double cos_p = cos(phi * DTOR);
27 double sin_p = sin(phi * DTOR);
28
29 r1[0] = +sin_t * cos_p;
30 r1[1] = +sin_t * sin_p;
31 r1[2] = +cos_t;
32
33 t1[0] = +cos_t * cos_p;
34 t1[1] = +cos_t * sin_p;
35 t1[2] = -sin_t;
36
37 p1[0] = -sin_p;
38 p1[1] = +cos_p;
39 p1[2] = 0;
40
41 d_complex_t pnx = d_complex(0, 0);
42 d_complex_t pny = d_complex(0, 0);
43 d_complex_t pnz = d_complex(0, 0);
44 d_complex_t plx = d_complex(0, 0);
45 d_complex_t ply = d_complex(0, 0);
46 d_complex_t plz = d_complex(0, 0);
47
48 int64_t nsurface = 2 * ((Nx * Ny) + (Ny * Nz) + (Nz * Nx));
49
50 for (int n = 0; n < nsurface; n++) {
51 surface_t *p = &surface[n];
52
53 // ZJ = n X ZH
54 d_complex_t cjx = d_sub(d_rmul(p->ny, p->hz[ifreq]), d_rmul(p->nz, p->hy[ifreq]));
55 d_complex_t cjy = d_sub(d_rmul(p->nz, p->hx[ifreq]), d_rmul(p->nx, p->hz[ifreq]));
56 d_complex_t cjz = d_sub(d_rmul(p->nx, p->hy[ifreq]), d_rmul(p->ny, p->hx[ifreq]));
57
58 // M = -n X E
59 d_complex_t cmx = d_sub(d_rmul(p->ny, p->ez[ifreq]), d_rmul(p->nz, p->ey[ifreq]));
60 d_complex_t cmy = d_sub(d_rmul(p->nz, p->ex[ifreq]), d_rmul(p->nx, p->ez[ifreq]));
61 d_complex_t cmz = d_sub(d_rmul(p->nx, p->ey[ifreq]), d_rmul(p->ny, p->ex[ifreq]));
62 cmx = d_rmul(-1, cmx);
63 cmy = d_rmul(-1, cmy);
64 cmz = d_rmul(-1, cmz);
65
66 // exp(jkr * r) * dS
67 double rr = (r1[0] * p->px) + (r1[1] * p->py) + (r1[2] * p->pz);
68 d_complex_t expds = d_rmul(p->ds, d_exp(kwave * rr));
69
70 // ZN += ZJ * exp(jkr * r) * dS
71 pnx = d_add(pnx, d_mul(cjx, expds));
72 pny = d_add(pny, d_mul(cjy, expds));
73 pnz = d_add(pnz, d_mul(cjz, expds));
74
75 // L += M * exp(jkr * r) * dS
76 plx = d_add(plx, d_mul(cmx, expds));
77 ply = d_add(ply, d_mul(cmy, expds));
78 plz = d_add(plz, d_mul(cmz, expds));
79 }
80
81 // ZN-theta, ZN-phi
82 d_complex_t pnt = d_add3(d_rmul(t1[0], pnx), d_rmul(t1[1], pny), d_rmul(t1[2], pnz));
83 d_complex_t pnp = d_add3(d_rmul(p1[0], pnx), d_rmul(p1[1], pny), d_rmul(p1[2], pnz));
84
85 // L-theta, L-phi
86 d_complex_t plt = d_add3(d_rmul(t1[0], plx), d_rmul(t1[1], ply), d_rmul(t1[2], plz));
87 d_complex_t plp = d_add3(d_rmul(p1[0], plx), d_rmul(p1[1], ply), d_rmul(p1[2], plz));
88
89 // F-theta, F-phi
90 *etheta = d_add(pnt, plp);
91 *ephi = d_sub(pnp, plt);
92 }