Ｅ波問題（線電流源）

以下に示すフォートラン言語によるプログラムでは線電流源を 振幅１として座標の原点に置いている。 また、折れ線近似した散乱体（完全導体）の形状をデータとして与えて、 各折れ線の中点での電流密度と放射パターン（遠方放射界）を計算している。

プログラムでは、 $X_n=J_n\Delta_n$を未知数として、 連立方程式、

$$\sum_{n=1}^{N} G_{m,n} X_n = F_m$$

を解いている。 また、寸法はすべて波長で規格化し、

$$\begin{eqnarray*} G_{m,n} &=& \left\{ \begin{array}{ll} H_{0}^{(2)}(2\pi R_{... ... &n=m \end{array}\right. \\ F_m &=& -H_{0}^{(2)}(2\pi r_{m}) \end{eqnarray*}$$

により計算している。

このとき、G(m,n)は対称行列となる。

8行のサブルーチン struct の中で折れ線近似した散乱体の形状を読み込み、 各折れ線の中点を要素点として、その座標$x_n$, $y_n$)を 配列要素elem(1,n), elem(2,n)に代入し、 各要素幅をelem(3,n)に代入する。 要素点数は変数nsizeに返される。

9-25行で連立方程式の係数$G_(n,m)$を配列要素G(m,n)に格納している。 $G_{n,m}$が対称行列であることを利用して計算を進めている。 また、連立方程式の右辺式$F_m$は配列要素G(m,nsize+1)に格納している。

27行でsubroutine cswepmの内部でガウスの消去法を用いて連立方程式を 解いている。 方程式の解は配列要素cg(n,nsize), (n=1,2,$\cdots$,n)に返される。 28-31行には連立方程式を解くことができなかった場合の処理を記述している。

33-45行で電流分布を計算し、ファイルi.datに書き出している。

41 write(10,1000) elem(1,n),elem(2,n),t,aj,ang

により、i.datの各行には、要素点の$x$,$y$座標、 散乱体表面に沿って取った座標(t)、 電流密度の絶対値(aj[A/波長])、 その位相角(ang)を書き出している。

47-63行で放射パターンを計算している。 求めた放射パターンはファイルpattern.datに書き出される。

48 do 50 akaku = 0., 360., 0.25

により、$0$から$360{\circ}$まで$0.25^{\circ}$きざみで計算される。 pattern.datの各行には、観測角akakuと dBを単位とした放射界の強さaeが書き出される。

67-76行の関数dbでは、 複素数で与えられた電界よりdBを単位とした電界の強さを計算している。

78-96行のサブルーチンstructは散乱体形状を読み込んで 要素点の座標などを算出している。

プログラム bem.f

1 implicit complex (c)
2 parameter(nmax=128, gamma=0.57721 56649 01532 86)
3 complex cg(nmax,nmax+1)
4 real elem(3,nmax)
5
6 pi = acos(-1.0)
7 pi2 = 2*pi
8 call struct(elem, nmax, nsize)
9 c G(n,m)
10 do 20 m = 1, nsize
11 do 10 n = 1, m-1
12 10 cg(n,m) = cg(m,n)
13 cg(m,m) = cmplx(1., -2./pi*
14 + (gamma-1.+alog(0.25*pi2*elem(3,m))))
15 do 12 n = m+1, nsize
16 rk = pi2*sqrt((elem(1,n)-elem(1,m))**2
17 + +(elem(2,n)-elem(2,m))**2)
18 cg(n,m) = cmplx(ajbesr(0.,rk),-ybesr(0.,rk))
19 12 continue
20 20 continue
21 c F(m)
22 do 22 m = 1, nsize
23 rk = pi2*sqrt(elem(1,m)**2 + elem(2,m)**2)
24 cg(m,nsize+1) = -cmplx(ajbesr(0.,rk),-ybesr(0.,rk))
25 22 continue
26
27 call cswepm(cg,nmax,nsize,1,1.e-5,ier)
28 if(ier.ne.0) then
29 write(*,*) 'ier=', ier
30 stop
31 endif
32
33 open(unit=10,file='i.dat',form='formatted')
34 t = 0
35 do 30 n = 1, nsize
36 h = 0.5*elem(3,1)
37 t = t + h
38 aj = abs(cg(n,nsize+1))/elem(3,n)
39 ang = 180./pi *
40 + atan2(aimag(cg(n,nsize+1)),real(cg(n,nsize+1)))
41 write(10,1000) elem(1,n),elem(2,n),t,aj,ang
42 1000 format(3f10.6,e15.6,f8.2)
43 t = t + h
44 30 continue
45 close(10)
46
47 open(unit=10,file='pattern.dat',form='formatted')
48 do 50 akaku = 0., 360., 0.25
49 phi = akaku/180. * pi
50 ce = 0.0
51 do 40 n = 1, nsize
52 r = sqrt(elem(1,n)**2 + elem(2,n)**2)
53 p = atan2(elem(2,n), elem(1,n)) - phi
54 a = pi2 * r * cos(p)
55 ce = ce + cg(n,nsize+1)*
56 + cmplx(cos(a),sin(a))
57 40 continue
58 ce = 1 + ce
59 ae = db(ce)
60 write(10,1100) akaku, ae
61 1100 format(2f20.6)
62 50 continue
63 close(10)
64 stop
65 end
66
67 function db(c)
68 complex c
69 a = abs(c)
70 if(c .eq. 0.) then
71 db = -200.
72 else
73 db = 20.*alog10(a)
74 endif
75 return
76 end
77
78 subroutine struct(elem, nmax, nsize)
79 real elem(3,nmax)
80 read(*,*) x,y
81 n = 0
82 2 continue
83 x0 = x
84 y0 = y
85 read(*,*,end=100) x,y
86 n = n + 1
87 elem(1,n) = 0.5*(x+x0)
88 elem(2,n) = 0.5*(y+y0)
89 elem(3,n) = sqrt((x-x0)**2 + (y-y0)**2)
90 goto 2
91
92 100 continue
93 write(*,*) '要素数=', n
94 nsize = n
95 return
96 end