Go言語で緯度経度差から2地点間の距離を計算する

ヒュベニの公式

必要に迫られたため簡便な公式で実装。

機会があればより高精度な計算式でも実装してみたい
http://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/algorithm/bl2st/bl2st.htm

/*
   参考にしたサイト：
   http://www.kashmir3d.com/kash/manual/std_siki.htm
   http://yamadarake.jp/trdi/report000001.html
*/

package main

import (
	"fmt"
	"math"
	"strconv"
)

const (
	EQUATORIAL_RADIUS    = 6378137.0            // 赤道半径 GRS80
	POLAR_RADIUS         = 6356752.314          // 極半径 GRS80
	ECCENTRICITY         = 0.081819191042815790 // 第一離心率 GRS80
)

type Coodinate struct {
	Latitude  float64
	Longitude float64
}

func main() {
	a := Coodinate{35.65500, 139.74472} // 東京
	b := Coodinate{36.10056, 140.09111} // 筑波

	hd := hubenyDistance(a, b)
	distance := strconv.FormatFloat(math.Floor(hd), 'f', 0, 64)

	fmt.Printf("距離: " + distance + " m\n")
}

func hubenyDistance(src Coodinate, dst Coodinate) float64 {
	dx := degree2radian(dst.Longitude - src.Longitude)
	dy := degree2radian(dst.Latitude - src.Latitude)
	my := degree2radian((src.Latitude + dst.Latitude) / 2)

	W := math.Sqrt(1 - (Power2(ECCENTRICITY) * Power2(math.Sin(my)))) // 卯酉線曲率半径の分母
	m_numer := EQUATORIAL_RADIUS * (1 - Power2(ECCENTRICITY))         // 子午線曲率半径の分子

	M := m_numer / math.Pow(W, 3) // 子午線曲率半径
	N := EQUATORIAL_RADIUS / W    // 卯酉線曲率半径

	d := math.Sqrt(Power2(dy*M) + Power2(dx*N*math.Cos(my)))

	return d
}

func degree2radian(x float64) float64 {
	return x * math.Pi / 180
}

func Power2(x float64) float64 {
	return math.Pow(x, 2)
}