Gregory J. Durnan, PhD, CEng, FIMarEST

jMiniNEC — 1984年のNOSC製アンテナソルバーを1行ずつJavaに移植した、その数学的背景

今回は自分に「とことん深く掘り下げる」と約束していました。このサイトの別の場所には、「昔ながらのMoM(モーメント法)コードがあって、大体こんなことをしている」という、軽い読み物風のバージョンがあります。今回はそれではありません。これはMININECの根底にある実際の数学——具体的には私自身のMININEC、つまり1990年代初頭に私が書いたQB45プログラムを、最近QB64に書き直し、さらに今回、手続き単位・バグ単位で(といっても以前より3つ少ないバグで)Javaのクリーンなヘッドレスクラスへと移植した、jMiniNECと呼んでいるものについてです。

始める前に、範囲についてひとこと補足しておきます。MININECには複雑な系譜があり、このコードが実際に提供する以上のことを約束したくないからです。「フル版」のMININEC(3)——1986年のNOSC技術文書938版、Logan and Rockwayによるもの——は、フレネル反射係数、誘電率、導電率、接地スクリーン、近傍界プローブなどを備えた、損失のある実大地に対応しています。それはこの記事の内容ではありません。私のQB45プログラムはもっと軽量な系統で、1980年のオリジナルのMININEC(NOSC TD 516、Julian, Logan and Rockway)に、完全鏡像大地オプションを加えたものに近いものです。Sommerfeld積分もなく、損失のある大地もなく、接地スクリーンもありません。あるのはワイヤセグメント、細線カーネル、パルス基底関数、点マッチング、そして反転すべき行列だけです。このシンプルさこそが、このコードが数百行のBASICに収まり、この段落よりも少ないメモリしかないマシンで動作した理由です。

このコードは実際どこから来たのか

MININECはサンディエゴにあったNaval Ocean Systems Center(NOSC)で書かれました。これは、現在アメリカ各地に点在するNavy Operational Support Centers(こちらは海軍予備役の管理事務所で、アンテナとは何の関係もありません)と混同しないでください。MININECのNOSCは正真正銘の電磁気学研究所で、1977年3月1日にNaval Electronics Laboratory CenterとNaval Undersea Centerの統合によって設立され、その名前のまま1990年代初頭まで存在し、その後一連の海軍研究所の統廃合を経て、現在のNaval Information Warfare Center Pacificへと引き継がれました。

元となる文献はA. J. Julian、J. C. Logan、J. W. RockwayによるMININEC: A Mini-Numerical Electromagnetics Code、NOSC技術文書516(1980年)です。これはまさに、私のQB45コード自身のバナーに「REFERENCE-JULIAN ET AL, NOSC TD 516」とクレジットされているプログラムであり、40年以上経った今それを確認できたのは嬉しいことです。その発想は、当時としては実に実践的なものでした。NEC(Numerical Electromagnetics Code)は当時も今も、この種の問題に対する本格的で包括的なツールですが、メインフレームが必要でした。MININECはNECの手法を、マイクロコンピュータ上でBASICで、現場で動かせるまでに徹底的に削ぎ落としたものです。「Mini」という名前は伊達ではありません。

物理的背景:ポックリントン方程式、のようなもの

細線モデルを使うMoMコードはどれも同じ出発点を持ちます。マクスウェル方程式、薄い完全導体のワイヤ、そしてワイヤ表面上で接線方向の電界がゼロにならなければならないという境界条件です。これは1897年にポックリントンが導いたもので、MININECも同じ考え方の一形態を使っていますが、ポックリントンの微分されたカーネルを直接使うのではなく、ベクトルポテンシャルとスカラーポテンシャルを介して表現しています——これが計算コストの安さに大きく関わってきます。

ワイヤに沿ったベクトルポテンシャル \(\bar{A}\) とスカラーポテンシャル \(\phi\) は次のように表されます

\[ \bar{A} = \frac{\mu}{4\pi} \int_c I(s')\,\hat{s}(s')\,k(s-s')\,ds', \qquad \phi = \frac{1}{4\pi\varepsilon} \int_c q(s')\,k(s-s')\,ds', \]

ここで \(q(s)\) は線電荷密度で、連続の式 \(q(s) = -\frac{1}{j\omega}\frac{dI}{ds}\) によって電流と関係づけられ、\(k(s-s')\) は細線カーネルです

\[ k(s-s') = \frac{1}{2\pi}\int_{-\pi}^{\pi} \frac{e^{-jkr}}{r}\,d\phi, \qquad r = \left[(s-s')^2 + 4a^2\sin^2\tfrac{\phi}{2}\right]^{1/2}, \]

これは、厳密な自由空間グリーン関数をワイヤの円周(半径 \(a\))方向に平均化することで、問題全体を1次元——ワイヤ軸に沿った弧長 \(s\) の関数——に保ち、完全な2次元の面積分を避けるものです。この平均化こそが「細線近似」であり、MININEC(そしてNECや、この種のあらゆるコード)を扱いやすくしている唯一の仮定です。

これらすべてを結びつける場の方程式は次の通りです

\[ -\bar{E}_{inc}\cdot\hat{s} = -j\omega\bar{A}\cdot\hat{s} - \hat{s}\cdot\nabla\phi , \]

これが表しているのは、ワイヤに沿った入射接線電界が、接線方向のベクトルポテンシャル項(「誘導性」の部分)とスカラーポテンシャルの勾配(「容量性」の部分)によって釣り合っているということです。ここから先はすべて、この1本の方程式——依然として厳密で連続だが、手計算では興味深いどんな問題も解けない——を、\(N\) 個の未知電流に関する \(N\) 本の線形方程式に変換する話です。

アンテナをパルスへと分割する

MININECが用いる離散化は、考えられる中で最もシンプルな基底関数、すなわちパルスです。ワイヤは \(N\) 個の等しいセグメントに分割され、セグメント \(n\) 上の電流は一定であると仮定されます——セグメントごとに複素数1つです。形式的には次のようになります

\[ I(s) = \sum_{n=1}^{N} I_n P_n(s), \qquad P_n(s) = \begin{cases} 1, & s_{n-1/2} < s < s_{n+1/2} \\ 0, & \text{otherwise} \end{cases} \]

ここで \(s_{n\pm 1/2}\) はセグメントの中点です。テスト(「点マッチング」、あるいはコロケーション)は、場の方程式をあらゆる場所で同時に満たそうとするのではなく、各セグメント自身の中点でのみ強制的に成立させます。これはガラーキン法に近い近道——展開に使ったのと同じパルス関数でテストする——であり、まさにこれによって、行列の各要素を埋めるのに二重積分ではなく1回の積分で済むという、行列方程式の安価な構築が可能になります。

ただし、電荷の表現の仕方には、電流のことだけを考えていては見えてこない微妙な点があります。\(q(s)\) は \(I(s)\) の導関数であり、\(I(s)\) はパルスの階段状関数なので、その導関数はセグメント境界でのインパルス列になります。そこでMININECは電荷を、電流パルスと同じ幅だが半セグメント分ずらしたもう1組別のパルス群として表現します。電流パルスはセグメントの接合部に、電荷パルスはセグメントの中心に位置します。この2つのパルス群の間にある半セグメント分のずれこそが、以下のカーネルコードの多くが「P1 ± 0.5」のようなセグメント演算で書かれている理由です。これは単なる飾りではなく、電荷の表現がインデックスの取り方に染み出したものなのです。

この方式では、ワイヤの端点や、ワイヤ同士の接合部の扱いも必要です。自由端を持つワイヤでは、そこのパルスを単に省略します——これは振幅ゼロの半パルスを強制するのと同じことで、未接続の先端で電流をゼロにするという、物理的要請通りの結果になります。2本以上のワイヤの接合部は、パルスを重ね合わせることで処理されます。一方のワイヤの最後のフルパルスが次のワイヤにまたがることを許し、その重なったパルスの振幅が両方のワイヤの方程式で共有されることで、キルヒホッフの電流則が自動的に満たされます。明示的な連続性の制約式を一切書かずに済ませる、巧妙な工夫です。

この仕組みがJava移植版でどう生き残っているかを、defineWire()から直接見てみましょう——ワイヤの方向余弦を計算し、そのパルスが共有座標配列のどこに位置するかを求めています:

// direction cosines (unit vector along the wire) and per-segment length
x3 = x2 - x1;  y3 = y2 - y1;  z3 = z2 - z1;
d = Math.sqrt(x3 * x3 + y3 * y3 + z3 * z3);
dArr[wireIndex] = x3 / d;   // D(I): x-direction cosine
fArr[wireIndex] = y3 / d;   // F(I): y-direction cosine
gArr[wireIndex] = z3 / d;   // G(I): z-direction cosine
sArr[wireIndex] = d / s1;   // S(I): length of one segment

// how many NEW pulses this wire contributes, minus one at each free end
n1 = n + 1;
n = n1 + s1;
if (i1 == 0) n = n - 1;   // end 1 unconnected -> one fewer pulse
if (i2 == 0) n = n - 1;   // end 2 unconnected -> one fewer pulse

// C(*,1) / C(*,2): which wire each pulse's two half-segments belong to
for (j = n1; j <= n; j++) {
    cArr[ix(j)][1] = wireIndex;
    cArr[ix(j)][2] = wireIndex;
}
cArr[ix(n1)][1] = i1;   // overwritten with the connecting wire's number
cArr[ix(n)][2] = i2;    // (or left as +-wireIndex if this end is free/grounded)

このcArr(C())配列は、目立たないところで多くの仕事をしています。行列を埋める段階になる頃には、各パルスは両側それぞれどの実際のワイヤに属しているかを把握しており、これはまさに、別途制約方程式を解くことなく接合部での電流連続性を強制するために、上述の重ね合わせ方式が必要とする情報そのものです。

カーネルのちょっとした「裏事情」

ここが、MININECを単なる「小型版NEC」ではなくMININECたらしめている部分です。インピーダンス行列を埋めるということは、各観測点について、各ソースセグメント上でカーネル \(k(s-s')\) の積分を評価するということです:

\[ \int_{s_u}^{s_v} k(s-s')\,ds'. \]

問題は、\(k\) がまさに一番必要な場所——観測点がソースセグメントのちょうど真上にある自己項(\(r \to 0\))——で対数特異点を持つことです。ここは単純にガウス求積で済ませるわけにはいきません。サンプリングしているまさにその点で被積分関数が発散してしまうからです。

MININECのトリックは、特異点を解析的に差し引いて取り除き、残りの素性の良い部分を数値的に積分し、差し引いた解析的な部分を閉形式で後から足し戻すというものです。カーネルを3つの部分に分割します:

\[ k(s_m - s') = -\frac{1}{\pi a}\ln\!\left[\frac{|s_m - s'|}{8a}\right] \;+\; \frac{\beta\,F(\tfrac{\pi}{2},\beta) + \ln\!\left[\frac{|s_m-s'|}{8a}\right]}{\pi a} \;+\; \frac{1}{2\pi}\int_{-\pi}^{\pi} \frac{e^{-jkr}-1}{r}\,d\phi , \]

ここで第1項は素の対数特異点(閉形式で積分)、第2項は厳密な近傍界の振る舞い——第1種楕円積分 \(F(\tfrac{\pi}{2},\beta)\)——とその対数特異点との差(なめらかで発散せず、数値積分しても安全)、第3項は残り全部を処理する項(こちらもなめらか)です。この3つを足し合わせると、特異点は近似で誤魔化すことなくきれいに消え去ります。

1つだけ厄介なのは、\(F(\tfrac{\pi}{2},\beta)\) 自体が楕円積分であり、1980年代のBASICインタプリタの中で、行列の要素ごとに特殊関数をゼロから評価したい人などいないということです。そこでMININECは、よく知られた有理式+対数の近似式に頼っています:

\[ F\!\left(\tfrac{\pi}{2},\beta\right) \approx \left(a_0 + a_1 m + a_2 m^2 + a_3 m^3\right) \;-\; \left(b_0 + b_1 m + b_2 m^2 + b_3 m^3\right)\ln\!\left(\tfrac{1}{m}\right), \]

ここで \(m = 1 - \beta^2\) であり、\(a\)、\(b\) それぞれに5個ずつの固定係数があります。これは教科書的な近似式(Abramowitz & Stegun流)で、MININECはこの10個の定数を単にBASICソース中のDATAとして持っています。これらはJava移植版のwArrフィールドに収まっており、gosub9()(BASICのTARGET9サブルーチンに対応する移植版の名前)にある近傍界補正の仕組みは、すべてここから来ています:

if (i6 != h0) {                                 // "near" case only
    d3 = pointDist * pointDist - a2;            // (a2 = wire-radius^2)
    b = d3 / (d3 + 4 * a2);                     // beta^2, roughly
    v1 = wArr[1]; v2 = wArr[6];                 // a0, b0
    b1 = 1; b2 = -36;
    for (l2 = 2; l2 <= 5; l2++) {                // accumulate a1..a4, b1..b4
        b2 = b2 / (l2 - 1);
        if (b < Math.pow(10, b2)) break;         // early-out once terms are negligible
        b1 = b1 * b;
        v1 = v1 + wArr[ix(l2)] * b1;
        v2 = v2 + wArr[ix(l2 + 5)] * b1;
    }
    v = (v1 - v2 * Math.log(b)) * Math.sqrt(1 - b);   // the F(pi/2, beta) approximation
    d3 = Math.sqrt(d3) / 8 / aArr[ix(p4)];
    t3 = t3 + (Math.log(d3) + v) / p / aArr[ix(p4)] - 1 / pointDist;
}

行列を埋める

カーネル積分が求まれば、インピーダンス行列の各要素は次のようになります

\[ Z_{mn} = \frac{-1}{4\pi j\omega\varepsilon} \Bigg[ k^2(\bar{r}_{m+1/2} - \bar{r}_{m-1/2})\cdot (\hat{s}_{n+1/2}\,\psi_{m,n,n+1/2} + \hat{s}_{n-1/2}\,\psi_{m,n-1/2,n}) - \cdots \Bigg], \]

これは、スカラーポテンシャル積分 \(\psi\)(前節のカーネル積分)と、観測セグメントとソースセグメントの間の方向余弦の内積を組み合わせたものです。平たく言えば、ソースセグメント \(n\) に沿って流れる電流が、遅延近傍界/遠方界カーネルを介して、MININECがセグメント \(m\) でテストしている接線方向の場へどれだけ強く結合するかということです。

この方向余弦の内積は行列を埋めるコードの中に明示的に現れており、上記のベクトル代数をそのまま直訳したものになっているので、見ておく価値があります:

// T5,T6,T7 = the observation segment's effective direction vector,
// combining the two half-segments either side of pulse I's junction
t5 = f5 * dArr[ix(i2)] + f4 * dArr[ix(i1)];
t6 = f5 * fArr[ix(i2)] + f4 * fArr[ix(i1)];
t7 = f5 * gArr[ix(i2)] + f4 * gArr[ix(i1)];

大地の扱いこそ、MININECのシンプルさが本当に報われる部分です。完全導体の接地面に対しては、鏡像理論により、大地をワイヤ構造全体の鏡像コピーで置き換え、既知の符号で関連づけられた電流を持つ、2倍の未知数を解くことになります。素朴に実装すれば、本当に行列サイズが2倍になってしまうところです。ところがMININECの行列構築ループは、内側にK = 1 To Gループ(自由空間ではGは1、「大地近傍」では2)を回し、実際に \(2N \times 2N\) の系を組み立てて縮約するのではなく、鏡像による寄与を、鏡映幾何に応じて \(z\) 成分の符号を反転させたうえで、既存の \(N\times N\) 行列要素へ直接足し込みます。同じ物理を半分のメモリで実現しており、これこそ、接地の有無が誘電率や導電率ではなく単一のBASIC Input文で済んでいる理由です。このコードパスのどこにも損失のある大地の項は存在せず、あるのは鏡像だけです。

解く:複素数演算によるガウスの消去法

\([Z_{mn}][I_n] = [V_m]\) が埋まったら、誰かがそれを反転しなければなりません。このコードがやっているのは、いかにも1984年らしい直接的な方法です。当時のBASICにはネイティブの複素数型がなかったため、実部・虚部それぞれの並行した行列(オリジナルではR()/B()、ここではrArr/bArr)に対して、ピボット選択付きのその場(in-place)ガウス消去法を行っています。

ここで指摘しておく価値があるのは、ピボット探索が単なる部分ピボット選択(現在の列の中で最大の要素)ではなく、行・列の両方について、まだ消去されていない残りの部分行列全体から最大絶対値の要素を探し、それを所定の位置へ入れ替えている点です:

for (int jj = ix(h0); jj <= nn; jj++) {
    if (pArr[jj] == h0) continue;                 // already used as a pivot row
    for (int kk = ix(h0); kk <= nn; kk++) {
        if (pArr[kk] == h0) continue;              // already used as a pivot column
        s2 = rArr[jj][kk] * rArr[jj][kk] + bArr[jj][kk] * bArr[jj][kk]; // |Z|^2
        if ((s1 - s2) > 0) continue;
        i1 = jj; i2 = kk; s1 = s2;                 // new best pivot candidate
    }
}

これは、多くの教科書が教える行のみの部分ピボット選択よりも、むしろ完全ピボット選択に近いものです。借り物の計算機時間で無人稼働することを前提としたコードとしては、なかなか堅実な数値的配慮と言えます。条件の悪い行列が黙ってでたらめな結果を返すことの方が、消去のステップごとに余分な \(O(N^2)\) の探索コストがかかることよりもはるかに悪いからです。ピボットが見つかり入れ替えられた後の消去処理自体は、実数版アルゴリズムを複素数へと拡張した標準的なもので、ピボット行をピボット要素で割り(ここでもネイティブの複素数型がないため、複素数の逆数を手計算で \((a - bi)/(a^2+b^2)\) として求めます)、それを他のすべての行から適切な倍数だけ差し引きます。P()(ここではpArr)はどの行・列がすでにピボットとして使われたかを記録し、H()/I()はピボットの履歴を記録することで、消去が終わった後に解を元の行の順序へ戻せるようにしています。

放射パターンを求める

最後の段階では、求まった電流分布を遠方界の放射パターンへと変換します。パルス電流の集合 \(I_n\) に対し、方向 \((\theta,\phi)\) における遠方電界は(方向に依存しない定数を除けば)次のようになります

\[ \bar{E}(\bar{r}_0) \propto \sum_{n} I_n\, \big[(\hat{k}\cdot\hat{s}_n)\hat{k} - \hat{s}_n\big]\, e^{\,j\bar{K}\cdot\bar{r}_n}, \]

これは、すべての電流パルスについて、そのパルスの寄与を観測方向 \(\hat{k}\) に対する位相で重み付けして足し合わせたものです。各ワイヤセグメントの寄与は、振幅項(有限長のセグメント上で一定電流を積分した、本質的にsinc関数状のもの)と、\(\bar{K}\cdot\bar{r}_n\) による進行位相項に分解でき、この両方が移植版でもそのまま残っています:

// one pulse's contribution to the far-field sum
o = -(r1 * dArr[lIdx] + r2 * fArr[lIdx] + r3 * gArr[lIdx]) * f3v;
s1 = w * sArr[lIdx] / 2;
if (Math.abs(o) < .0000001) {
    o = (2 - o * o * s1 * s1 / 3) * s1;   // small-angle series: avoids 0/0
} else {
    o = 2 * Math.sin(o * s1) / o;          // the sinc-shaped element pattern
}
s2 = w * (xArr[pulseCoord] * r1 + yArr[pulseCoord] * r2 + zArr[pulseCoord] * r3);
s1 = Math.cos(s2);  s2 = Math.sin(s2);      // progressive phase, this pulse's position
b1 = (s1 * hArr[pulseI] - s2 * iArr[pulseI]) * o / f4v;   // complex multiply:
b2 = (s1 * iArr[pulseI] + s2 * hArr[pulseI]) * o / f4v;   // (current) x (phase) x (element pattern)

インピーダンス行列のときと同様、接地面がある場合も配列サイズを2倍にするのではなく、総和の項を2倍にします(直接波に加えて鏡像波)。ここでもやはり、完全導体による鏡像理論以外の何物でもありません。損失のある実際の土壌からのフレネル反射が欲しいという人は、それはこのコードではなくMININEC(3)の領域です。実は以前、この移植版を、元のソースに付属していた古いヘルプファイルと突き合わせていたときに、まさにこの違いに行き当たりました。そのヘルプファイルの計算例では実大地(誘電率13、導電率0.005 S/m)を使っており、これはこのプログラムが決して尋ねることのない値です。同じ形状をjMiniNECに完全鏡像大地オプションで通したところ、抵抗値はヘルプファイルの答えから0.2%以内、リアクタンスは約6%以内に収まりました。これはまさに予想通りの一致度です。抵抗値は近く、リアクタンスはややずれる——完全導体と損失のある実大地は似てはいても同一の接地ではないからです。バグではなく、良好な健全性チェックの結果です。

実際に動くのか?

アンテナ理論において口先だけの話は安いものなので、実際に動いている移植版をお見せします。30MHzの半波長ダイポール、10セグメント、中央のパルスに1Vで給電した場合です:

これはまさに期待通りの電流分布であり、この移植作業の途中ではしばらくの間、まったく得られていなかったものです——給電点でなめらかなピークを持ち、両端の自由端に向かって対称的に(ほぼ)ゼロまでテーパーしていく分布です。ここに至るまでには、単なる移植ミスではなく、1984年のBASICコード自体に潜んでいた3つの別々のバグを見つけて修正する必要がありました:

• 近傍界カーネルに紛れ込んでいた余分な係数で、座標平面上に寝ているワイヤがあるとNaNになってしまうもの、
• 遠方界の位相項で、配列参照であるべきところに素のスカラーが使われていたもの、
• そして最も見つけにくかったもの:1つの配列スロットが、ワイヤの恒久的な方向余弦とカーネルの一時的な距離計算用の作業領域という二重の役割を担っていて、計算の途中で自分自身を静かに上書きしてしまうもの。

この3つはすべて、移植版自身のコメント内に、元のBASICコードの正確な行番号付きで記録してあります。QB64版でも心当たりがある方のために。

上の放射パターンは教科書通りの半波長ダイポールの8の字形で、正面方向(broadside)利得はおよそ2.1dBiです(理想的な無限に細い半波長ダイポールの教科書値は2.14dBiで、このモデルが1mmという有限の半径でセグメント数もわずか10であることを考えれば十分近い値です)。給電点インピーダンスは \(67.71 - j17.35\,\Omega\) となり、この半径・セグメント数では真の半波長共振よりわずかに短い、というダイポールらしい妥当な範囲に収まっています。

最終的な答えだけでなく、生のインピーダンス行列——実際の \(Z_{mn}\) の各要素——についても、MININECの独立した、現在も活発にメンテナンスされているPython再実装であるpymininecと突き合わせました。7MHzにおける21.4mのワイヤ(そのプロジェクト自身のテストスイートからそのまま持ってきたテストケース)について、jMiniNECの自己項は \(-8.4796 - j0.009668\) となり、参照値の \(-8.4851 - j0.009669\) と一致、最初の非対角項も \(4.2375 - j0.009573\) 対 \(4.2404 - j0.009573\) で、有効数字3〜4桁の一致が見られ、行全体も素性の良い細線カーネルらしく、なめらかかつ単調に減衰していきました。1980年代のBASICプログラムとして、これ以上ないくらい「正しさが証明できた」に近い状態と言えるでしょう。

まとめ

以上が数学的な全体像です。細線積分方程式を、パルス基底関数と点マッチングで手なずけ、1984年からずっと10個のDATA定数の中に眠っていた楕円積分近似でカーネルの特異点を差し引き、行列サイズを2倍にする代わりに接地面を総和の2倍化として利用して行列を埋め、驚くほど注意深い複素ガウス消去法で解き、閉形式の遠方界総和によって放射パターンへと変換する——それだけです。しかもそのすべてが、今日このページを表示するのにも苦労するようなマシンの上で、かつてゆっくりと動いていたほど小さなものでした。

完全なヘッドレスJavaエンジン(入出力なし、組み込み利用を想定したもの)と、それを動かすコンソールドライバは、ご希望があればお渡しします——このページはコードの垂れ流しではなくアイデアについてのページにしたいので、全文を載せるには長すぎます。40年前の海軍製BASICを、バグも含めて忠実に移植し、その上で明確にバグだと言えるものだけを直すとどうなるのか——それを見てみたい方には、まさにそういう話です。

#GregoryJDurnan #DurnanOrg #GregDurnanRF #RFEngineering #MININEC #MethodOfMoments #AntennaTheory #NumericalElectromagnetics #AntennaDesign #ElectromagneticsResearch #NavalHistory