概要
以前、GoogleAppsScriptを用いて緯度経度から２点間の距離を求める関数を作成した。
Qiita記事：【GAS】緯度経度から２点間の距離を求める
長距離で誤差が大きくなったり、緯度経度に負の値がある場合に計算できないなどの問題を抱えつつ、個人的用途の都内での使用に限れば使えていたため問題を放置していた。
当時は「ヒュベニの公式」と言う数式で計算していたけど、他にも球面上の線長を求める方法があるとのことで新たに試してみたのが「測地線航海算法」。
そしたら諸々の問題が解決したっぽいので、改めて投稿^^)b
距離計算
ヒュベニ、航海算法ともに、ググると専門的なページがいっぱい出てくるので理屈や公式などは割愛。
GoogleMap測地系
球面上の線長を求めるには、楕円体の長半径と短半径が必要。地球の長半径（赤道半径）と短半径（極半径）には時代や測地系などによって種類があるが、ここではGoogleMapの距離計算に近づけるためにGoogleMapで採用されている定数を使用。
ちなみにGoogleMapでは長半径短半径ともに6371008(m)と、地球を真球として距離を求めている模様。
Qiita記事：【GoogleMap】地球を半径6371008mの球体として距離計算していた
計算結果！
ヒュベニの公式で明らかに結果がおかしかった「国会議事堂～ホワイトハウス」の距離も誤差無し！「東京タワー～スカイツリー」で誤差があるように見えるけど、有効数字的には問題なし。
やったぜ╭( ・ㅂ・)و ̑̑ ｸﾞｯ !
コード
A地点の緯度、A地点の経度、B地点の緯度、B地点の経度を引数として、２点間の距離(km)を返す関数を作成。
function geoSailing(latiA,longA,latiB,longB) {
//定数（GoogleMap系）
const A = 6371008; //長半径
const B = 6371008; //短半径
//定数（ベッセル楕円体 旧日本測地系）
//const A = 6377397.16; //長半径
//const B = 6356079; //短半径
//定数（GRS80 世界測地系）
//const A = 6378137; //長半径
//const B = 6356752.31414; //短半径
//扁平率（F)
var F = (A-B)/A;
//緯度経度をラジアンに変換
var latA = latiA * (Math.PI/180);
var lngA = longA * (Math.PI/180);
var latB = latiB * (Math.PI/180);
var lngB = longB * (Math.PI/180);
//化成緯度（Φ）
var phiA = Math.atan(B/A * Math.tan(latA));
var phiB = Math.atan(B/A * Math.tan(latB));
//球面上の距離（X）
var f1 = Math.sin(phiA) * Math.sin(phiB);
var f2 = Math.cos(phiA) * Math.cos(phiB);
var f3 = Math.cos(lngA - lngB);
var X = Math.acos(f1 + f2 * f3);
//距離補正量（Δρ）
var f4 = Math.sin(X)-X;
var f5 = Math.sin(phiA) + Math.sin(phiB);
var f6 = Math.cos(X/2) * Math.cos(X/2);
var f7 = Math.sin(X)+X;
var f8 = Math.sin(phiA) - Math.sin(phiB);
var f9 = Math.sin(X/2) * Math.sin(X/2);
var drho = F/8 * (f4 * f5*f5 / f6 - f7 * f8*f8 / f9);
//測地線長（ρ）
var rho = A*(X+drho); //距離(m)
var rhokm = rho/1000; //距離(km)
return rhokm;
}
算数苦手なので、お見苦しい途中式だらけに^^;
こちらのコードを１行に圧縮する記事(丿￣ο￣)丿 続編：【ワンライナー】緯度経度から１行で距離計算
検討
コード内の定数を別の測地系に変えれば、その測地系に合わせた計算ができる。
上記の地点間の距離を、国土地理院の距離計算サイトとGASの両方で計算。国土地理院、ヒュベニ、航海算法、すべてベッセルの定数を使用したところ、航海算法で完璧な計算結果！
利用する測地系に合わせた定数で航海算法を使うのが良さそう^^