// 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 piR = math.Pi * earthRadius twoPiR = 2 * piR ) // Haversine ... func Haversine(latA, lonA, latB, lonB float64) float64 { φ1 := latA * radians λ1 := lonA * radians φ2 := latB * radians λ2 := lonB * radians Δφ := φ2 - φ1 Δλ := λ2 - λ1 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 } // DistanceFromHaversine... func DistanceFromHaversine(haversine float64) float64 { return earthRadius * 2 * math.Asin(math.Sqrt(haversine)) } // DistanceTo return the distance in meters between two point. func DistanceTo(latA, lonA, latB, lonB float64) (meters float64) { a := Haversine(latA, lonA, latB, lonB) return DistanceFromHaversine(a) } // 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) }