子午線弧長とは何か — 投影計算の心臓部

前回は、楕円体のあらゆる定数が長半径a・扁平率fの2値から導けるという話を書きました。数式連載の第1回となる今回は、そのa・fを使って実際に計算される最初の重要な値——子午線弧長——を扱います。平面直角座標系のX座標の正体そのものです。

なぜ弧長が要るか — X座標の正体

平面直角座標系のX（南北方向の座標値）は、突き詰めると「赤道からその緯度までの、子午線に沿った距離」を、系ごとの原点からのオフセットで補正したものです。つまり緯度から平面座標を作る計算は、必ず「この緯度は赤道から何メートル離れているか」という子午線弧長の計算を土台にしています。

座標変換の流れ

ガウス・クリューゲル投影の式全体は複雑に見えますが、その中心にあるのは「緯度→子午線弧長」という、比較的シンプルな1変数の計算です。

円なら簡単・楕円だと積分になる

もし地球が完全な球であれば、弧長は「半径×角度（ラジアン）」で一発で求まります。しかし楕円体では、前回記事で紹介した子午線曲率半径M(φ)が緯度φによって変化するため、話はそう単純ではありません。

子午線弧長S(φ)は、赤道（φ=0）からその緯度φまで、M(φ)を緯度方向に足し合わせた（積分した）値になります。

S(φ) = ∫[0→φ] M(φ) dφ

問題は、この積分が初等関数（三角関数や多項式の組み合わせ）だけでは閉じた形に書けないことです。これは楕円積分と呼ばれるタイプの積分で、厳密解を代数的に表現できません。

前回記事：楕円体は a と f だけで決まる — 離心率・曲率半径の導出

級数展開の考え方 — 近似でどこまで迫れるか

閉じた形で解けないなら、実用上は「級数展開」で近似します。M(φ)を離心率の2乗e²を小さなパラメータとしてべき級数に展開し、それぞれの項を緯度で積分すると、sin2φ・sin4φ・sin6φ…といった多重角のサイン項の和として弧長を表現できることが分かります。

S(φ) ≈ a・(A₀φ − A₂sin2φ + A₄sin4φ − A₆sin6φ + …)

（A₀・A₂・A₄…はa・fから導かれる係数で、地理院の公式資料や測地学の教科書に具体的な導出が示されています。本記事では構造の理解を優先し、係数の桁数まで厳密に転記することはしません。）

級数は項数を増やすほど精度が上がりますが、GRS80のように扁平率が非常に小さい楕円体では、収束が速く、実務で使う程度の項数（4〜5項程度）でミリメートル精度に達すると言われています。「無限に続く積分」を「有限個のサイン項の和」に置き換えられる、というのがこの手法の核心です。

歴史コラム — メートルの定義と子午線測量

子午線弧長は、実は「メートル」という単位そのものの起源とも深く関わっています。1791年、フランスはパリを通る子午線の、北極から赤道までの弧長の1000万分の1を1メートルと定義しました。定義のために、実際にダンケルクからバルセロナまでの子午線弧長を測量する大事業が行われています。

日本でも、江戸時代に伊能忠敬が全国測量の過程で「緯度1度の距離」を実測しています。当時の測量技術でどこまで正確に子午線弧長に迫れたのか、という話は、それだけで1本の記事になるテーマなので、詳しくは姉妹サイトGeoPrism JPの学習記事に譲ります。

GeoCoreJP実装との対応

GeoCoreJPでは、前回記事で紹介したEllipsoidConstants（a・fおよびそこから導かれるe²）を受け取り、緯度から子午線弧長を計算する関数として実装しています。イメージは次のような構造です。

extension EllipsoidConstants {
    // 緯度latRad（ラジアン）における子午線弧長（メートル）
    func meridianArcLength(latRad: Double) -> Double {
        // e2, e4, e6 ... はe2から順次計算するべき級数の係数
        let e2 = self.e2
        let e4 = e2 * e2
        let e6 = e4 * e2

        // A0, A2, A4, A6 は a・e2 から導かれる係数（式は教科書・地理院資料を参照）
        let a0 = 1 - e2/4 - 3*e4/64 - 5*e6/256
        let a2 = 3*(e2 + e4/4 + 15*e6/128) / 8
        let a4 = 15*(e4 + 3*e6/4) / 256
        let a6 = 35*e6 / 3072

        return a * (a0 * latRad
                     - a2 * sin(2 * latRad)
                     + a4 * sin(4 * latRad)
                     - a6 * sin(6 * latRad))
    }
}

係数の詳細な導出過程はここでは省略しますが、実装のポイントは「a・fという2つの根本定数だけから、係数A₀〜A₆をすべて計算で導ける」という点です。前回記事の設計方針（根本の2値だけを信頼できる定数として持つ）が、ここでもそのまま活きています。