3.6 遠方界

リスト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	}