WIP: Great circle calculations

README.md

@@ -98,6 +98,12 @@ distance(x_lla, y_lla)                   # 401.54 meters
 ```
 (assuming the `wgs84` datum, which can be configured in `distance(x, y, datum)`).
 
+Or the great circle distance can equivalently be calculated:
+```julia
+z_lla = LatLon(51.5077, -0.1277) # Nelson's Column, London, UK
+g = GreatCircle(wgs84)           # Construct a great circle cache
+g(x_lla, z_lla).dist             # 16521.975 km
+```
 
 ## Basic Terminology
 
@@ -444,7 +450,9 @@ counterparts.
 
 ### Distance
 
-Currently, the only defined distance measure is the Cartesian distance,
+#### `distance`
+
+The simplest distance measure is the Cartesian distance,
 `distance(x, y, [datum = wgs84])`, which works for all combinations of types for
 `x` and `y` - except that the UTM zone and hemisphere must also be provided
 for `UTM` types, as in `distance(utm1, utm2, zone, hemisphere, [datum = wgs84])`
@@ -453,5 +461,29 @@ conversion to `ECEF`).
 
 This is the only function currently in
 *Geodesy* which takes a default datum, and *should* be relatively accurate for
-close points where Cartesian distances may be most important. Future work
-may focus on geodesics and related calculations (contributions welcome!).
+close points where Cartesian distances may be most important.
+
+#### `GreatCircle`
+
+Great circle distances, which are calculated along the shortest path along
+the surface of the ellipsoid between points, are calculated by creating a
+`GreatCircle` specfiying the datum to use, like `g = GreatCircle(wgs84)`.
+This creates a cache of values for the ellipsoid, and can be used in two ways:
+
+1. Compute the distance between two points: `g(p1, p2)`.
+2. Compute the end point from moving a set distance from a starting point
+   along a set azimuth: `g(p1, azi=45, dist=100_000)`.
+
+As well as distance, other properties of the great circle arc such as
+azimuth are found.
+
+Case (1) returns a named tuple of azimuth from `p1` to `p2` (`azi` °),
+backazimuth from `p2` to `p1` (`baz`, °), the distance (`dist`, m) and angular
+distance between points on the equivalent sphere (`angle`, °).
+
+Case (2) returns a named tuple of final longitude and latitude
+(`lon`, and `lat`, °), backazimuth from the final point (`baz`, °),
+as well as distance and angular distance (`dist`, m, and `angle`, °).
+
+The great circle calculations in *Geodesy* are ported from *GeographicLib*
+and have the equivalent accuracy as the original library.

src/Geodesy.jl

@@ -48,7 +48,10 @@ export
     ITRF, GDA94,
 
     # UTM helpers
-    utm_zone
+    utm_zone,
+
+    # Great circle constructors
+    GreatCircle
 
 include("ellipsoids.jl")
 include("datums.jl")
@@ -58,6 +61,7 @@ include("polar_stereographic.jl")
 include("transformations.jl")
 include("conversion.jl")
 include("distances.jl")
+include("geodesics.jl")
 include("utm.jl")
 include("datum_transformations.jl")
 

src/GeographicLib/Accumulators.jl

@@ -0,0 +1,87 @@
+module Accumulators
+
+import ..Math
+
+"""Like math.fsum, but allows a running sum"""
+mutable struct Accumulator{T<:AbstractFloat}
+    _s::T
+    _t::T
+    function Accumulator(_s::T1, _t::T2) where {T1,T2}
+    	T = promote_type(float(T1), float(T2))
+    	Accumulator{T}(T(_s), T(_t))
+    end
+    Accumulator{T}(_s, _t) where T = new{T}(T(_s), T(_t))
+end
+
+"""Constructor"""
+Accumulator(y=0.0) = Accumulator(y, 0.0)
+Accumulator(acc::Accumulator) = Accumulator(acc._s, acc._t)
+
+Base.:(==)(a::Accumulator, b::Accumulator) = a._s == b._s && a._t == b._t
+
+"""Set value from argument"""
+Set!(self::Accumulator, y::Accumulator) = ((self._s, self._t) = (y._s, y._t); self)
+Set!(self::Accumulator, y) = ((self._s, self._t) = (y, 0); self)
+
+
+"""Add a value"""
+function Add!(self, y)
+  # Here's Shewchuk's solution...
+  # hold exact sum as [s, t, u]
+  y, u = Math.sum(y, self._t)             # Accumulate starting at
+  self._s, self._t = Math.sum(y, self._s) # least significant end
+  # Start is _s, _t decreasing and non-adjacent.  Sum is now (s + t + u)
+  # exactly with s, t, u non-adjacent and in decreasing order (except
+  # for possible zeros).  The following code tries to normalize the
+  # result.  Ideally, we want _s = round(s+t+u) and _u = round(s+t+u -
+  # _s).  The follow does an approximate job (and maintains the
+  # decreasing non-adjacent property).  Here are two "failures" using
+  # 3-bit floats:
+  #
+  # Case 1: _s is not equal to round(s+t+u) -- off by 1 ulp
+  # [12, -1] - 8 -> [4, 0, -1] -> [4, -1] = 3 should be [3, 0] = 3
+  #
+  # Case 2: _s+_t is not as close to s+t+u as it shold be
+  # [64, 5] + 4 -> [64, 8, 1] -> [64,  8] = 72 (off by 1)
+  #                    should be [80, -7] = 73 (exact)
+  #
+  # "Fixing" these problems is probably not worth the expense.  The
+  # representation inevitably leads to small errors in the accumulated
+  # values.  The additional errors illustrated here amount to 1 ulp of
+  # the less significant word during each addition to the Accumulator
+  # and an additional possible error of 1 ulp in the reported sum.
+  #
+  # Incidentally, the "ideal" representation described above is not
+  # canonical, because _s = round(_s + _t) may not be true.  For
+  # example, with 3-bit floats:
+  #
+  # [128, 16] + 1 -> [160, -16] -- 160 = round(145).
+  # But [160, 0] - 16 -> [128, 16] -- 128 = round(144).
+  #
+  if self._s == 0            # This implies t == 0,
+    self._s = u               # so result is u
+  else
+    self._t += u              # otherwise just accumulate u to t.
+  end
+  self
+end
+
+"""Return sum + y"""
+function Sum(self, y = 0.0)
+  if y == 0.0
+    return self._s
+  else
+    b = Accumulator(self)
+    Add!(b, y)
+    return b._s
+  end
+end
+
+"""Negate sum"""
+function Negate!(self)
+  self._s *= -1
+  self._t *= -1
+  self
+end
+
+end # module

src/GeographicLib/Constants.jl

@@ -0,0 +1,8 @@
+module Constants
+
+"Semimajor radius of WGS84 ellipsoid in m"
+const WGS84_a = 6378137.0
+"Flattening of EGS84 ellipsoid"
+const WGS84_f = 1/298.257223563
+
+end # module

src/GeographicLib/GeodesicCapability.jl

@@ -0,0 +1,27 @@
+module GeodesicCapability
+
+const CAP_NONE = 0
+const CAP_C1   = 1 << 0
+const CAP_C1p  = 1 << 1
+const CAP_C2   = 1 << 2
+const CAP_C3   = 1 << 3
+const CAP_C4   = 1 << 4
+const CAP_ALL  = Int(0x1F)
+const CAP_MASK = CAP_ALL
+const OUT_ALL  = Int(0x7F80)
+const OUT_MASK = Int(0xFF80)             # Includes LONG_UNROLL
+const EMPTY         = 0
+const LATITUDE      = 1 << 7  | CAP_NONE
+const LONGITUDE     = 1 << 8  | CAP_C3
+const AZIMUTH       = 1 << 9  | CAP_NONE
+const DISTANCE      = 1 << 10 | CAP_C1
+const STANDARD      = LATITUDE | LONGITUDE | AZIMUTH | DISTANCE
+const DISTANCE_IN   = 1 << 11 | CAP_C1 | CAP_C1p
+const REDUCEDLENGTH = 1 << 12 | CAP_C1 | CAP_C2
+const GEODESICSCALE = 1 << 13 | CAP_C1 | CAP_C2
+const AREA          = 1 << 14 | CAP_C4
+const LONG_UNROLL   = 1 << 15
+const LONG_NOWRAP   = LONG_UNROLL       # For backwards compatibility only
+const ALL           = OUT_ALL | CAP_ALL
+
+end # module