2018-10-11 00:25:40 +03:00
|
|
|
|
// Copyright 2018 Joshua J Baker. All rights reserved.
|
|
|
|
|
// Use of this source code is governed by an MIT-style
|
|
|
|
|
// license that can be found in the LICENSE file.
|
|
|
|
|
|
|
|
|
|
package geo
|
|
|
|
|
|
|
|
|
|
import (
|
|
|
|
|
"math"
|
|
|
|
|
)
|
|
|
|
|
|
|
|
|
|
const (
|
|
|
|
|
earthRadius = 6371e3
|
|
|
|
|
radians = math.Pi / 180
|
|
|
|
|
degrees = 180 / math.Pi
|
2018-10-27 19:23:29 +03:00
|
|
|
|
piR = math.Pi * earthRadius
|
|
|
|
|
twoPiR = 2 * piR
|
2018-10-11 00:25:40 +03:00
|
|
|
|
)
|
|
|
|
|
|
2018-10-27 19:23:29 +03:00
|
|
|
|
// Haversine ...
|
|
|
|
|
func Haversine(latA, lonA, latB, lonB float64) float64 {
|
2018-10-11 00:25:40 +03:00
|
|
|
|
φ1 := latA * radians
|
|
|
|
|
λ1 := lonA * radians
|
|
|
|
|
φ2 := latB * radians
|
|
|
|
|
λ2 := lonB * radians
|
|
|
|
|
Δφ := φ2 - φ1
|
|
|
|
|
Δλ := λ2 - λ1
|
2018-10-27 19:23:29 +03:00
|
|
|
|
sΔφ2 := math.Sin(Δφ / 2)
|
|
|
|
|
sΔλ2 := math.Sin(Δλ / 2)
|
|
|
|
|
return sΔφ2*sΔφ2 + math.Cos(φ1)*math.Cos(φ2)*sΔλ2*sΔλ2
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
// NormalizeDistance ...
|
|
|
|
|
func NormalizeDistance(meters float64) float64 {
|
|
|
|
|
m1 := math.Mod(meters, twoPiR)
|
|
|
|
|
if m1 <= piR {
|
|
|
|
|
return m1
|
|
|
|
|
}
|
|
|
|
|
return twoPiR - m1
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
// DistanceToHaversine ...
|
|
|
|
|
func DistanceToHaversine(meters float64) float64 {
|
|
|
|
|
// convert the given distance to its haversine
|
|
|
|
|
sin := math.Sin(0.5 * meters / earthRadius)
|
|
|
|
|
return sin * sin
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
// DistanceTo return the distance in meteres between two point.
|
|
|
|
|
func DistanceTo(latA, lonA, latB, lonB float64) (meters float64) {
|
|
|
|
|
a := Haversine(latA, lonA, latB, lonB)
|
2018-10-11 00:25:40 +03:00
|
|
|
|
c := 2 * math.Atan2(math.Sqrt(a), math.Sqrt(1-a))
|
|
|
|
|
return earthRadius * c
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
// DestinationPoint return the destination from a point based on a
|
|
|
|
|
// distance and bearing.
|
|
|
|
|
func DestinationPoint(lat, lon, meters, bearingDegrees float64) (
|
|
|
|
|
destLat, destLon float64,
|
|
|
|
|
) {
|
|
|
|
|
// see http://williams.best.vwh.net/avform.htm#LL
|
|
|
|
|
δ := meters / earthRadius // angular distance in radians
|
|
|
|
|
θ := bearingDegrees * radians
|
|
|
|
|
φ1 := lat * radians
|
|
|
|
|
λ1 := lon * radians
|
|
|
|
|
φ2 := math.Asin(math.Sin(φ1)*math.Cos(δ) +
|
|
|
|
|
math.Cos(φ1)*math.Sin(δ)*math.Cos(θ))
|
|
|
|
|
λ2 := λ1 + math.Atan2(math.Sin(θ)*math.Sin(δ)*math.Cos(φ1),
|
|
|
|
|
math.Cos(δ)-math.Sin(φ1)*math.Sin(φ2))
|
|
|
|
|
λ2 = math.Mod(λ2+3*math.Pi, 2*math.Pi) - math.Pi // normalise to -180..+180°
|
|
|
|
|
return φ2 * degrees, λ2 * degrees
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
// BearingTo returns the (initial) bearing from point 'A' to point 'B'.
|
|
|
|
|
func BearingTo(latA, lonA, latB, lonB float64) float64 {
|
|
|
|
|
// tanθ = sinΔλ⋅cosφ2 / cosφ1⋅sinφ2 − sinφ1⋅cosφ2⋅cosΔλ
|
|
|
|
|
// see mathforum.org/library/drmath/view/55417.html for derivation
|
|
|
|
|
|
|
|
|
|
φ1 := latA * radians
|
|
|
|
|
φ2 := latB * radians
|
|
|
|
|
Δλ := (lonB - lonA) * radians
|
|
|
|
|
y := math.Sin(Δλ) * math.Cos(φ2)
|
|
|
|
|
x := math.Cos(φ1)*math.Sin(φ2) - math.Sin(φ1)*math.Cos(φ2)*math.Cos(Δλ)
|
|
|
|
|
θ := math.Atan2(y, x)
|
|
|
|
|
|
|
|
|
|
return math.Mod(θ*degrees+360, 360)
|
|
|
|
|
}
|