はじめに
～あらすじ～
むかしむかしあるところにコードを理屈が分からないままに使うのが嫌なエンジニアがおりました。
ある日彼はとあるサイトで2点間の直線距離を求めるjavascript・コードをみつけました。
コードを読んでも意味が理解できなかった彼はくやしくなって、数学の勉強を（さわりだけ）やってみました。
……ひとに見せるのが目的じゃあない雑記メモなので、とくに解説はしません。
備考
※各引数の内訳
current_lat = 現在地の緯度
current_long = 現在地の経度
target_lat = 目的地の緯度
target_long = 目的地の経度
precision = 精度(くわしくは後述)
※「化成緯度」について
地球楕円体上のある地点を通る赤道面の法線と地球楕円体と
赤道で接する球面との交点を考える。その交点と地球楕円体中心
とを結ぶ直線が赤道面となす角が化成緯度。
http://motochan.sakura.ne.jp/public_html/GISContents/16.htm より
……図を見て「ああこんなもんなんだ」と理解した。筆者は。たぶん。
コード
説明とかは半分くらいコードの中に書いてある。
// 現在地と目的地との距離を計算
function geo_calculate(current_lat, current_long, target_lat, target_long, precision){
// ※公式「測地線航海算法」（Lambert-Andoyerの公式）を用いています。
let distance = 0;
// ２点間の緯度の差と経度の差が両方とも0.00001(＝10万分の一)未満なら、２点間は同じ位置であるとみなし、計算しない
if((Math.abs(current_lat - target_lat) < 0.00001) && (Math.abs(current_long - target_long) < 0.00001)) {
distance = 0;
}else{
// 緯度に円周率を掛けて180で割る
current_lat = current_lat * Math.PI / 180;
current_long = current_long * Math.PI / 180;
target_lat = target_lat * Math.PI / 180;
target_long = target_long * Math.PI / 180;
const QUATOR_RADIUS = 6378140; // <-地球赤道半径(単位はm)
const POSITION_RADIUS = 6356755; // <-地球極半径(単位はm)
const ASPECT_RATIO = (QUATOR_RADIUS - POSITION_RADIUS) / QUATOR_RADIUS; // <-地球の扁平率
// 測地緯度（ここでは、current_latとtarget_latのこと）を化成緯度に変換。
// P1は現在地の化成緯度、P2は目的地の化成緯度。
const P1 = Math.atan((POSITION_RADIUS / QUATOR_RADIUS) * Math.tan(current_lat));
const P2 = Math.atan((POSITION_RADIUS / QUATOR_RADIUS) * Math.tan(target_lat));
// 球面上の距離の計算
const X = Math.acos(Math.sin(P1) * Math.sin(P2) + Math.cos(P1) * Math.cos(P2) * Math.cos(current_long - target_long));
// 距離補正量を算出
const L = (ASPECT_RATIO / 8) * ((Math.sin(X) - X) * Math.pow((Math.sin(P1) + Math.sin(P2)), 2) / Math.pow(Math.cos(X / 2), 2) - (Math.sin(X) - X) * Math.pow(Math.sin(P1) - Math.sin(P2), 2) / Math.pow(Math.sin(X), 2));
// 上で算出した距離補正量をもとに、測地線の長さ（≒地球上の２つの場所の最短距離）を算出する
distance = QUATOR_RADIUS * (X + L);
// 算出された値の精度
// この「精度」は「正確さ」ではなく、「値を表現する細かさ」のことなので注意！
const decimal_no = Math.pow(10, precision);
// ここで出たdistanceの単位はメートル。
distance = Math.round(decimal_no * distance / 1) / decimal_no;
}
return distance;
}
参考文献
http://www2.nc-toyama.ac.jp/~mkawai/lecture/sailing/geodetic/geosail.html
http://motochan.sakura.ne.jp/public_html/GISContents/16.htm
https://qiita.com/damyarou/items/9cb633e844c78307134a
https://ja.wikipedia.org/wiki/%E6%B8%AC%E5%9C%B0%E7%B7%9A