Vincenty's Formulae

はじめに

過去に仕事で必要になったので、大圏航路（大圏コース）にかかる計算法などを調べ、ソフトウェアとして実装した。
それまでは、理論すら知らなかった。

ここでは自身の経験に基づき、計算法の紹介と、コンピュータプログラムとして実装する際に、 どういったことに気をつけるべきかの私見を残したいと思う。
ネットの時代なので、もっとより良い知見は他にもあるはずなので、ここに書かれたことを鵜呑みにしないでもらいたい。

当然、ここに書かれたことに基づき読者に何らかの損害が発生したとしても、一切の責は負わない。

なお、出てくる Vincety's Formulae とは何かや使われる用語などは、必要に応じてウィキペディア などでより詳しい情報を得られたし。 ただし、2014年9月時点では、Vincenty's Formulaeは日本語版のウィキペディアには存在していない。英語版の Wikipediaを参照されたし。
(2023年6月時点では、日本語版のウィキペディアにも載っていますね。)

大圏航路とは

地球表面上の２地点を結ぶ地球表面をなぞる最短コースのことである。
地球表面をなぞる３次元的には曲線（実は円弧に近い）もので、地中に潜り込んで２地点を結ぶ直線ではないことに注意されたし。

地球のモデル

地球は、半径約6370キロメートル、外周（大円）約40,000キロメートルのほぼ球体をしているが、真球(perfect sphere)ではない。
赤道付近が膨らみ、北極と南極を上下からちょっと押し潰したような形をしている。 （楕円形を短径を軸に三次元空間で回転させたときに描かれる形、即ち回転楕円体(spheroid)に近い。） 特に、「地球楕円体(Earth ellipsoid)」という用語が存在する。

実際の地球には、山あり谷あり、重力場の揺らぎあり、などなどの要因からもろもろの凹凸があるわけであるが、 ここで紹介している大圏航路の計算ではそれらは無視しており、完全なる回転楕円体にモデル化した上での話である。

以下、地球のモデルをいくつか挙げておく。

Vincenty's Formulae

何を解くのか？ というと次の２つの問題を解くものである。

• ある１地点を与える。さらにその地点からの方位角及び距離を与える。そのときの、到達地点を求める。
  これをVincentyのpaperでは、direct問題と称する。
• ある２地点を与える。そのときの２地点間の大圏コースの距離を求める。また、一方の地点から他方の地点への方位角をも求める。
  これをVincentyのpaperでは、inverse問題と称する。

地球が完全な球体であれば、球面三角法という解法が存在するのであるが、あいにく地球は完全な球体ではなく回転楕円体（地球楕円体）である。
回転楕円体においては、球面三角法と計算量が同等の解法は存在しない。

Vincentyの計算法の特徴は次にあげられる。

• 積分（微小区間に区切って計算したものの足し算）をするのではない。
  ある計算式にある初期値を与えて計算した結果を次の初期値として計算するといった反復（iteration）を行い、計算結果がある値に収束したら、 収束した値を解とするという計算法である。（解の公式を持たない方程式の解を計算機のパワーを利用して求めるニュートン法のような解法である、 といったら「反復」と「収束」の計算方法のイメージが湧くであろうか。）
• 反復回数が少ないため、トータルの計算量及び計算時間が少なくてすむ。
• 反復させる１回の計算式の数も少ない。プログラミング言語的に言えば、ループの中の処理も少ないということだ。
  そもそも、Vincentyは1975年にpaperを書いている。当時の計算機のパワーは現代のそれに比べてはるかに小さいものであった。 記憶できるプログラムの量も極めて小さいという制限があるなかで、いかに有効に計算できるかを追及した賜物である。
• 解の精度がとても高い。
  direct問題で、18000km先の地点の解の誤差は、1mm（ミリメートル）もない。
  地球の裏側へ到達するような計算をしたときにでも到達地点の誤差が1mmもないということである。 ビックリするとともに、この精度が実用上必要になる局面は思いつかない。 だって、先に述べたように実際の地球は完全なる回転楕円体でもないのだから。

ただし、inverse問題の計算において、地球のほぼ裏側同士に位置する2地点を与えたときには、反復を繰り返しても解が収束しないという欠点がある。

計算には計算誤差対策を打とう

（計算機は無限の桁の数値を扱うことはできない。そのため、計算機とプログラムでの実装では必ず計算誤差が生じる。
どのような計算誤差を考察し、プログラミングではどういった工夫をすべきか、について語りたい。）

double型変数の特徴を理解しておこう

（コンピュータのプログラミングで、実数を扱うときに整数型変数ではなくて、float型（単精度浮動小数点数）やdouble型（倍精度浮動小数点数）を使う。
この型はどのようにコンピュータの中で表現されているのか、それ故にどのような特徴を持つのか、さらにその特徴を理解すると、 プログラミング上で考慮すべき点は何なのかについて語りたい。
だって、「実数を入れる器は、floatかdoubleなんだよ。以上。」では、なんだか寂しいじゃん。それに、思ったようにプログラムは動いてくれないかもよ。）

コンピュータの中では、二進数が使われる。これは、ディジタル回路の究極の最小単位は、二値（電圧がかかっている状態とかかっていない状態）であることに基づく。
二値の一桁で表すことのできる情報の区別は、当然二個であり、この単位をビット（bit）という。
2ビット使うと4個の状態を区別でき、3ビット使うと8個の状態、...8ビット使うと156個の状態を区別して扱うことができる。

正の整数は、必要なビット数の領域に二進数で格納する。1ビットは、0か1の値しか持てず、それを1桁として上位の桁に繰り上げていく。だから二進数なのだ。

同様に、小数も二進数での表現になる。(1/2)の桁、(1/4)の桁、(1/8)の桁、・・・の値を持つかどうかで表現することになる。
例えば、十進小数で 0.1 や 0.9 という値は、(十進で)
0.1₍₁₀₎ = (1/2) x 0 + (1/4) x 0 + (1/8) x 0 + (1/16) x 1 + (1/32) x 1 + (1/64) x 0 + ・・・
0.9₍₁₀₎ = (1/2) x 1 + (1/4) x 1 + (1/8) x 1 + (1/16) x 0 + (1/32) x 0 + (1/64) x1 + ・・・
となるので、二進小数では、
0.1₍₁₀₎ = 0.000110・・・₍₂₎
0.9₍₁₀₎ = 0.111001・・・₍₂₎
といった感じで、有限桁では正確には保持することはできず、近似値を持つことになる。
より詳しくは、「IEEE754-2008」などを参照されたし。

さらに徐々につづく。。。

三角関数 atan2() とは

（現代のコンピュータの世界では、三角関数（sin, cos, tan...）は、ライブラリとして用意されている。
同様に逆関数（逆三角関数）も用意されているのだが、asin（アークサイン）、acos（アークコサイン）、atan（アークタンジェント）の他に、 atan2というアークタンジェントなんだけどatanとはちょっと違うモノもある。
atan2()がどういった局面で有効な関数であるのかを考察してみたい。）

C++/Javaによるサンプルプログラム

（プログラミング言語C++を使って、かつC++風（オブジェクト指向風）に、サンプルプログラムを書いてみたい。
過去に書いたのは会社の資産であり、おいらの資産ではないので、ここに載せることは出来ない。それ以前に会社から持ち出すことも出来ない。 自身のプログラミング技術を錆びつかせないためにも、一から書いてみようと思う。）

2023/6/12 C++言語ではなくて、Java言語での実装例を公開する。
ちなみに、計算ロジックのサンプルコードは2023年時点でネット上のいたるところに存在している。また、はやりのchatGPTだって、コードを書き出してくれる。 計算結果を取得するだけのサイトもたくさんある。国土地理院のサイトでも計算させることができる。
それなのに、ここで公開しようと思うに至った、としあきが重点を置いたものとは・・・

• コードの再利用性・拡張容易性とは具体的な実装においてどうすることなのか
  本プログラムではメインプログラム内に、どの地球楕円体を使い、 どの計算アルゴリズムを使用しているのかのコードを含まない。
  将来新たなアルゴリズムや地球楕円体が登場したときに、その追加のための追加のコードは別途必要になるが、 既存のコードの改修やリコンパイルは不要、即ち既存のコードはそのままで再利用できていることになるし、拡張からも独立したままとなってい。
• オブジェクト指向言語における多態性（ポリモーフィズム）と、動的（実行時）バインディング（動的ディスパッチ） という言語の機能を利用して、上記を実現している。
  この機能を持たない言語では、どうしても、コードの改修やリコンパイルが必要になるし、コード間の依存が残ることになる。（その違いを意識して読んでみてください）

ファイルたちはパックしたファイルにて提供するので、ダウンロードして解凍してご利用ください。 なお、JavaですのでJavaの実行環境のインストールは別途必要だ。
ご自身で一部書き替えてお試しするためには実行環境だけではなく開発環境(JDK)のインストールが必要だ。(どちらもOracleのページから無料で入手できる) パックされているファイルの概要は以下となっている。 （ソースコードの解読やクラス構成を決めるときに考察したことなどを含めて、解読の一助に。）

• 設計方針など説明　　・・・どういった方針で作成してきたかを説明している。
• サブフォルダの説明　・・・サブフォルダの目的などを説明している。
• Javaソースコードファイルたち(*.java)
• Javaクラスファイルたち(*.class)
• javadocによる説明ファイル(docサブフォルダに格納)

一部中途半端と思える箇所もあるだろう。それらを含めて、「僕（私）だったら、ここはこう書くかな」「こういった拡張性を考えたら、この部分はこうしておくほうが良いかな」 などを考えるきっかけになり、そして、そう考えることが読者の技術向上になることは間違いないことと思う。