diff --git a/README.md b/README.md index e53dee6..0476082 100644 --- a/README.md +++ b/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. diff --git a/src/Geodesy.jl b/src/Geodesy.jl index f5a6e38..9121fd9 100644 --- a/src/Geodesy.jl +++ b/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") diff --git a/src/GeographicLib/Accumulators.jl b/src/GeographicLib/Accumulators.jl new file mode 100644 index 0000000..0af442f --- /dev/null +++ b/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 diff --git a/src/GeographicLib/Constants.jl b/src/GeographicLib/Constants.jl new file mode 100644 index 0000000..66f2e8e --- /dev/null +++ b/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 diff --git a/src/GeographicLib/GeodesicCapability.jl b/src/GeographicLib/GeodesicCapability.jl new file mode 100644 index 0000000..13c5eee --- /dev/null +++ b/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 diff --git a/src/GeographicLib/GeodesicLines.jl b/src/GeographicLib/GeodesicLines.jl new file mode 100644 index 0000000..893618c --- /dev/null +++ b/src/GeographicLib/GeodesicLines.jl @@ -0,0 +1,486 @@ +module GeodesicLines + +import ..Math, ..GeodesicCapability, ..Geodesics, ..Result + +export GeodesicLine, SetArc, SetDistance, Position, ArcPosition + +struct GeodesicLine + a::Float64 + f::Float64 + _b::Float64 + _c2::Float64 + _f1::Float64 + """the capabilities (readonly)""" + caps::Int + """the latitude of the first point in degrees (readonly)""" + lat1::Float64 + """the longitude of the first point in degrees (readonly)""" + lon1::Float64 + """the azimuth at the first point in degrees (readonly)""" + azi1::Float64 + """the sine of the azimuth at the first point (readonly)""" + salp1::Float64 + """the cosine of the azimuth at the first point (readonly)""" + calp1::Float64 + _dn1::Float64 + _salp0::Float64 + _calp0::Float64 + _ssig1::Float64 + _csig1::Float64 + _somg1::Float64 + _comg1::Float64 + _k2::Float64 + _A1m1::Float64 + _C1a::Vector{Float64} + _B11::Float64 + _stau1::Float64 + _ctau1::Float64 + _C1pa::Vector{Float64} + _A2m1::Float64 + _C2a::Vector{Float64} + _B21::Float64 + _C3a::Vector{Float64} + _A3c::Float64 + _B31::Float64 + _C4a::Vector{Float64} + _A4::Float64 + _B41::Float64 + """the distance between point 1 and point 3 in meters (readonly)""" + s13::Float64 + """the arc length between point 1 and point 3 in degrees (readonly)""" + a13::Float64 +end + +# Treat the GeodesicLine struct as a scalar +Base.Broadcast.broadcastable(line::GeodesicLine) = Ref(line) + +""" + GeodesicLine(geod::Geodesic, lat1, lon1, azi1; caps, salp1, calp1) -> line + +Create a `GeodesicLine` with starting latitude `lat1`° and longitude `lon1`°, +and azimuth `azi1`°. + +Control the capabilities of the `line` with `caps`. Optionally specify +the sine and cosine of the azimuth at point 1, respectively `salp1` and `calp1`. +""" +function GeodesicLine(geod::Geodesics.Geodesic, lat1, lon1, azi1; + caps = GeodesicCapability.STANDARD | GeodesicCapability.DISTANCE_IN, + salp1 = Math.nan, calp1 = Math.nan) + self_a = geod.a + self_f = geod.f + self__b = geod._b + self__c2 = geod._c2 + self__f1 = geod._f1 + self_caps = (caps | Geodesics.LATITUDE | Geodesics.AZIMUTH | + Geodesics.LONG_UNROLL) + + # Guard against underflow in salp0 + self_lat1 = Math.LatFix(lat1) + self_lon1 = lon1 + if Math.isnan(salp1) || Math.isnan(calp1) + self_azi1 = Math.AngNormalize(azi1) + self_salp1, self_calp1 = Math.sincosd(Math.AngRound(azi1)) + else + self_azi1 = azi1 + self_salp1 = salp1 + self_calp1 = calp1 + end + + # real cbet1, sbet1 + sbet1, cbet1 = Math.sincosd(Math.AngRound(lat1)) + sbet1 *= self__f1 + # Ensure cbet1 = +epsilon at poles + sbet1, cbet1 = Math.norm(sbet1, cbet1) + cbet1 = max(Geodesics.tiny_, cbet1) + self__dn1 = sqrt(1 + geod._ep2 * Math.sq(sbet1)) + + # Evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0), + self__salp0 = self_salp1 * cbet1 # alp0 in [0, pi/2 - |bet1|] + # Alt: calp0 = hypot(sbet1, calp1 * cbet1). The following + # is slightly better (consider the case salp1 = 0). + self__calp0 = hypot(self_calp1, self_salp1 * sbet1) + # Evaluate sig with tan(bet1) = tan(sig1) * cos(alp1). + # sig = 0 is nearest northward crossing of equator. + # With bet1 = 0, alp1 = pi/2, we have sig1 = 0 (equatorial line). + # With bet1 = pi/2, alp1 = -pi, sig1 = pi/2 + # With bet1 = -pi/2, alp1 = 0 , sig1 = -pi/2 + # Evaluate omg1 with tan(omg1) = sin(alp0) * tan(sig1). + # With alp0 in (0, pi/2], quadrants for sig and omg coincide. + # No atan2(0,0) ambiguity at poles since cbet1 = +epsilon. + # With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi. + self__ssig1 = sbet1 + self__somg1 = self__salp0 * sbet1 + self__csig1 = self__comg1 = ((sbet1 != 0 || self_calp1 != 0) ? cbet1 * self_calp1 : 1) + # sig1 in (-pi, pi] + self__ssig1, self__csig1 = Math.norm(self__ssig1, self__csig1) + # No need to normalize + # self__somg1, self__comg1 = Math.norm(self__somg1, self__comg1) + + self__k2 = Math.sq(self__calp0) * geod._ep2 + eps = self__k2 / (2 * (1 + sqrt(1 + self__k2)) + self__k2) + + if (self_caps & Geodesics.CAP_C1) > 0 + self__A1m1 = Geodesics._A1m1f(eps) + self__C1a = Vector{Float64}(undef, Geodesics.nC1_ + 1) #list(range(Geodesics.nC1_ + 1)) + Geodesics._C1f(eps, self__C1a) + self__B11 = Geodesics._SinCosSeries(true, self__ssig1, self__csig1, self__C1a) + s = sin(self__B11) + c = cos(self__B11) + # tau1 = sig1 + B11 + self__stau1 = self__ssig1 * c + self__csig1 * s + self__ctau1 = self__csig1 * c - self__ssig1 * s + # Not necessary because C1pa reverts C1a + # _B11 = -_SinCosSeries(true, _stau1, _ctau1, _C1pa) + else + self__A1m1 = self__B11 = self__stau1 = self__ctau1 = Math.nan + self__C1a = Float64[] + end + + if (self_caps & Geodesics.CAP_C1p) > 0 + self__C1pa = Vector{Float64}(undef, Geodesics.nC1p_ + 1) #list(range(Geodesics.nC1p_ + 1)) + Geodesics._C1pf(eps, self__C1pa) + else + self__C1pa = Float64[] + end + + if (self_caps & Geodesics.CAP_C2) > 0 + self__A2m1 = Geodesics._A2m1f(eps) + self__C2a = Vector{Float64}(undef, Geodesics.nC2_ + 1) #list(range(Geodesics.nC2_ + 1)) + Geodesics._C2f(eps, self__C2a) + self__B21 = Geodesics._SinCosSeries(true, self__ssig1, self__csig1, self__C2a) + else + self__A2m1 = self__B21 = Math.nan + self__C2a = Float64[] + end + + if (self_caps & Geodesics.CAP_C3) > 0 + self__C3a = Vector{Float64}(undef, Geodesics.nC3_) #list(range(Geodesics.nC3_)) + Geodesics._C3f(geod, eps, self__C3a) + self__A3c = -self_f * self__salp0 * Geodesics._A3f(geod, eps) + self__B31 = Geodesics._SinCosSeries(true, self__ssig1, self__csig1, self__C3a) + else + self__C3a = Float64[] + self__A3c = self__B31 = Math.nan + end + + if (self_caps & Geodesics.CAP_C4) > 0 + self__C4a = Vector{Float64}(undef, Geodesics.nC4_) #list(range(Geodesics.nC4_)) + Geodesics._C4f(geod, eps, self__C4a) + # Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0) + self__A4 = Math.sq(self_a) * self__calp0 * self__salp0 * geod._e2 + self__B41 = Geodesics._SinCosSeries(false, self__ssig1, self__csig1, self__C4a) + else + self__C4a = Float64[] + self__A4 = self__B41 = Math.nan + end + + self_s13 = Math.nan + self_a13 = Math.nan + + GeodesicLine(self_a, self_f, self__b, self__c2, self__f1, self_caps, + self_lat1, self_lon1, self_azi1, self_salp1, self_calp1, + self__dn1, self__salp0, self__calp0, self__ssig1, self__csig1, + self__somg1, self__comg1, self__k2, self__A1m1, self__C1a, + self__B11, self__stau1, self__ctau1, self__C1pa, self__A2m1, + self__C2a, self__B21, self__C3a, self__A3c, self__B31, self__C4a, + self__A4, self__B41, self_s13, self_a13) +end + +"""Private: General solution of position along geodesic""" +function _GenPosition(self::GeodesicLine, arcmode, s12_a12, outmask) + s12_a12 = Float64(s12_a12) + + a12 = lat2 = lon2 = azi2 = s12 = m12 = M12 = M21 = S12 = Math.nan + outmask &= self.caps & Geodesics.OUT_MASK + if ! (arcmode || + (self.caps & (Geodesics.OUT_MASK & Geodesics.DISTANCE_IN)) > 0) + # Uninitialized or impossible distance calculation requested + return a12, lat2, lon2, azi2, s12, m12, M12, M21, S12 + end + + # Avoid warning about uninitialized B12. + B12 = 0.0 + AB1 = 0.0 + if arcmode + # Interpret s12_a12 as spherical arc length + sig12 = deg2rad(s12_a12) + ssig12, csig12 = Math.sincosd(s12_a12) + else + # Interpret s12_a12 as distance + tau12 = s12_a12 / (self._b * (1 + self._A1m1)) + s = sin(tau12) + c = cos(tau12) + # tau2 = tau1 + tau12 + B12 = - Geodesics._SinCosSeries(true, + self._stau1 * c + self._ctau1 * s, + self._ctau1 * c - self._stau1 * s, + self._C1pa) + sig12 = tau12 - (B12 - self._B11) + ssig12 = sin(sig12) + csig12 = cos(sig12) + if abs(self.f) > 0.01 + # Reverted distance series is inaccurate for |f| > 1/100, so correct + # sig12 with 1 Newton iteration. The following table shows the + # approximate maximum error for a = WGS_a() and various f relative to + # GeodesicExact. + # erri = the error in the inverse solution (nm) + # errd = the error in the direct solution (series only) (nm) + # errda = the error in the direct solution (series + 1 Newton) (nm) + # + # f erri errd errda + # -1/5 12e6 1.2e9 69e6 + # -1/10 123e3 12e6 765e3 + # -1/20 1110 108e3 7155 + # -1/50 18.63 200.9 27.12 + # -1/100 18.63 23.78 23.37 + # -1/150 18.63 21.05 20.26 + # 1/150 22.35 24.73 25.83 + # 1/100 22.35 25.03 25.31 + # 1/50 29.80 231.9 30.44 + # 1/20 5376 146e3 10e3 + # 1/10 829e3 22e6 1.5e6 + # 1/5 157e6 3.8e9 280e6 + ssig2 = self._ssig1 * csig12 + self._csig1 * ssig12 + csig2 = self._csig1 * csig12 - self._ssig1 * ssig12 + B12 = Geodesics._SinCosSeries(true, ssig2, csig2, self._C1a) + serr = ((1 + self._A1m1) * (sig12 + (B12 - self._B11)) - s12_a12 / self._b) + sig12 = sig12 - serr / sqrt(1 + self._k2 * Math.sq(ssig2)) + ssig12 = sin(sig12) + csig12 = cos(sig12) + # Update B12 below + end + end + + # real omg12, lam12, lon12 + # real ssig2, csig2, sbet2, cbet2, somg2, comg2, salp2, calp2 + # sig2 = sig1 + sig12 + ssig2 = self._ssig1 * csig12 + self._csig1 * ssig12 + csig2 = self._csig1 * csig12 - self._ssig1 * ssig12 + dn2 = sqrt(1 + self._k2 * Math.sq(ssig2)) + if outmask & (Geodesics.DISTANCE | Geodesics.REDUCEDLENGTH | Geodesics.GEODESICSCALE) > 0 + if arcmode || abs(self.f) > 0.01 + B12 = Geodesics._SinCosSeries(true, ssig2, csig2, self._C1a) + end + AB1 = (1 + self._A1m1) * (B12 - self._B11) + end + # sin(bet2) = cos(alp0) * sin(sig2) + sbet2 = self._calp0 * ssig2 + # Alt: cbet2 = hypot(csig2, salp0 * ssig2) + cbet2 = hypot(self._salp0, self._calp0 * csig2) + if cbet2 == 0 + # I.e., salp0 = 0, csig2 = 0. Break the degeneracy in this case + cbet2 = csig2 = Geodesics.tiny_ + end + # tan(alp0) = cos(sig2)*tan(alp2) + salp2 = self._salp0 + calp2 = self._calp0 * csig2 # No need to normalize + + if (outmask & Geodesics.DISTANCE) > 0 + s12 = arcmode ? self._b * ((1 + self._A1m1) * sig12 + AB1) : s12_a12 + end + + if (outmask & Geodesics.LONGITUDE) > 0 + # tan(omg2) = sin(alp0) * tan(sig2) + somg2 = self._salp0 * ssig2 + comg2 = csig2 # No need to normalize + E = Math.copysign(1, self._salp0) # East or west going? + # omg12 = omg2 - omg1 + omg12 = ((outmask & Geodesics.LONG_UNROLL) > 0 ? + E * (sig12 + - (atan( ssig2, csig2) - + atan( self._ssig1, self._csig1)) + + (atan(E * somg2, comg2) - + atan(E * self._somg1, self._comg1))) : + atan(somg2 * self._comg1 - comg2 * self._somg1, + comg2 * self._comg1 + somg2 * self._somg1)) + lam12 = omg12 + self._A3c * ( + sig12 + (Geodesics._SinCosSeries(true, ssig2, csig2, self._C3a) - self._B31)) + lon12 = rad2deg(lam12) + lon2 = ((outmask & Geodesics.LONG_UNROLL) > 0 ? + self.lon1 + lon12 : + Math.AngNormalize(Math.AngNormalize(self.lon1) + + Math.AngNormalize(lon12))) + end + + if (outmask & Geodesics.LATITUDE) > 0 + lat2 = Math.atan2d(sbet2, self._f1 * cbet2) + end + + if (outmask & Geodesics.AZIMUTH) > 0 + azi2 = Math.atan2d(salp2, calp2) + end + + if (outmask & (Geodesics.REDUCEDLENGTH | Geodesics.GEODESICSCALE)) > 0 + B22 = Geodesics._SinCosSeries(true, ssig2, csig2, self._C2a) + AB2 = (1 + self._A2m1) * (B22 - self._B21) + J12 = (self._A1m1 - self._A2m1) * sig12 + (AB1 - AB2) + if (outmask & Geodesics.REDUCEDLENGTH) > 0 + # Add parens around (_csig1 * ssig2) and (_ssig1 * csig2) to ensure + # accurate cancellation in the case of coincident points. + m12 = self._b * (( dn2 * (self._csig1 * ssig2) - + self._dn1 * (self._ssig1 * csig2)) + - self._csig1 * csig2 * J12) + end + if (outmask & Geodesics.GEODESICSCALE) > 0 + t = (self._k2 * (ssig2 - self._ssig1) * + (ssig2 + self._ssig1) / (self._dn1 + dn2)) + M12 = csig12 + (t * ssig2 - csig2 * J12) * self._ssig1 / self._dn1 + M21 = csig12 - (t * self._ssig1 - self._csig1 * J12) * ssig2 / dn2 + end + end + + if (outmask & Geodesics.AREA) > 0 + B42 = Geodesics._SinCosSeries(false, ssig2, csig2, self._C4a) + # real salp12, calp12 + if self._calp0 == 0 || self._salp0 == 0 + # alp12 = alp2 - alp1, used in atan2 so no need to normalize + salp12 = salp2 * self.calp1 - calp2 * self.salp1 + calp12 = calp2 * self.calp1 + salp2 * self.salp1 + else + # tan(alp) = tan(alp0) * sec(sig) + # tan(alp2-alp1) = (tan(alp2) -tan(alp1)) / (tan(alp2)*tan(alp1)+1) + # = calp0 * salp0 * (csig1-csig2) / (salp0^2 + calp0^2 * csig1*csig2) + # If csig12 > 0, write + # csig1 - csig2 = ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1) + # else + # csig1 - csig2 = csig1 * (1 - csig12) + ssig12 * ssig1 + # No need to normalize + salp12 = self._calp0 * self._salp0 * ( + csig12 <= 0 ? + self._csig1 * (1 - csig12) + ssig12 * self._ssig1 : + ssig12 * (self._csig1 * ssig12 / (1 + csig12) + self._ssig1)) + calp12 = (Math.sq(self._salp0) + + Math.sq(self._calp0) * self._csig1 * csig2) + end + S12 = (self._c2 * atan(salp12, calp12) + + self._A4 * (B42 - self._B41)) + end + + a12 = arcmode ? s12_a12 : rad2deg(sig12) + a12, lat2, lon2, azi2, s12, m12, M12, M21, S12 +end + +""" + Position(line::GeodesicLine, s12, outmask=STANDARD) -> result::Result + +Find the position on the line given `s12`, the distance from the +first point to the second in metres. + +The default value of `outmask` is `STANDARD`, i.e., the `lat1`, +`lon1`, `azi1`, `lat2`, `lon2`, `azi2`, `s12`, `a12` entries are +returned. The `GeodesicLine` object must have been constructed with +the `DISTANCE_IN` capability. +""" +function Position(self::GeodesicLine, s12; outmask = GeodesicCapability.STANDARD) + result = Result() + result.lat1 = self.lat1 + result.lon1 = (outmask & Geodesics.LONG_UNROLL) > 0 ? + self.lon1 : + Math.AngNormalize(self.lon1) + result.azi1 = self.azi1 + result.s12 = s12 + a12, lat2, lon2, azi2, s12, m12, M12, M21, S12 = _GenPosition(self, false, s12, outmask) + outmask &= Geodesics.OUT_MASK + result.a12 = a12 + (outmask & Geodesics.LATITUDE) > 0 && (result.lat2 = lat2) + (outmask & Geodesics.LONGITUDE) > 0 && (result.lon2 = lon2) + (outmask & Geodesics.AZIMUTH) > 0 && (result.azi2 = azi2) + (outmask & Geodesics.REDUCEDLENGTH) > 0 && (result.m12 = m12) + if (outmask & Geodesics.GEODESICSCALE) > 0 + result.M12 = M12 + result.M21 = M21 + end + (outmask & Geodesics.AREA) > 0 && (result.S12 = S12) + result +end + +""" + ArcPosition(line::GeodesicLine, a12, outmask=STANDARD) -> result::Result + +Find the position on the line given `a12`, the spherical arc length from +the first point to the second in degrees. + +The default value of `outmask` is `STANDARD`, i.e., the `lat1`, +`lon1`, `azi1`, `lat2`, `lon2`, `azi2`, `s12`, `a12` entries are +returned. +""" +function ArcPosition(self::GeodesicLine, a12; outmask = GeodesicCapability.STANDARD) + result = Result() + result.lat1 = self.lat1 + result.lon1 = (outmask & Geodesics.LONG_UNROLL) > 0 ? + self.lon1 : + Math.AngNormalize(self.lon1) + result.azi1 = self.azi1 + result.a12 = a12 + a12, lat2, lon2, azi2, s12, m12, M12, M21, S12 = _GenPosition(self, true, a12, outmask) + outmask &= Geodesics.OUT_MASK + (outmask & Geodesics.DISTANCE) > 0 && (result.s12 = s12) + (outmask & Geodesics.LATITUDE) > 0 && (result.lat2 = lat2) + (outmask & Geodesics.LONGITUDE) > 0 && (result.lon2 = lon2) + (outmask & Geodesics.AZIMUTH) > 0 && (result.azi2 = azi2) + (outmask & Geodesics.REDUCEDLENGTH) > 0 && (result.m12 = m12) + if (outmask & Geodesics.GEODESICSCALE) > 0 + result.M12 = M12 + result.M21 = M21 + end + (outmask & Geodesics.AREA) > 0 && (result.S12 = S12) + result +end + +""" + SetDistance(line::GeodesicLine, s13) -> line′::GeodesicLine + +Specify the position of point 3 in terms of distance +from point 1 to point 3 in meters + +Return a new `GeodesicLine` with `s13` and `a13` set. +""" +function SetDistance(self::GeodesicLine, s13) + self = GeodesicLine(self.a, self.f, self._b, self._c2, self._f1, self.caps, + self.lat1, self.lon1, self.azi1, self.salp1, self.calp1, + self._dn1, self._salp0, self._calp0, self._ssig1, self._csig1, + self._somg1, self._comg1, self._k2, self._A1m1, self._C1a, + self._B11, self._stau1, self._ctau1, self._C1pa, self._A2m1, + self._C2a, self._B21, self._C3a, self._A3c, self._B31, self._C4a, + self._A4, self._B41, + s13, self.a13) + a13, _, _, _, _, _, _, _, _ = _GenPosition(self, false, s13, 0) + GeodesicLine(self.a, self.f, self._b, self._c2, self._f1, self.caps, + self.lat1, self.lon1, self.azi1, self.salp1, self.calp1, + self._dn1, self._salp0, self._calp0, self._ssig1, self._csig1, + self._somg1, self._comg1, self._k2, self._A1m1, self._C1a, + self._B11, self._stau1, self._ctau1, self._C1pa, self._A2m1, + self._C2a, self._B21, self._C3a, self._A3c, self._B31, self._C4a, + self._A4, self._B41, + s13, a13) +end + +""" + SetArc(line::GeodesicLine, a13) -> line′::GeodesicLine + +Specify the position of point 3 in terms of spherical arc length +from point 1 to point 3 in degrees. + +Return a new `GeodesicLine` with `a13` and `s13` set. +""" +function SetArc(self::GeodesicLine, a13) + self = GeodesicLine(self.a, self.f, self._b, self._c2, self._f1, self.caps, + self.lat1, self.lon1, self.azi1, self.salp1, self.calp1, + self._dn1, self._salp0, self._calp0, self._ssig1, self._csig1, + self._somg1, self._comg1, self._k2, self._A1m1, self._C1a, + self._B11, self._stau1, self._ctau1, self._C1pa, self._A2m1, + self._C2a, self._B21, self._C3a, self._A3c, self._B31, self._C4a, + self._A4, self._B41, + self.s13, a13) + _, _, _, _, s13, _, _, _, _ = _GenPosition(self, true, a13, Geodesics.DISTANCE) + GeodesicLine(self.a, self.f, self._b, self._c2, self._f1, self.caps, + self.lat1, self.lon1, self.azi1, self.salp1, self.calp1, + self._dn1, self._salp0, self._calp0, self._ssig1, self._csig1, + self._somg1, self._comg1, self._k2, self._A1m1, self._C1a, + self._B11, self._stau1, self._ctau1, self._C1pa, self._A2m1, + self._C2a, self._B21, self._C3a, self._A3c, self._B31, self._C4a, + self._A4, self._B41, + s13, a13) +end + +end # module diff --git a/src/GeographicLib/Geodesics.jl b/src/GeographicLib/Geodesics.jl new file mode 100644 index 0000000..36ee182 --- /dev/null +++ b/src/GeographicLib/Geodesics.jl @@ -0,0 +1,1094 @@ +module Geodesics + +using StaticArrays: @SVector + +import ..Math, ..Constants, ..GeodesicCapability, ..Result + +export Geodesic, ArcDirect, Direct, Inverse + +const GEOGRAPHICLIB_GEODESIC_ORDER = 6 +const nA1_ = GEOGRAPHICLIB_GEODESIC_ORDER +const nC1_ = GEOGRAPHICLIB_GEODESIC_ORDER +const nC1p_ = GEOGRAPHICLIB_GEODESIC_ORDER +const nA2_ = GEOGRAPHICLIB_GEODESIC_ORDER +const nC2_ = GEOGRAPHICLIB_GEODESIC_ORDER +const nA3_ = GEOGRAPHICLIB_GEODESIC_ORDER +const nA3x_ = nA3_ +const nC3_ = GEOGRAPHICLIB_GEODESIC_ORDER +const nC3x_ = (nC3_ * (nC3_ - 1)) ÷ 2 +const nC4_ = GEOGRAPHICLIB_GEODESIC_ORDER +const nC4x_ = (nC4_ * (nC4_ + 1)) ÷ 2 +const maxit1_ = 20 +const maxit2_ = maxit1_ + Math.digits + 10 + +const tiny_ = sqrt(Math.minval) +const tol0_ = Math.epsilon +const tol1_ = 200 * tol0_ +const tol2_ = sqrt(tol0_) +const tolb_ = tol0_ * tol2_ +const xthresh_ = 1000 * tol2_ + +const CAP_NONE = GeodesicCapability.CAP_NONE +const CAP_C1 = GeodesicCapability.CAP_C1 +const CAP_C1p = GeodesicCapability.CAP_C1p +const CAP_C2 = GeodesicCapability.CAP_C2 +const CAP_C3 = GeodesicCapability.CAP_C3 +const CAP_C4 = GeodesicCapability.CAP_C4 +const CAP_ALL = GeodesicCapability.CAP_ALL +const CAP_MASK = GeodesicCapability.CAP_MASK +const OUT_ALL = GeodesicCapability.OUT_ALL +const OUT_MASK = GeodesicCapability.OUT_MASK + +"""No capabilities, no output.""" +const EMPTY = GeodesicCapability.EMPTY +"""Calculate latitude `lat2`.""" +const LATITUDE = GeodesicCapability.LATITUDE +"""Calculate longitude `lon2`.""" +const LONGITUDE = GeodesicCapability.LONGITUDE +"""Calculate azimuths `azi1` and `azi2`.""" +const AZIMUTH = GeodesicCapability.AZIMUTH +"""Calculate distance `s12`.""" +const DISTANCE = GeodesicCapability.DISTANCE +"""All of the above.""" +const STANDARD = GeodesicCapability.STANDARD +"""Allow distance `s12` to be used as input in the direct geodesic problem.""" +const DISTANCE_IN = GeodesicCapability.DISTANCE_IN +"""Calculate reduced length `m12`.""" +const REDUCEDLENGTH = GeodesicCapability.REDUCEDLENGTH +"""Calculate geodesic scales `M12` and `M21`.""" +const GEODESICSCALE = GeodesicCapability.GEODESICSCALE +"""Calculate area `S12`.""" +const AREA = GeodesicCapability.AREA +"""All of the above.""" +const ALL = GeodesicCapability.ALL +"""Unroll longitudes, rather than reducing them to the range [-180°,180°].""" +const LONG_UNROLL = GeodesicCapability.LONG_UNROLL + +struct Geodesic + a::Float64 + f::Float64 + _f1::Float64 + _e2::Float64 + _ep2::Float64 + _n::Float64 + _b::Float64 + _c2::Float64 + _etol2::Float64 + _A3x::Vector{Float64} + _C3x::Vector{Float64} + _C4x::Vector{Float64} +end + +# Treat the Geodesic struct as a scalar +Base.Broadcast.broadcastable(geod::Geodesic) = Ref(geod) + +Base.:(==)(g1::Geodesic, g2::Geodesic) = + all(getfield(g1, f) == getfield(g2, f) for f in fieldnames(Geodesic)) + +""" + Geodesic(a, f) -> geodesic + +Set up an ellipsoid for geodesic calculations. `a` is the semimajor radius of the +ellipsoid, whilst flattening is given by `f`. +""" +function Geodesic(a, f) + f1 = 1 - f + e2 = f*(2 - f) + ep2 = e2/Math.sq(f1) # e2 / (1 - e2) + n = f/(2 - f) + b = a*f1 + # authalic radius squared + term = if e2 == 0 + 1.0 + elseif e2 > 0 + Math.atanh(sqrt(e2))/sqrt(abs(e2)) + else + atan(sqrt(-e2))/sqrt(abs(e2)) + end + c2 = (Math.sq(a) + Math.sq(b) * term) / 2 + # The sig12 threshold for "really short". Using the auxiliary sphere + # solution with dnm computed at (bet1 + bet2) / 2, the relative error in + # the azimuth consistency check is sig12^2 * abs(f) * min(1, 1-f/2) / 2. + # (Error measured for 1/100 < b/a < 100 and abs(f) >= 1/1000. For a given + # f and sig12, the max error occurs for lines near the pole. If the old + # rule for computing dnm = (dn1 + dn2)/2 is used, then the error increases + # by a factor of 2.) Setting this equal to epsilon gives sig12 = etol2. + # Here 0.1 is a safety factor (error decreased by 100) and max(0.001, + # abs(f)) stops etol2 getting too large in the nearly spherical case. + etol2 = 0.1 * tol2_ / sqrt(max(0.001, abs(f)) * min(1.0, 1-f/2) / 2) + !(Math.isfinite(a) && a > 0) && + throw(ArgumentError("Major radius is not positive")) + !(Math.isfinite(b) && b > 0) && + throw(ArgumentError("Minor radius is not positive")) + A3x = Vector{Float64}(undef, nA3x_) + C3x = Vector{Float64}(undef, nC3x_) + C4x = Vector{Float64}(undef, nC4x_) + self = Geodesic(a, f, f1, e2, ep2, n, b, c2, etol2, A3x, C3x, C4x) + _A3coeff(self) + _C3coeff(self) + _C4coeff(self) + self +end + +""" + Inverse(geodesic, lat1, lon1, lat2, lon2, outmask=STANDARD) -> result::Result + +Solve the inverse geodesic problem and return a `Result` containing the +parameters of interest. + +Input arguments: +- `lat1`: latitude of the first point in degrees +- `lon1`: longitude of the first point in degrees +- `lat2`: latitude of the second point in degrees +- `lon2`: longitude of the second point in degrees +- `outmask`: a mask setting which output values are computed (see note below) + +Compute geodesic between (`lat1`, `lon1`) and (`lat2`, `lon2`). +The default value of `outmask` is `Geodesics.STANDARD`, i.e., the `lat1`, +`lon1`, `azi1`, `lat2`, `lon2`, `azi2`, `s12`, `a12` entries are returned. + +### Output mask + +May be any combination of: +`Geodesics.EMPTY`, `Geodesics.LATITUDE`, `Geodesics.LONGITUDE`, +`Geodesics.AZIMUTH`, `Geodesics.DISTANCE`, `Geodesics.STANDARD`, +`Geodesics.DISTANCE_IN`, `Geodesics.REDUCEDLENGTH`, `Geodesics.GEODESICSCALE`, +`Geodesics.AREA`, `Geodesics.ALL` or `Geodesics.LONG_UNROLL`. +See the docstring for each for more information. + +Flags are combined by bitwise or-ing values together, e.g. +`Geodesics.AZIMUTH | Geodesics.DISTANCE`. +""" +function Inverse(self::Geodesic, lat1::T1, lon1::T2, lat2::T3, lon2::T4, + outmask = GeodesicCapability.STANDARD) where {T1,T2,T3,T4} + + T = float(promote_type(Float64, T1, T2, T3, T4)) + lat1, lon1, lat2, lon2 = T.((lat1, lon1, lat2, lon2)) + + a12, s12, salp1, calp1, salp2, calp2, m12, M12, M21, S12 = _GenInverse(self, + lat1, lon1, lat2, lon2, outmask) + outmask &= OUT_MASK + if (outmask & LONG_UNROLL) > 0 + lon12, e = Math.AngDiff(lon1, lon2) + lon2 = (lon1 + lon12) + e + else + lon2 = Math.AngNormalize(lon2) + end + result = Result() + result.lat1 = Math.LatFix(lat1) + result.lon1 = (outmask & LONG_UNROLL) > 0 ? lon1 : Math.AngNormalize(lon1) + result.lat2 = Math.LatFix(lat2) + result.lon2 = lon2 + result.a12 = a12 + if (outmask & DISTANCE) > 0 + result.s12 = s12 + end + if (outmask & AZIMUTH) > 0 + result.azi1 = atand(salp1, calp1) + result.azi2 = atand(salp2, calp2) + end + if (outmask & REDUCEDLENGTH) > 0 + result.m12 = m12 + end + if (outmask & GEODESICSCALE) > 0 + result.M12 = M12 + result.M21 = M21 + end + if (outmask & AREA) > 0 + result.S12 = S12 + end + result +end + +"""Private: Evaluate a trig series using Clenshaw summation.""" +function _SinCosSeries(sinp, sinx::T1, cosx::T2, c) where {T1,T2} + # Evaluate + # y = sinp ? sum(c[i] * sin( 2*i * x), i, 1, n) : + # sum(c[i] * cos((2*i+1) * x), i, 0, n-1) + # using Clenshaw summation. N.B. c[0] is unused for sin series + # Approx operation count = (n + 5) mult and (2 * n + 2) add + T = float(promote_type(T1, T2, eltype(c))) + # T = identity + k = length(c) # Point to one beyond last element + n = k - Int(sinp) + ar = 2 * (cosx - sinx) * (cosx + sinx) # 2 * cos(2 * x) + y1 = T(0) # accumulators for sum + if n & 1 == 1 + k -= 1 + y0 = T(c[k+1]) + else + y0 = T(0) + end + # Now n is even + n = n ÷ 2 + while n != 0 # while n--: + n -= 1 + # Unroll loop x 2, so accumulators return to their original role + k -= 1 + y1 = ar * y0 - y1 + c[k+1] + k -= 1 + y0 = ar * y1 - y0 + c[k+1] + end + sinp ? 2 * sinx * cosx * y0 : # sin(2 * x) * y0 + cosx * (y0 - y1) # cos(x) * (y0 - y1) +end + +"""Private: solve astroid equation.""" +function _Astroid(x, y) + # Solve k^4+2*k^3-(x^2+y^2-1)*k^2-2*y^2*k-y^2 = 0 for positive root k. + # This solution is adapted from Geocentric::Reverse. + p = Math.sq(x) + q = Math.sq(y) + r = (p + q - 1) / 6 + if !(q == 0 && r <= 0) + # Avoid possible division by zero when r = 0 by multiplying equations + # for s and t by r^3 and r, resp. + S = p * q / 4 # S = r^3 * s + r2 = Math.sq(r) + r3 = r * r2 + # The discrimant of the quadratic equation for T3. This is zero on + # the evolute curve p^(1/3)+q^(1/3) = 1 + disc = S * (S + 2 * r3) + u = r + if disc >= 0 + T3 = S + r3 + # Pick the sign on the sqrt to maximize abs(T3). This minimizes loss + # of precision due to cancellation. The result is unchanged because + # of the way the T is used in definition of u. + T3 += T3 < 0 ? -sqrt(disc) : sqrt(disc) # T3 = (r * t)^3 + # N.B. cbrt always returns the real root. cbrt(-8) = -2. + T = Math.cbrt(T3) # T = r * t + # T can be zero; but then r2 / T -> 0. + u += T + (T != 0 ? (r2 / T) : 0.0) + else + # T is complex, but the way u is defined the result is real. + ang = atan(sqrt(-disc), -(S + r3)) + # There are three possible cube roots. We choose the root which + # avoids cancellation. Note that disc < 0 implies that r < 0. + u += 2 * r * cos(ang / 3) + end + v = sqrt(Math.sq(u) + q) # guaranteed positive + # Avoid loss of accuracy when u < 0. + uv = u < 0 ? q/(v - u) : u + v # u+v, guaranteed positive + w = (uv - q) / (2 * v) # positive? + # Rearrange expression for k to avoid loss of accuracy due to + # subtraction. Division by 0 not possible because uv > 0, w >= 0. + k = uv / (sqrt(uv + Math.sq(w)) + w) # guaranteed positive + else # q == 0 && r <= 0 + # y = 0 with |x| <= 1. Handle this case directly. + # for y small, positive root is k = abs(y)/sqrt(1-x^2) + k = 0.0 + end + k +end + +"""Private: return A1-1.""" +function _A1m1f(eps) + coeff = @SVector[1, 4, 64, 0, 256] + m = nA1_ ÷ 2 + t = Math.polyval(m, coeff, 0, Math.sq(eps)) / coeff[m + 2] + (t + eps)/(1 - eps) +end + +"""Private: return C1.""" +function _C1f(eps, c) + coeff = @SVector [-1, 6, -16, 32, -9, 64, -128, 2048, 9, -16, 768, 3, -5, 512, + -7, 1280, -7, 2048] + eps2 = Math.sq(eps) + d = eps + o = 0 + for l in 1:nC1_ # l is index of C1p[l] + m = (nC1_ - l) ÷ 2 # order of polynomial in eps^2 + c[l+1] = d * Math.polyval(m, coeff, o, eps2) / coeff[o + m + 2] + o += m + 2 + d *= eps + end + c +end + +"""Private: return C1'""" +function _C1pf(eps, c) + coeff = @SVector [205, -432, 768, 1536, 4005, -4736, 3840, 12288, -225, 116, 384, + -7173, 2695, 7680, 3467, 7680, 38081, 61440] + eps2 = Math.sq(eps) + d = eps + o = 0 + for l in 1:nC1p_ # l is index of C1p[l] + m = (nC1p_ - l) ÷ 2 # order of polynomial in eps^2 + c[l+1] = d * Math.polyval(m, coeff, o, eps2) / coeff[o + m + 2] + o += m + 2 + d *= eps + end + c +end + +# TODO: Double check that this is correct. Differs from Python library by +# 1.0e-17 for eps = 0.01 +# 8.8e-10 for eps = 0.1 +"""Private: return A2-1""" +function _A2m1f(eps) + coeff = @SVector [-11, -28, -192, 0, 256] + m = nA2_ ÷ 2 + t = Math.polyval(m, coeff, 0, Math.sq(eps)) / coeff[m + 2] + (t - eps) / (1 + eps) +end + +"""Private: return C2""" +function _C2f(eps, c) + coeff = @SVector [1, 2, 16, 32, 35, 64, 384, 2048, 15, 80, 768, 7, 35, 512, + 63, 1280, 77, 2048] + eps2 = Math.sq(eps) + d = eps + o = 0 + for l in 1:nC2_ # l is index of C2[l] + m = (nC2_ - l) ÷ 2 # order of polynomial in eps^2 + c[l+1] = d * Math.polyval(m, coeff, o, eps2) / coeff[o + m + 2] + o += m + 2 + d *= eps + end + c +end + +"""Private: return coefficients for A3""" +function _A3coeff(self::Geodesic) + coeff = (-3, 128, -2, -3, 64, -1, -3, -1, 16, 3, -1, -2, 8, 1, -1, 2, 1, 1) + o = k = 0 + for j in (nA3_-1):-1:0 # coeff of eps^j + m = min(nA3_ - j - 1, j) # order of polynomial in n + self._A3x[k+1] = Math.polyval(m, coeff, o, self._n) / coeff[o + m + 2] + k += 1 + o += m + 2 + end + self._A3x +end + +"""Private: return coefficients for C3""" +function _C3coeff(self::Geodesic) + coeff = [3, 128, 2, 5, 128, -1, 3, 3, 64, -1, 0, 1, 8, -1, 1, 4, 5, 256, + 1, 3, 128, -3, -2, 3, 64, 1, -3, 2, 32, 7, 512, -10, 9, 384, 5, + -9, 5, 192, 7, 512, -14, 7, 512, 21, 2560] + o = k = 0 + for l in 1:(nC3_ - 1) # l is index of C3[l] + for j in (nC3_ - 1):-1:l # coeff of eps^j + m = min(nC3_ - j - 1, j) # order of polynomial in n + self._C3x[k+1] = Math.polyval(m, coeff, o, self._n) / coeff[o + m + 2] + k += 1 + o += m + 2 + end + end + self._C3x +end + +"""Private: return coefficients for C4""" +function _C4coeff(self::Geodesic) + coeff = [97, 15015, 1088, 156, 45045, -224, -4784, 1573, 45045, -10656, 14144, + -4576, -858, 45045, 64, 624, -4576, 6864, -3003, 15015, 100, 208, 572, + 3432, -12012, 30030, 45045, 1, 9009, -2944, 468, 135135, 5792, 1040, + -1287, 135135, 5952, -11648, 9152, -2574, 135135, -64, -624, 4576, + -6864, 3003, 135135, 8, 10725, 1856, -936, 225225, -8448, 4992, -1144, + 225225, -1440, 4160, -4576, 1716, 225225, -136, 63063, 1024, -208, + 105105, 3584, -3328, 1144, 315315, -128, 135135, -2560, 832, 405405, + 128, 99099] + o = k = 0 + for l in 0:(nC4_ - 1) # l is index of C4[l] + for j in (nC4_ - 1):-1:l # coeff of eps^j + m = nC4_ - j - 1 # order of polynomial in n + self._C4x[k+1] = Math.polyval(m, coeff, o, self._n) / coeff[o + m + 2] + k += 1 + o += m + 2 + end + end + self._C4x +end + +"""Private: return A3""" +_A3f(self::Geodesic, eps) = Math.polyval(nA3_ - 1, self._A3x, 0, eps) + +"""Private: return C3""" +function _C3f(self::Geodesic, eps, c) + # Evaluate C3 + # Elements c[1] thru c[nC3_ - 1] are set + mult = oneunit(eps) + o = 0 + for l in 1:(nC3_ - 1) # l is index of C3[l] + m = nC3_ - l - 1 # order of polynomial in eps + mult *= eps + c[l+1] = mult * Math.polyval(m, self._C3x, o, eps) + o += m + 1 + end + c +end + +"""Private: return C4""" +function _C4f(self::Geodesic, eps, c) + # Evaluate C4 coeffs by Horner's method + # Elements c[0] thru c[nC4_ - 1] are set + mult = oneunit(eps) + o = 0 + for l in 0:(nC4_ - 1) # l is index of C4[l] + m = nC4_ - l - 1 # order of polynomial in eps + c[l+1] = mult * Math.polyval(m, self._C4x, o, eps) + o += m + 1 + mult *= eps + end + c +end + +"""Private: return a bunch of lengths""" +function _Lengths(self::Geodesic, eps, sig12, + ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2, outmask, + # Scratch areas of the right size + C1a, C2a) + # Return s12b, m12b, m0, M12, M21, where + # m12b = (reduced length)/_b; s12b = distance/_b, + # m0 = coefficient of secular term in expression for reduced length. + outmask &= OUT_MASK + # outmask & DISTANCE: set s12b + # outmask & REDUCEDLENGTH: set m12b & m0 + # outmask & GEODESICSCALE: set M12 & M21 + + s12b = m12b = m0 = M12 = M21 = Math.nan + if (outmask & (DISTANCE | REDUCEDLENGTH | GEODESICSCALE)) > 0 + A1 = _A1m1f(eps) + _C1f(eps, C1a) + if (outmask & (REDUCEDLENGTH | GEODESICSCALE)) > 0 + A2 = _A2m1f(eps) + _C2f(eps, C2a) + m0x = A1 - A2 + A2 = 1 + A2 + end + A1 = 1 + A1 + end + if (outmask & DISTANCE) > 0 + B1 = _SinCosSeries(true, ssig2, csig2, C1a) - _SinCosSeries(true, ssig1, csig1, C1a) + # Missing a factor of _b + s12b = A1 * (sig12 + B1) + if (outmask & (REDUCEDLENGTH | GEODESICSCALE)) > 0 + B2 = _SinCosSeries(true, ssig2, csig2, C2a) - _SinCosSeries(true, ssig1, csig1, C2a) + J12 = m0x * sig12 + (A1 * B1 - A2 * B2) + end + elseif (outmask & (REDUCEDLENGTH | GEODESICSCALE)) > 0 + # Assume here that nC1_ >= nC2_ + for l in 1:(nC2_ - 1) + C2a[l+1] = A1 * C1a[l+1] - A2 * C2a[l+1] + J12 = m0x * sig12 + _SinCosSeries(true, ssig2, csig2, C2a) - + _SinCosSeries(true, ssig1, csig1, C2a) + end + end + if (outmask & REDUCEDLENGTH) > 0 + m0 = m0x + # Missing a factor of _b. + # Add parens around (csig1 * ssig2) and (ssig1 * csig2) to ensure + # accurate cancellation in the case of coincident points. + m12b = (dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) - + csig1 * csig2 * J12) + end + if (outmask & GEODESICSCALE) > 0 + csig12 = csig1 * csig2 + ssig1 * ssig2 + t = self._ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2) + M12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1 + M21 = csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2 + end + return s12b, m12b, m0, M12, M21 +end + +"""Private: Find a starting value for Newton's method.""" +function _InverseStart(self, sbet1, cbet1, dn1, sbet2, cbet2, dn2, + lam12, slam12, clam12, + # Scratch areas of the right size + C1a, C2a) + # Return a starting point for Newton's method in salp1 and calp1 (function + # value is -1). If Newton's method doesn't need to be used, return also + # salp2 and calp2 and function value is sig12. + sig12 = -1.0 + salp2 = calp2 = dnm = Math.nan # Return values + # bet12 = bet2 - bet1 in [0, pi); bet12a = bet2 + bet1 in (-pi, 0] + sbet12 = sbet2 * cbet1 - cbet2 * sbet1 + cbet12 = cbet2 * cbet1 + sbet2 * sbet1 + # Volatile declaration needed to fix inverse cases + # 88.202499451857 0 -88.202499451857 179.981022032992859592 + # 89.262080389218 0 -89.262080389218 179.992207982775375662 + # 89.333123580033 0 -89.333123580032997687 179.99295812360148422 + # which otherwise fail with g++ 4.4.4 x86 -O3 + sbet12a = sbet2 * cbet1 + sbet12a += cbet2 * sbet1 + + shortline = cbet12 >= 0 && sbet12 < 0.5 && cbet2 * lam12 < 0.5 + if shortline + sbetm2 = Math.sq(sbet1 + sbet2) + # sin((bet1+bet2)/2)^2 + # = (sbet1 + sbet2)^2 / ((sbet1 + sbet2)^2 + (cbet1 + cbet2)^2) + sbetm2 /= sbetm2 + Math.sq(cbet1 + cbet2) + dnm = sqrt(1 + self._ep2 * sbetm2) + omg12 = lam12 / (self._f1 * dnm) + somg12 = sin(omg12) + comg12 = cos(omg12) + else + somg12 = slam12 + comg12 = clam12 + end + + salp1 = cbet2 * somg12 + calp1 = if comg12 >= 0 + sbet12 + cbet2 * sbet1 * Math.sq(somg12) / (1 + comg12) + else + sbet12a - cbet2 * sbet1 * Math.sq(somg12) / (1 - comg12) + end + + ssig12 = hypot(salp1, calp1) + csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12 + + if shortline && ssig12 < self._etol2 + # really short lines + salp2 = cbet1 * somg12 + calp2 = sbet12 - cbet1 * sbet2 * (comg12 >= 0 ? + Math.sq(somg12) / (1 + comg12) : + 1 - comg12) + salp2, calp2 = Math.norm(salp2, calp2) + # Set return value + sig12 = atan(ssig12, csig12) + elseif (abs(self._n) >= 0.1 || # Skip astroid calc if too eccentric + csig12 >= 0 || + ssig12 >= 6 * abs(self._n) * pi * Math.sq(cbet1)) + # Nothing to do, zeroth order spherical approximation is OK + else + # Scale lam12 and bet2 to x, y coordinate system where antipodal point + # is at origin and singular point is at y = 0, x = -1. + # real y, lamscale, betscale + # Volatile declaration needed to fix inverse case + # 56.320923501171 0 -56.320923501171 179.664747671772880215 + # which otherwise fails with g++ 4.4.4 x86 -O3 + # volatile real x + lam12x = atan(-slam12, -clam12) + if self.f >= 0 # In fact f == 0 does not get here + # x = dlong, y = dlat + k2 = Math.sq(sbet1) * self._ep2 + eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2) + lamscale = self.f * cbet1 * _A3f(self, eps) * pi + betscale = lamscale * cbet1 + x = lam12x / lamscale + y = sbet12a / betscale + else # _f < 0 + # x = dlat, y = dlong + cbet12a = cbet2 * cbet1 - sbet2 * sbet1 + bet12a = atan(sbet12a, cbet12a) + # real m12b, m0, dummy + # In the case of lon12 = 180, this repeats a calculation made in + # Inverse. + dummy, m12b, m0, dummy, dummy = _Lengths(self, + self._n, pi + bet12a, sbet1, -cbet1, dn1, sbet2, cbet2, dn2, + cbet1, cbet2, REDUCEDLENGTH, C1a, C2a) + x = -1 + m12b / (cbet1 * cbet2 * m0 * pi) + betscale = (x < -0.01 ? sbet12a / x : + -self.f * Math.sq(cbet1) * pi) + lamscale = betscale / cbet1 + y = lam12x / lamscale + end + + if y > -tol1_ && x > -1 - xthresh_ + # strip near cut + if self.f >= 0 + salp1 = min(1.0, -x) + calp1 = - sqrt(1 - Math.sq(salp1)) + else + calp1 = max((x > -tol1_ ? 0.0 : -1.0), x) + salp1 = sqrt(1 - Math.sq(calp1)) + end + else + # Estimate alp1, by solving the astroid problem. + # + # Could estimate alpha1 = theta + pi/2, directly, i.e., + # calp1 = y/k; salp1 = -x/(1+k); for _f >= 0 + # calp1 = x/(1+k); salp1 = -y/k; for _f < 0 (need to check) + # + # However, it's better to estimate omg12 from astroid and use + # spherical formula to compute alp1. This reduces the mean number of + # Newton iterations for astroid cases from 2.24 (min 0, max 6) to 2.12 + # (min 0 max 5). The changes in the number of iterations are as + # follows: + # + # change percent + # 1 5 + # 0 78 + # -1 16 + # -2 0.6 + # -3 0.04 + # -4 0.002 + # + # The histogram of iterations is (m = number of iterations estimating + # alp1 directly, n = number of iterations estimating via omg12, total + # number of trials = 148605): + # + # iter m n + # 0 148 186 + # 1 13046 13845 + # 2 93315 102225 + # 3 36189 32341 + # 4 5396 7 + # 5 455 1 + # 6 56 0 + # + # Because omg12 is near pi, estimate work with omg12a = pi - omg12 + k = _Astroid(x, y) + omg12a = lamscale * (self.f >= 0 ? -x * k/(1 + k) : + -y * (1 + k)/k ) + somg12 = sin(omg12a) + comg12 = -cos(omg12a) + # Update spherical estimate of alp1 using omg12 instead of lam12 + salp1 = cbet2 * somg12 + calp1 = sbet12a - cbet2 * sbet1 * Math.sq(somg12) / (1 - comg12) + end + end + # Sanity check on starting guess. Backwards check allows NaN through. + if !(salp1 <= 0) + salp1, calp1 = Math.norm(salp1, calp1) + else + salp1 = 1.0 + calp1 = 0.0 + end + sig12, salp1, calp1, salp2, calp2, dnm +end + +"""Private: Solve hybrid problem""" +function _Lambda12(self, sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1, + slam120, clam120, diffp, + # Scratch areas of the right size + C1a, C2a, C3a) + if sbet1 == 0 && calp1 == 0 + # Break degeneracy of equatorial line. This case has already been + # handled. + calp1 = -tiny_ + end + + # sin(alp1) * cos(bet1) = sin(alp0) + salp0 = salp1 * cbet1 + calp0 = hypot(calp1, salp1 * sbet1) # calp0 > 0 + + # real somg1, comg1, somg2, comg2, lam12 + # tan(bet1) = tan(sig1) * cos(alp1) + # tan(omg1) = sin(alp0) * tan(sig1) = tan(omg1)=tan(alp1)*sin(bet1) + ssig1 = sbet1; somg1 = salp0 * sbet1 + csig1 = comg1 = calp1 * cbet1 + ssig1, csig1 = Math.norm(ssig1, csig1) + # Math.norm(somg1, comg1); -- don't need to normalize! + + # Enforce symmetries in the case abs(bet2) = -bet1. Need to be careful + # about this case, since this can yield singularities in the Newton + # iteration. + # sin(alp2) * cos(bet2) = sin(alp0) + salp2 = cbet2 != cbet1 ? salp0 / cbet2 : salp1 + # calp2 = sqrt(1 - sq(salp2)) + # = sqrt(sq(calp0) - sq(sbet2)) / cbet2 + # and subst for calp0 and rearrange to give (choose positive sqrt + # to give alp2 in [0, pi/2]). + calp2 = if cbet2 != cbet1 || abs(sbet2) != -sbet1 + sqrt(Math.sq(calp1 * cbet1) + + (cbet1 < -sbet1 ? + (cbet2 - cbet1) * (cbet1 + cbet2) : + (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2 + else + abs(calp1) + end + # tan(bet2) = tan(sig2) * cos(alp2) + # tan(omg2) = sin(alp0) * tan(sig2). + ssig2 = sbet2; somg2 = salp0 * sbet2 + csig2 = comg2 = calp2 * cbet2 + ssig2, csig2 = Math.norm(ssig2, csig2) + # Math.norm(somg2, comg2); -- don't need to normalize! + + # sig12 = sig2 - sig1, limit to [0, pi] + sig12 = atan(max(0.0, csig1 * ssig2 - ssig1 * csig2), + csig1 * csig2 + ssig1 * ssig2) + + # omg12 = omg2 - omg1, limit to [0, pi] + somg12 = max(0.0, comg1 * somg2 - somg1 * comg2) + comg12 = comg1 * comg2 + somg1 * somg2 + # eta = omg12 - lam120 + eta = atan(somg12 * clam120 - comg12 * slam120, + comg12 * clam120 + somg12 * slam120) + + # real B312 + k2 = Math.sq(calp0) * self._ep2 + eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2) + _C3f(self, eps, C3a) + B312 = (_SinCosSeries(true, ssig2, csig2, C3a) - + _SinCosSeries(true, ssig1, csig1, C3a)) + domg12 = -self.f * _A3f(self, eps) * salp0 * (sig12 + B312) + lam12 = eta + domg12 + + if diffp + if calp2 == 0 + dlam12 = - 2 * self._f1 * dn1 / sbet1 + else + dummy, dlam12, dummy, dummy, dummy = _Lengths(self, + eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2, + REDUCEDLENGTH, C1a, C2a) + dlam12 *= self._f1 / (calp2 * cbet2) + end + else + dlam12 = Math.nan + end + + lam12, salp2, calp2, sig12, ssig1, csig1, ssig2, csig2, eps, domg12, dlam12 +end + +"""Private: General version of the inverse problem""" +function _GenInverse(self, lat1::T1, lon1::T2, lat2::T3, lon2::T4, outmask) where {T1,T2,T3,T4} + T = float(promote_type(Float64, T1, T2, T3, T4)) + lat1, lon1, lat2, lon2 = T.((lat1, lon1, lat2, lon2)) + a12 = s12 = m12 = M12 = M21 = S12 = Math.nan # return vals + + outmask &= OUT_MASK + # Compute longitude difference (AngDiff does this carefully). Result is + # in [-180, 180] but -180 is only for west-going geodesics. 180 is for + # east-going and meridional geodesics. + lon12, lon12s = Math.AngDiff(lon1, lon2) + # Make longitude difference positive. + lonsign = lon12 >= 0 ? 1 : -1 + # If very close to being on the same half-meridian, then make it so. + lon12 = lonsign * Math.AngRound(lon12) + lon12s = Math.AngRound((180 - lon12) - lonsign * lon12s) + lam12 = deg2rad(lon12) + if lon12 > 90 + slam12, clam12 = Math.sincosd(lon12s) + clam12 = -clam12 + else + slam12, clam12 = Math.sincosd(lon12) + end + + # If really close to the equator, treat as on equator. + lat1 = Math.AngRound(Math.LatFix(lat1)) + lat2 = Math.AngRound(Math.LatFix(lat2)) + # Swap points so that point with higher (abs) latitude is point 1 + # If one latitude is a nan, then it becomes lat1. + swapp = abs(lat1) < abs(lat2) ? -1 : 1 + if swapp < 0 + lonsign *= -1 + lat2, lat1 = lat1, lat2 + end + # Make lat1 <= 0 + latsign = lat1 < 0 ? 1 : -1 + lat1 *= latsign + lat2 *= latsign + # Now we have + # + # 0 <= lon12 <= 180 + # -90 <= lat1 <= 0 + # lat1 <= lat2 <= -lat1 + # + # longsign, swapp, latsign register the transformation to bring the + # coordinates to this canonical form. In all cases, 1 means no change was + # made. We make these transformations so that there are few cases to + # check, e.g., on verifying quadrants in atan2. In addition, this + # enforces some symmetries in the results returned. + + # real phi, sbet1, cbet1, sbet2, cbet2, s12x, m12x + + sbet1, cbet1 = Math.sincosd(lat1) + sbet1 *= self._f1 + # Ensure cbet1 = +epsilon at poles + sbet1, cbet1 = Math.norm(sbet1, cbet1) + cbet1 = max(tiny_, cbet1) + + sbet2, cbet2 = Math.sincosd(lat2); sbet2 *= self._f1 + # Ensure cbet2 = +epsilon at poles + sbet2, cbet2 = Math.norm(sbet2, cbet2) + cbet2 = max(tiny_, cbet2) + + # If cbet1 < -sbet1, then cbet2 - cbet1 is a sensitive measure of the + # |bet1| - |bet2|. Alternatively (cbet1 >= -sbet1), abs(sbet2) + sbet1 is + # a better measure. This logic is used in assigning calp2 in Lambda12. + # Sometimes these quantities vanish and in that case we force bet2 = +/- + # bet1 exactly. An example where is is necessary is the inverse problem + # 48.522876735459 0 -48.52287673545898293 179.599720456223079643 + # which failed with Visual Studio 10 (Release and Debug) + + if cbet1 < -sbet1 + if cbet2 == cbet1 + sbet2 = sbet2 < 0 ? sbet1 : -sbet1 + end + else + if abs(sbet2) == -sbet1 + cbet2 = cbet1 + end + end + + dn1 = sqrt(1 + self._ep2 * Math.sq(sbet1)) + dn2 = sqrt(1 + self._ep2 * Math.sq(sbet2)) + + # real a12, sig12, calp1, salp1, calp2, salp2 + # index zero elements of these arrays are unused + C1a = Vector{Float64}(undef, nC1_ + 1) + C2a = Vector{Float64}(undef, nC2_ + 1) + C3a = Vector{Float64}(undef, nC3_) + + meridian = lat1 == -90 || slam12 == 0 + + if meridian + + # Endpoints are on a single full meridian, so the geodesic might lie on + # a meridian. + + calp1 = clam12 + salp1 = slam12 # Head to the target longitude + calp2 = 1.0 + salp2 = 0.0 # At the target we're heading north + + # tan(bet) = tan(sig) * cos(alp) + ssig1 = sbet1; csig1 = calp1 * cbet1 + ssig2 = sbet2; csig2 = calp2 * cbet2 + + # sig12 = sig2 - sig1 + sig12 = atan(max(0.0, csig1 * ssig2 - ssig1 * csig2), + csig1 * csig2 + ssig1 * ssig2) + + s12x, m12x, dummy, M12, M21 = _Lengths(self, + self._n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2, + outmask | DISTANCE | REDUCEDLENGTH, C1a, C2a) + + # Add the check for sig12 since zero length geodesics might yield m12 < + # 0. Test case was + # + # echo 20.001 0 20.001 0 | GeodSolve -i + # + # In fact, we will have sig12 > pi/2 for meridional geodesic which is + # not a shortest path. + if sig12 < 1 || m12x >= 0 + if sig12 < 3 * tiny_ + sig12 = m12x = s12x = 0.0 + end + m12x *= self._b + s12x *= self._b + a12 = rad2deg(sig12) + else + # m12 < 0, i.e., prolate and too close to anti-podal + meridian = false + end + end + # end if meridian: + + # somg12 > 1 marks that it needs to be calculated + somg12 = 2.0 + comg12 = 0.0 + omg12 = 0.0 + if (! meridian && + sbet1 == 0 && # and sbet2 == 0 + # Mimic the way Lambda12 works with calp1 = 0 + (self.f <= 0 || lon12s >= self.f * 180)) + + # Geodesic runs along equator + calp1 = calp2 = 0.0 + salp1 = salp2 = 1.0 + s12x = self.a * lam12 + sig12 = omg12 = lam12 / self._f1 + m12x = self._b * sin(sig12) + if (outmask & GEODESICSCALE) > 0 + M12 = M21 = cos(sig12) + end + a12 = lon12 / self._f1 + + elseif ! meridian + + # Now point1 and point2 belong within a hemisphere bounded by a + # meridian and geodesic is neither meridional or equatorial. + + # Figure a starting point for Newton's method + sig12, salp1, calp1, salp2, calp2, dnm = _InverseStart(self, + sbet1, cbet1, dn1, sbet2, cbet2, dn2, lam12, slam12, clam12, C1a, C2a) + + if sig12 >= 0 + # Short lines (InverseStart sets salp2, calp2, dnm) + s12x = sig12 * self._b * dnm + m12x = (Math.sq(dnm) * self._b * sin(sig12 / dnm)) + if (outmask & GEODESICSCALE) > 0 + M12 = M21 = cos(sig12 / dnm) + end + a12 = rad2deg(sig12) + omg12 = lam12 / (self._f1 * dnm) + else + + # Newton's method. This is a straightforward solution of f(alp1) = + # lambda12(alp1) - lam12 = 0 with one wrinkle. f(alp) has exactly one + # root in the interval (0, pi) and its derivative is positive at the + # root. Thus f(alp) is positive for alp > alp1 and negative for alp < + # alp1. During the course of the iteration, a range (alp1a, alp1b) is + # maintained which brackets the root and with each evaluation of f(alp) + # the range is shrunk if possible. Newton's method is restarted + # whenever the derivative of f is negative (because the new value of + # alp1 is then further from the solution) or if the new estimate of + # alp1 lies outside (0,pi); in this case, the new starting guess is + # taken to be (alp1a + alp1b) / 2. + # real ssig1, csig1, ssig2, csig2, eps + numit = 0 + tripn = tripb = false + # Bracketing range + salp1a = tiny_ + calp1a = 1.0 + salp1b = tiny_ + calp1b = -1.0 + + while numit < maxit2_ + # the WGS84 test set: mean = 1.47, sd = 1.25, max = 16 + # WGS84 and random input: mean = 2.85, sd = 0.60 + (v, salp2, calp2, sig12, ssig1, csig1, ssig2, csig2, + eps, domg12, dv) = _Lambda12(self, + sbet1, cbet1, dn1, sbet2, cbet2, dn2, + salp1, calp1, slam12, clam12, numit < maxit1_, + C1a, C2a, C3a) + # 2 * tol0 is approximately 1 ulp for a number in [0, pi]. + # Reversed test to allow escape with NaNs + if tripb || !(abs(v) >= (tripn ? 8 : 1) * tol0_) + break + end + # Update bracketing values + if v > 0 && (numit > maxit1_ || + calp1/salp1 > calp1b/salp1b) + salp1b = salp1; calp1b = calp1 + elseif v < 0 && (numit > maxit1_ || + calp1/salp1 < calp1a/salp1a) + salp1a = salp1; calp1a = calp1 + end + + numit += 1 + if numit < maxit1_ && dv > 0 + dalp1 = -v/dv + sdalp1 = sin(dalp1) + cdalp1 = cos(dalp1) + nsalp1 = salp1 * cdalp1 + calp1 * sdalp1 + if nsalp1 > 0 && abs(dalp1) < pi + calp1 = calp1 * cdalp1 - salp1 * sdalp1 + salp1 = nsalp1 + salp1, calp1 = Math.norm(salp1, calp1) + # In some regimes we don't get quadratic convergence because + # slope -> 0. So use convergence conditions based on epsilon + # instead of sqrt(epsilon). + tripn = abs(v) <= 16 * tol0_ + continue + end + end + # Either dv was not positive or updated value was outside + # legal range. Use the midpoint of the bracket as the next + # estimate. This mechanism is not needed for the WGS84 + # ellipsoid, but it does catch problems with more eccentric + # ellipsoids. Its efficacy is such for + # the WGS84 test set with the starting guess set to alp1 = 90deg: + # the WGS84 test set: mean = 5.21, sd = 3.93, max = 24 + # WGS84 and random input: mean = 4.74, sd = 0.99 + salp1 = (salp1a + salp1b)/2 + calp1 = (calp1a + calp1b)/2 + salp1, calp1 = Math.norm(salp1, calp1) + tripn = false + tripb = (abs(salp1a - salp1) + (calp1a - calp1) < tolb_ || + abs(salp1 - salp1b) + (calp1 - calp1b) < tolb_) + end + + lengthmask = outmask + if (outmask & (REDUCEDLENGTH | GEODESICSCALE)) > 0 + lengthmask |= DISTANCE + else + lengthmask |= EMPTY + end + s12x, m12x, dummy, M12, M21 = _Lengths(self, + eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2, + lengthmask, C1a, C2a) + + m12x *= self._b + s12x *= self._b + a12 = rad2deg(sig12) + if (outmask & AREA) > 0 + # omg12 = lam12 - domg12 + sdomg12 = sin(domg12) + cdomg12 = cos(domg12) + somg12 = slam12 * cdomg12 - clam12 * sdomg12 + comg12 = clam12 * cdomg12 + slam12 * sdomg12 + end + end + end + # end elif not meridian + + if (outmask & DISTANCE) > 0 + s12 = 0.0 + s12x # Convert -0 to 0 + end + + if (outmask & REDUCEDLENGTH) > 0 + m12 = 0.0 + m12x # Convert -0 to 0 + end + + if (outmask & AREA) > 0 + # From Lambda12: sin(alp1) * cos(bet1) = sin(alp0) + salp0 = salp1 * cbet1 + calp0 = hypot(calp1, salp1 * sbet1) # calp0 > 0 + # real alp12 + if calp0 != 0 && salp0 != 0 + # From Lambda12: tan(bet) = tan(sig) * cos(alp) + ssig1 = sbet1; csig1 = calp1 * cbet1 + ssig2 = sbet2; csig2 = calp2 * cbet2 + k2 = Math.sq(calp0) * self._ep2 + eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2) + # Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0). + A4 = Math.sq(self.a) * calp0 * salp0 * self._e2 + ssig1, csig1 = Math.norm(ssig1, csig1) + ssig2, csig2 = Math.norm(ssig2, csig2) + C4a = Vector{Float64}(undef, nC4_) + _C4f(self, eps, C4a) + B41 = _SinCosSeries(false, ssig1, csig1, C4a) + B42 = _SinCosSeries(false, ssig2, csig2, C4a) + S12 = A4 * (B42 - B41) + else + # Avoid problems with indeterminate sig1, sig2 on equator + S12 = 0.0 + end + + if ! meridian && somg12 > 1 + somg12 = sin(omg12) + comg12 = cos(omg12) + end + + if (!meridian && + # omg12 < 3/4 * pi + comg12 > -0.7071 && # Long difference not too big + sbet2 - sbet1 < 1.75) # Lat difference not too big + # Use tan(Gamma/2) = tan(omg12/2) + # * (tan(bet1/2)+tan(bet2/2))/(1+tan(bet1/2)*tan(bet2/2)) + # with tan(x/2) = sin(x)/(1+cos(x)) + domg12 = 1 + comg12; dbet1 = 1 + cbet1; dbet2 = 1 + cbet2 + alp12 = 2 * atan( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ), + domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) ) + else + # alp12 = alp2 - alp1, used in atan2 so no need to normalize + salp12 = salp2 * calp1 - calp2 * salp1 + calp12 = calp2 * calp1 + salp2 * salp1 + # The right thing appears to happen if alp1 = +/-180 and alp2 = 0, viz + # salp12 = -0 and alp12 = -180. However this depends on the sign + # being attached to 0 correctly. The following ensures the correct + # behavior. + if salp12 == 0 && calp12 < 0 + salp12 = tiny_ * calp1 + calp12 = -1.0 + end + alp12 = atan(salp12, calp12) + end + S12 += self._c2 * alp12 + S12 *= swapp * lonsign * latsign + # Convert -0 to 0 + S12 += 0.0 + end + + # Convert calp, salp to azimuth accounting for lonsign, swapp, latsign. + if swapp < 0 + salp2, salp1 = salp1, salp2 + calp2, calp1 = calp1, calp2 + if (outmask & GEODESICSCALE) > 0 + M21, M12 = M12, M21 + end + end + + salp1 *= swapp * lonsign + calp1 *= swapp * latsign + salp2 *= swapp * lonsign + calp2 *= swapp * latsign + + a12, s12, salp1, calp1, salp2, calp2, m12, M12, M21, S12 +end + +end # module diff --git a/src/GeographicLib/GeographicLib.jl b/src/GeographicLib/GeographicLib.jl new file mode 100644 index 0000000..ddf3f02 --- /dev/null +++ b/src/GeographicLib/GeographicLib.jl @@ -0,0 +1,348 @@ +""" +# GeographicLib + +Compute distances and angles of geodesics (shortest paths) on a flattened sphere. + +See the online documentation at https://anowacki.github.io/GeographicLib.jl/stable +for full usage instructions. + +## Overview + +Users are recommended to use the Julia-style interface: + +### Constructors +- `Geodesic` creates an instance of an ellipsoid and associated parameters for later computation. +- `GeodesicLine` creates a great circle, including information about the ellipsoid, which allows + for multiple great circle calculations to be performed efficiently. +- `Polygon` creates a polygon made of great circles. + +### `forward` and `forward_deg` +- `forward` computes the end point given a starting point, azimuth and distance in + m (or whichever units the ellipsoid is defined in); +- `forward_deg` is the same but the distance from the starting point is defined + by an angular distance in °. + +Both of these functions can use either a `Geodesic` or `GeodesicLine` for the calculations. + +### `inverse` +- `inverse` computes the distance and forward and backazimuths between two points + on the flattened sphere. + +### `waypoints` +- `waypoints` will compute a set of evenly-spaced points along a `GeodesicLine`. + +### `Polygon`s +- `add_point!` will add a new vertex to a `Polygon` +- `add_edge!` adds a new vertex to a `Polygon` defined by an azimuth and + distance from the last vertex in the polygon + + +## Underlying library +One can also access the implementation of the original library, which has been +literally transcribed into Julia. These are exposed via: + +- `GeographicLib.AddPoint()` +- `GeographicLib.AddEdge()` +- `GeographicLib.ArcDirect()` +- `GeographicLib.ArcPosition()` +- `GeographicLib.Direct()` +- `GeographicLib.Compute()` +- `GeographicLib.Inverse()` +- `GeographicLib.SetArc()` +- `GeographicLib.SetDistance()` + +These operate as the original Python library `geographiclib`, except the first +argument to each is an ellipsoid as a `GeographicLib.Geodesics.Geodesic` object. +These can be constructed via the exported `Geodesic()` constructor. + +## Attribution + +GeographicLib.jl is a literal Julia transcription of the Python implementation of +the GeographicLib library by Charles F. F. Karney. + +The algorithms are derived in + + Charles F. F. Karney, + Algorithms for geodesics, J. Geodesy 87, 43-55 (2013), + https://doi.org/10.1007/s00190-012-0578-z + Addenda: https://geographiclib.sourceforge.io/geod-addenda.html + +The original library is licensed under the following terms: + +> Copyright (c) Charles Karney (2011-2017) and licensed +> under the MIT/X11 License. For more information, see +> https://geographiclib.sourceforge.io/ + +This Julia version is copyright (c) Andy Nowacki (2019) +and licensed under the MIT License. See the file `LICENSE.md` for details. +""" +module GeographicLib + +include("results.jl") +include("Constants.jl") +include("Math.jl") +include("GeodesicCapability.jl") +include("Geodesics.jl") +include("GeodesicLines.jl") +include("Accumulators.jl") + +import .Constants, .Math, .GeodesicCapability, .Geodesics, .GeodesicLines, + .Accumulators + +"Ellipsoid of the WGS84 system, with a semimajor radius of $(Constants.WGS84_a) m +and flattening $(Constants.WGS84_f)." +const WGS84 = Geodesics.Geodesic(Constants.WGS84_a, Constants.WGS84_f) + +include("direct.jl") +include("Polygons.jl") + +import .Polygons + +using .Geodesics: Inverse, Geodesic +using .GeodesicLines: GeodesicLine, Position, ArcPosition, SetDistance, SetArc +using .Polygons: Polygon + +include("inverse_line.jl") +include("show.jl") + +export Geodesic, GeodesicLine, Polygon, + add_edge!, add_point!, forward, forward_deg, inverse, properties, waypoints + + +""" + forward([ellipsoid::Geodesic=WGS84,] lon, lat, azi, dist) -> lon′, lat′, baz, dist, angle + +Compute the final position when travelling from longitude `lon`°, latitude `lat`°, +along an azimuth of `azi`° for `dist` m (or whichever units the `ellipsoid` +is defined using). The final coordinates are (`lon′`, `lat′`)°, +the backazimuth from the final point to the start is `baz`°, and the angular distance is +`angle`°. + +If `ellipsoid` is not supplied, then [`WGS84`](@ref) is used by default. +""" +function forward(g::Geodesic, lon, lat, azi, dist) + r = Direct(g, lat, lon, azi, dist, Geodesics.ALL) + (lon=r.lon2, lat=r.lat2, baz=Math.AngNormalize(r.azi2 + 180), dist=dist, angle=r.a12) +end + +forward(lon, lat, azi, dist) = forward(WGS84, lon, lat, azi, dist) + +""" + forward(a, f, lon, lat, azi, dist) -> lon′, lat′, baz, dist, angle + +Compute the final point assuming an ellipsoid with semimajor axis `a` m +and flattening `f`. + +Note that precomputing the ellipsoid with [`Geodesic`](@ref) and then reusing this +if multiple [`forward`](@ref) or [`forward_deg`](@ref) calculations are +needed will be more efficient. +""" +forward(a, f, lon, lat, azi, dist) = forward(Geodesic(a, f), lon, lat, azi, dist) + +""" + forward(line::GeodesicLine, dist) -> lon′, lat′, baz, dist, angle + +Compute the final point when travelling along a pre-computed great circle `line` for +a distance of `dist` m. `angle` is the angular distance in °. +""" +function forward(line::GeodesicLine, dist) + r = GeodesicLines.Position(line, dist) + (lon=r.lon2, lat=r.lat2, baz=Math.AngNormalize(r.azi2 + 180), dist=dist, angle=r.a12) +end + + +""" + forward_deg([ellipsoid::Geodesic=WGS84,] lon, lat, azi, angle) -> lon′, lat′, baz, dist, angle + +Compute the final position when travelling from longitude `lon`°, latitude `lat`°, +along an azimuth of `azi`° for an angular distance of `angle`°. +The final coordinates are (`lon′`, `lat′`)°, +the backazimuth from the final point to the start is `baz`°, and the distance is +`dist` m. + +If `ellipsoid` is not supplied, then [`WGS84`](@ref) is used by default. +""" +function forward_deg(g::Geodesic, lon, lat, azi, angle) + r = ArcDirect(g, lat, lon, azi, angle, Geodesics.ALL) + (lon=r.lon2, lat=r.lat2, baz=Math.AngNormalize(r.azi2 + 180), dist=r.s12, angle=r.a12) +end + +forward_deg(lon, lat, azi, angle) = forward_deg(WGS84, lon, lat, azi, angle) + +""" + + forward_deg(a, f, lon, lat, azi, angle) -> lon′, lat′, baz, dist, angle + +Compute the final point assuming an ellipsoid with semimajor radius `a` m +and flattening `f`. + +Note that precomputing the ellipsoid with [`Geodesic`](@ref) and then reusing this +if multiple `forward_deg` calculations are needed will be more efficient. +""" +forward_deg(a, f, lon, lat, azi, angle) = forward_deg(Geodesic(a, f), lon, lat, azi, angle) + +""" + forward_deg(line::GeodesicLine, angle) -> lon, lat, baz, dist + +Compute the final point when travelling along a pre-computed great circle `line` for +an angular distance of `angle` °. `dist` is the distance in m. +""" +function forward_deg(line::GeodesicLine, angle) + r = GeodesicLines.ArcPosition(line, angle) + (lon=r.lon2, lat=r.lat2, baz=Math.AngNormalize(r.azi2 + 180), dist=r.s12, angle=r.a12) +end + +""" + inverse([ellipsoid::Geodesic=WGS84,] lon1, lat1, lon2, lat2) -> azi, baz, dist, angle + +Compute for forward azimuth `azi`°, backazimuth `baz`°, surface distance `dist` m +and angular distance `angle`° when travelling from (`lon1`, `lat1`)° to a +second point (`lon2`, `lat2`)°. + +If `ellipsoid` is not supplied, then [`WGS84`](@ref) is used by default. +""" +function inverse(g::Geodesic, lon1, lat1, lon2, lat2) + r = Inverse(g, lat1, lon1, lat2, lon2) + (azi=r.azi1, baz=Math.AngNormalize(r.azi2 + 180), dist=r.s12, angle=r.a12) +end + +inverse(lon1, lat1, lon2, lat2) = inverse(WGS84, lon1, lat1, lon2, lat2) + +""" + inverse(a, f, lon1, lat1, lon2, lat2) -> azi, baz, dist, angle + +Compute the final point assuming an ellipsoid with semimajor radius `a` m +and flattening `f`. + +Note that precomputing the ellipsoid with `Geodesic(a, f)` and then reusing this +if multiple `inverse` calculations are needed will be more efficient. +""" +inverse(a, f, lon1, lat1, lon2, lat2) = inverse(Geodesic(a, f), lon1, lat1, lon2, lat2) + +""" + GeodesicLine([ellipsoid::Geodesic.WGS84,] lon1, lat1; azi, lon2, lat2, angle, dist) + +Construct a `GeodesicLine`, which may be used to efficiently compute many distances +along a great circle. Set the coordinates of the starting point `lon1`° and `lat1`°. +There are two ways to define the great circle this way: + +1. Set a start point, azimuth, and optionally distance. This requires the keyword + arguments `azi1` (all °) and optionally either `angle` (angular distance, °) + or distance `dist` (m). +2. Set a start point and end point. This requires the keyword arguments + `lon2` and `lat2` (all °). + +If `ellipsoid` is not supplied, then [`WGS84`](@ref) is used by default. + +See [`forward`](@ref) and [`forward_deg`](@ref) for details of computing points along a `GeodesicLine`. +""" +function GeodesicLine(geod::Geodesic, lon1, lat1; azi=nothing, lon2=nothing, lat2=nothing, + angle=nothing, dist=nothing) + if azi === nothing && any(x->x===nothing, (lon2, lat2)) + throw(ArgumentError("cannot define a GeodesicLine without either azi or (lon2, lat2)")) + end + if azi !== nothing + (lon2 !== nothing || lat2 !== nothing) && + throw(ArgumentError("cannot define both azimuth and end position")) + (angle !== nothing && dist !== nothing) && + throw(ArgumentError("cannot define both angle and dist")) + if angle !== nothing + ArcDirectLine(geod, lat1, lon1, azi, angle, Geodesics.ALL) + elseif dist !== nothing + DirectLine(geod, lat1, lon1, azi, dist, Geodesics.ALL) + else + GeodesicLine(geod, lat1, lon1, azi, caps=Geodesics.ALL) + end + elseif lon2 !== nothing && lat2 !== nothing + (angle === nothing && dist == nothing) || + throw(ArgumentError("cannot define a dist as well as an end point")) + InverseLine(geod, lat1, lon1, lat2, lon2, Geodesics.ALL) + else + throw(ArgumentError("invalid combination of keyword arguments")) + end +end + +GeodesicLine(lon1, lat1; kwargs...) = GeodesicLine(WGS84, lon1, lat1; kwargs...) + +""" + waypoints(line::GeodesicLine; n, dist, angle) -> points + +Return a set of `points` along a great circle `line` (which must have been specified with a +distance or as between two points). There are three options: + +- Specify `n`, the number of points (including the endpoints) +- Specift `dist`, the distance between each point in m +- Specify `angle`, the angular distance between each point in °. + +In all cases, the start and end points are always included. When giving either `dist` or `angle`, +the penultimate point may not be respectively `dist` m or `angle`° away from the final point. + +The output is a vector of named tuples as returned by [`forward`](@ref) +or [`forward_deg`](@ref). +""" +function waypoints(line::GeodesicLine; n=nothing, dist=nothing, angle=nothing) + isnan(line.a13) && + throw(ArgumentError("cannot compute waypoints for a line constructed without a distance or endpoint")) + count(x->x!==nothing, (n, dist, angle)) == 1 || + throw(ArgumentError("only one keyword parameter can be given")) + if n !== nothing + n >= 2 || throw(ArgumentError("number of points must be 2 or more")) + forward.(line, range(0, stop=line.s13, length=n)) + elseif dist !== nothing + n = ceil(Int, line.s13/dist) + [forward(line, min(dist*i, line.s13)) for i in 0:n] + elseif angle !== nothing + n = ceil(Int, line.a13/angle) + [forward_deg(line, min(angle*i, line.a13)) for i in 0:n] + end +end + +""" + Polygon([ellipsoid::Geodesic=WGS84,] lons, lat) -> polygon + +Construct a `polygon` from arrays of coordinates `lons` and `lats`, +optionally specifying an `ellipsoid`, which defaults to [`WGS84`](@ref). + +Calculate the polygon's area and perimeter with [`properties`](@ref). +""" +function Polygon(ellipsoid::Geodesic, lons::AbstractArray, lats::AbstractArray) + axes(lons) == axes(lats) || + throw(ArgumentError("lon and lat arrays must be the same size")) + polygon = Polygon() + for (lon, lat) in zip(lons, lats) + add_point!(polygon, lon, lat) + end + polygon +end + +Polygon(lons::AbstractArray, lats::AbstractArray) = Polygon(WGS84, lons, lats) + +""" + add_point!(polygon, lon, lat) -> polygon + +Add a point to a geodesic polygon in order to later compute its perimeter and +area with [`properties`](@ref). +""" +add_point!(polygon::Polygon, lon, lat) = Polygons.AddPoint!(polygon, lat, lon) + +""" + add_edge!(polygon, azi, dist) -> polygon + +Add a point to a polygon by specifying the azimuth `azi` (°) and +distance `dist` (m) from the previous point, in order to later compute its +perimeter and area with [`properties`](@ref). +""" +add_edge!(polygon::Polygon, azi, dist) = Polygons.AddEdge!(polygon, azi, dist) + +""" + properties(polygon::Polygon) -> npoints, perimeter, area + +Return the number of points `npoints`, `perimeter` (m) and `area` (m²) +of the geodesic `polygon`. +""" +function properties(polygon::Polygon) + n, perimeter, area = Polygons.Compute(polygon) + (n=n, perimeter=perimeter, area=area) +end + +end # module diff --git a/src/GeographicLib/Math.jl b/src/GeographicLib/Math.jl new file mode 100644 index 0000000..e5a8b59 --- /dev/null +++ b/src/GeographicLib/Math.jl @@ -0,0 +1,144 @@ +module Math + +const digits = 53 +const epsilon = eps() +const minval = 2.0^-1022 +const maxval = 2.0^1023 * (2 - epsilon) +const inf = Inf +const nan = NaN + +"""Square a number""" +sq(x) = x*x + +"""Real cube root of a number""" +function cbrt(x) + y = abs(x)^(1/3.0) + x >= 0 ? y : -y +end + +"""log(1 + x) accurate for small x (missing from python 2.5.2)""" +log1p(x) = Base.log1p(x) + +"""atanh(x) (missing from python 2.5.2)""" +function atanh(x) + y = abs(x) + y = log1p(2 * y/(1 - y))/2 + x < 0 ? -y : y +end + +"""return x with the sign of y (missing from python 2.5.2)""" +copysign(x, y) = Base.copysign(x, y) + +"""Private: Normalize a two-vector""" +norm(x, y) = (r = hypot(x, y); (x/r, y/r)) + +"""Error free transformation of a sum. +Note that t can be the same as one of the first two arguments.""" +function sum(u, v) + s = u + v + up = s - v + vpp = s - up + up -= u + vpp -= v + t = -(up + vpp) + s, t +end + +"""Evaluate a polynomial""" +function polyval(N, p, s, x) + y = float(N < 0 ? 0 : p[s+1]) + while N > 0 + N -= 1 + s += 1 + y = y*x + p[s+1] + end + y +end + +""" +Private: Round an angle so that small values underflow to zero. +""" +function AngRound(x) + # The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57 + # for reals = 0.7 pm on the earth if x is an angle in degrees. (This + # is about 1000 times more resolution than we get with angles around 90 + # degrees.) We use this to avoid having to deal with near singular + # cases when x is non-zero but tiny (e.g., 1.0e-200). + z = 1/16.0 + y = abs(x) + # The compiler mustn't "simplify" z - (z - y) to y + y < z && (y = z - (z - y)) + x < 0 ? zero(x) - y : y +end + +"""reduce angle in (-180,180]""" +function AngNormalize(x) + y = x % 360 + y = x == 0 ? x : y + if y <= -180 + y + 360 + elseif y <= 180 + y + else + y - 360 + end +end + +"""replace angles outside [-90,90] by NaN""" +LatFix(x) = abs(x) > 90 ? nan : x + +"Compute y - x and reduce to [-180, 180]° accurately" +function AngDiff(x, y) + d, t = sum(AngNormalize(-x), AngNormalize(y)) + d = AngNormalize(d) + sum(d == 180 && t > 0 ? -180 : d, t) +end + +"""Compute sine and cosine of x in degrees.""" +function sincosd(x) + r = x % 360 + q = isnan(r) ? nan : floor(Int, r / 90 + 0.5) + r -= 90 * q + r = deg2rad(r) + s = sin(r) + c = cos(r) + q = mod(q, 4) + if q == 1 + s, c = c, -s + elseif q == 2 + s, c = -s, -c + elseif q == 3 + s, c = -c, s + end + s, c = x == 0 ? (x, c) : (0.0+s, 0.0+c) +end + +"""compute atan2(y, x) with the result in degrees""" +function atan2d(y::T1, x::T2) where {T1,T2} + T = float(promote_type(T1, T2)) + if abs(y) > abs(x) + q = 2 + x, y = y, x + else + q = 0 + end + if x < 0 + q += 1 + x = -x + end + ang = rad2deg(atan(y, x)) + if q == 1 + ang = (y >= 0 ? 180 : -180) - ang + elseif q == 2 + ang = 90 - ang + elseif q == 3 + ang = -90 + ang + end + T(ang) +end + +isfinite(x) = Base.isfinite(x) + +isnan(x) = Base.isnan(x) + +end # module diff --git a/src/GeographicLib/Polygons.jl b/src/GeographicLib/Polygons.jl new file mode 100644 index 0000000..e63bd25 --- /dev/null +++ b/src/GeographicLib/Polygons.jl @@ -0,0 +1,366 @@ +module Polygons + +import ..Math, ..Geodesics, ..Accumulators, .._GenDirect, ..WGS84 + +export Polygon, AddPoint!, AddEdge!, Compute + +#"""Area of a geodesic polygon""" +mutable struct Polygon + """The geodesic object (readonly)""" + earth::Geodesics.Geodesic + """Is this a polyline? (readonly)""" + polyline::Bool + """The total area of the ellipsoid in meter^2 (readonly)""" + area0::Float64 + _mask::Int + _areasum::Accumulators.Accumulator{Float64} + _perimetersum::Accumulators.Accumulator{Float64} + """The current number of points in the polygon (readonly)""" + num::Int + _lon0::Float64 + _lat0::Float64 + """The current latitude in degrees (readonly)""" + lat1::Float64 + """The current longitude in degrees (readonly)""" + lon1::Float64 + _crossings::Int +end + +Base.:(==)(p1::Polygon, p2::Polygon) = + all(x -> (isbits(x[1]) && isnan(x[1]) && isnan(x[2])) || x[1] == x[2], + (getfield(p1, f), getfield(p2, f)) for f in fieldnames(Polygon)) + +# Treat the Geodesic struct as a scalar +Base.Broadcast.broadcastable(p::Polygon) = Ref(p) + +""" + Polygon([ellipsoid::Geodesic=WGS84,] polyline=false) + +Construct a `Polygon`, which contains a set of points on a certain +`ellipsoid`. + +With this construction, the `Polygon` contains no points. + +If `polyline` is true, then the `Polygon` will not accumulate area +and instead only its perimeter can be calculated. +""" +function Polygon(ellipsoid::Geodesics.Geodesic, polyline::Bool=false) + area0 = 4π*ellipsoid._c2 + _mask = (Geodesics.LATITUDE | Geodesics.LONGITUDE | + Geodesics.DISTANCE | (polyline ? Geodesics.EMPTY : + Geodesics.AREA | Geodesics.LONG_UNROLL)) + _areasum = Accumulators.Accumulator() + _perimetersum = Accumulators.Accumulator() + num = 0 + _lat0 = _lon0 = lat1 = lon1 = Math.nan + _crossings = 0 + Polygon(ellipsoid, polyline, area0, _mask, _areasum, _perimetersum, + num, _lon0, _lat0, lat1, lon1, _crossings) +end + +Polygon(polyline::Bool=false) = Polygon(WGS84, polyline) + +"""Count crossings of prime meridian for AddPoint.""" +function _transit(lon1, lon2) + # Return 1 or -1 if crossing prime meridian in east or west direction. + # Otherwise return zero. + # Compute lon12 the same way as Geodesic::Inverse. + lon1 = Math.AngNormalize(lon1) + lon2 = Math.AngNormalize(lon2) + lon12, _ = Math.AngDiff(lon1, lon2) + cross = (lon1 <= 0 && lon2 > 0 && lon12 > 0 ? 1 : + (lon2 <= 0 && lon1 > 0 && lon12 < 0 ? -1 : 0)) + cross +end + +"""Count crossings of prime meridian for AddEdge.""" +function _transitdirect(lon1, lon2) + # We want to compute exactly + # int(floor(lon2 / 360)) - int(floor(lon1 / 360)) + # Since we only need the parity of the result we can use std::remquo but + # this is buggy with g++ 4.8.3 and requires C++11. So instead we do + lon1 = lon1 % 720.0 + lon2 = lon2 % 720.0 + (((lon2 >= 0 && lon2 < 360) || lon2 < -360) ? 0 : 1) - + (((lon1 >= 0 && lon1 < 360) || lon1 < -360) ? 0 : 1) +end + +"""Reset to empty polygon.""" +function Clear!(self) + self.num = 0 + self._crossings = 0 + if !self.polyline + Accumulators.Set!(self._areasum, 0) + end + Accumulators.Set!(self._perimetersum, 0) + self._lat0 = self._lon0 = self.lat1 = self.lon1 = Math.nan + self +end + +""" + AddPoint!(polygon, lat, lon) -> polygon + +Add the next vertex to the polygon + +`lat`: the latitude of the point in degrees +`lon`: the longitude of the point in degrees + +This adds an edge from the current vertex to the new vertex. +""" +function AddPoint!(self::Polygon, lat, lon) + if self.num == 0 + self._lat0 = self.lat1 = lat + self._lon0 = self.lon1 = lon + else + _, s12, _, _, _, _, _, _, _, S12 = Geodesics._GenInverse(self.earth, + self.lat1, self.lon1, lat, lon, self._mask) + Accumulators.Add!(self._perimetersum, s12) + if !self.polyline + Accumulators.Add!(self._areasum, S12) + self._crossings += _transit(self.lon1, lon) + end + self.lat1 = lat + self.lon1 = lon + end + self.num += 1 + self +end + +""" + AddEdge!(polygon, azimuth, distance) -> polygon + +Add the next edge to the polygon + +`azi`: the azimuth at the current the point in degrees +`s`: the length of the edge in meters + +This specifies the new vertex in terms of the edge from the current +vertex. +""" +function AddEdge!(self::Polygon, azi, s) + if self.num != 0 + _, lat, lon, _, _, _, _, _, S12 = _GenDirect(self.earth, + self.lat1, self.lon1, azi, false, s, self._mask) + Accumulators.Add!(self._perimetersum, s) + if !self.polyline + Accumulators.Add!(self._areasum, S12) + self._crossings += _transitdirect(self.lon1, lon) + end + self.lat1 = lat + self.lon1 = lon + self.num += 1 + else + error("cannot add an edge to a polygon with no points") + end + self +end + +# return number, perimeter, area +""" + Compute(polygon, reverse=false, sign=true) -> n, perimeter, area + +Compute the properties of the polygon + +- `reverse`: if true then clockwise (instead of + counter-clockwise) traversal counts as a positive area +- `sign`: if true then return a signed result for the area if the + polygon is traversed in the "wrong" direction instead of returning + the area for the rest of the ellispoid + +Return a tuple of number, perimeter (meters), area (meters^2). + +If the object is a polygon (and not a polygon), the perimeter +includes the length of a final edge connecting the current point to +the initial point. If the object is a polyline, then area is nan. + +More points can be added to the polygon after this call. + +""" +function Compute(self::Polygon, reverse = false, sign = true) + if self.polyline + area = Math.nan + end + if self.num < 2 + perimeter = 0.0 + if !self.polyline + area = 0.0 + end + return self.num, perimeter, area + end + + if self.polyline + perimeter = Accumulators.Sum(self._perimetersum) + return self.num, perimeter, area + end + + _, s12, _, _, _, _, _, _, _, S12 = Geodesics._GenInverse(self.earth, + self.lat1, self.lon1, self._lat0, self._lon0, self._mask) + perimeter = Accumulators.Sum(self._perimetersum, s12) + tempsum = Accumulators.Accumulator(self._areasum) + Accumulators.Add!(tempsum, S12) + crossings = self._crossings + _transit(self.lon1, self._lon0) + if (crossings & 1) > 0 + Accumulators.Add!(tempsum, (Accumulators.Sum(tempsum) < 0 ? 1 : -1) * self.area0/2) + end + # area is with the clockwise sense. If !reverse convert to + # counter-clockwise convention. + if !reverse + Accumulators.Negate!(tempsum) + end + # If sign put area in (-area0/2, area0/2], else put area in [0, area0) + if sign + if Accumulators.Sum(tempsum) > self.area0/2 + Accumulators.Add!(tempsum, -self.area0 ) + elseif Accumulators.Sum(tempsum) <= -self.area0/2 + Accumulators.Add!(tempsum, self.area0 ) + end + else + if Accumulators.Sum(tempsum) >= self.area0 + Accumulators.Add!(tempsum, -self.area0 ) + elseif Accumulators.Sum(tempsum) < 0 + Accumulators.Add!(tempsum, self.area0 ) + end + end + + area = 0.0 + Accumulators.Sum(tempsum) + self.num, perimeter, area +end + + +"""Compute the properties for a tentative additional vertex + +:param lat: the latitude of the point in degrees +:param lon: the longitude of the point in degrees +:param reverse: if true then clockwise (instead of + counter-clockwise) traversal counts as a positive area +:param sign: if true then return a signed result for the area if the + polygon is traversed in the "wrong" direction instead of returning + the area for the rest of the earth +:return: a tuple of number, perimeter (meters), area (meters^2) + +""" +function TestPoint(self::Polygon, lat, lon, reverse = false, sign = true) + if self.polyline + area = Math.nan + end + if self.num == 0 + perimeter = 0.0 + if !self.polyline + area = 0.0 + end + return 1, perimeter, area + end + + perimeter = Accumulators.Sum(self._perimetersum) + tempsum = self.polyline ? 0.0 : Accumulators.Sum(self._areasum) + crossings = self._crossings + num = self.num + 1 + for i in (self.polyline ? (0,) : (0, 1)) + _, s12, _, _, _, _, _, _, _, S12 = Geodesics._GenInverse(self.earth, + i == 0 ? self.lat1 : lat, + i == 0 ? self.lon1 : lon, + i != 0 ? self._lat0 : lat, + i != 0 ? self._lon0 : lon, + self._mask) + perimeter += s12 + if !self.polyline + tempsum += S12 + crossings += _transit(i == 0 ? self.lon1 : lon, + i != 0 ? self._lon0 : lon) + end + end + + if self.polyline + return num, perimeter, area + end + + if (crossings & 1) > 0 + tempsum += (tempsum < 0 ? 1 : -1) * self.area0/2 + end + # area is with the clockwise sense. If !reverse convert to + # counter-clockwise convention. + if !reverse + tempsum *= -1 + end + # If sign put area in (-area0/2, area0/2], else put area in [0, area0) + if sign + if tempsum > self.area0/2 + tempsum -= self.area0 + elseif tempsum <= -self.area0/2 + tempsum += self.area0 + end + else + if tempsum >= self.area0 + tempsum -= self.area0 + elseif tempsum < 0 + tempsum += self.area0 + end + end + + area = 0.0 + tempsum + num, perimeter, area +end + +# return num, perimeter, area +"""Compute the properties for a tentative additional edge + +:param azi: the azimuth at the current the point in degrees +:param s: the length of the edge in meters +:param reverse: if true then clockwise (instead of + counter-clockwise) traversal counts as a positive area +:param sign: if true then return a signed result for the area if the + polygon is traversed in the "wrong" direction instead of returning + the area for the rest of the earth +:return: a tuple of number, perimeter (meters), area (meters^2) + +""" +function TestEdge(self::Polygon, azi, s, reverse = False, sign = True) + if self.num == 0 # we don't have a starting point! + return 0, Math.nan, Math.nan + end + num = self.num + 1 + perimeter = self._perimetersum.Sum() + s + if self.polyline + return num, perimeter, Math.nan + end + + tempsum = Accumulators.Sum(self._areasum) + crossings = self._crossings + _, lat, lon, _, _, _, _, _, S12 = _GenDirect(self.earth, + self.lat1, self.lon1, azi, false, s, self._mask) + tempsum += S12 + crossings += _transitdirect(self.lon1, lon) + _, s12, _, _, _, _, _, _, _, S12 = Geodesics._GenInverse(self.earth, + lat, lon, self._lat0, self._lon0, self._mask) + perimeter += s12 + tempsum += S12 + crossings += _transit(lon, self._lon0) + + if (crossings & 1) > 0 + tempsum += (tempsum < 0 ? 1 : -1) * self.area0/2 + end + # area is with the clockwise sense. If !reverse convert to + # counter-clockwise convention. + if !reverse + tempsum *= -1 + end + # If sign put area in (-area0/2, area0/2], else put area in [0, area0) + if sign + if tempsum > self.area0/2 + tempsum -= self.area0 + elseif tempsum <= -self.area0/2 + tempsum += self.area0 + end + else + if tempsum >= self.area0 + tempsum -= self.area0 + elseif tempsum < 0 + tempsum += self.area0 + end + end + + area = 0.0 + tempsum + num, perimeter, area +end + +end # module diff --git a/src/GeographicLib/direct.jl b/src/GeographicLib/direct.jl new file mode 100644 index 0000000..0aa85b9 --- /dev/null +++ b/src/GeographicLib/direct.jl @@ -0,0 +1,130 @@ +""" + Direct(geodesic, lat1, lon1, azi1, s12, outmask=Geodesics.STANDARD) -> result::Geodesics.Result + +Solve the direct geodesic problem and return a `Geodesics.Result` +containing the parameters of interest + +Input parameters: +`lat1`: latitude of the first point in degrees +`lon1`: longitude of the first point in degrees +`azi1`: azimuth at the first point in degrees +`s12`: the distance from the first point to the second in meters +`outmask`: a mask setting which output values are computed (see note below) + +Compute geodesic starting at (`lat1`, `lon1`) with azimuth `azi1` +and length `s12`. The default value of `outmask` is STANDARD, i.e., +the `lat1`, `lon1`, `azi1`, `lat2`, `lon2`, `azi2`, `s12`, `a12` +entries are returned. + +### Output mask + +May be any combination of: +`Geodesics.EMPTY`, `Geodesics.LATITUDE`, `Geodesics.LONGITUDE`, +`Geodesics.AZIMUTH`, `Geodesics.DISTANCE`, `Geodesics.STANDARD`, +`Geodesics.DISTANCE_IN`, `Geodesics.REDUCEDLENGTH`, `Geodesics.GEODESICSCALE`, +`Geodesics.AREA`, `Geodesics.ALL` or `Geodesics.LONG_UNROLL`. +See the docstring for each for more information. + +Flags are combined by bitwise or-ing values together, e.g. +`Geodesics.AZIMUTH | Geodesics.DISTANCE`. +""" +function Direct(self::Geodesics.Geodesic, lat1::T1, lon1::T2, azi1::T3, s12::T4, + outmask = GeodesicCapability.STANDARD) where {T1,T2,T3,T4} + + T = float(promote_type(Float64, T1, T2, T3, T4)) + lat1, lon1, azi1, s12 = T.((lat1, lon1, azi1, s12)) + outmask = Int(outmask) + + a12, lat2, lon2, azi2, s12, m12, M12, M21, S12 = _GenDirect(self, + lat1, lon1, azi1, false, s12, outmask) + outmask &= Geodesics.OUT_MASK + result = Geodesics.Result() + result.lat1 = Math.LatFix(lat1) + result.lon1 = (outmask & Geodesics.LONG_UNROLL) > 0 ? + lon1 : + Math.AngNormalize(lon1) + result.azi1 = Math.AngNormalize(azi1) + result.s12 = s12 + result.a12 = a12 + (outmask & Geodesics.LATITUDE) > 0 && (result.lat2 = lat2) + (outmask & Geodesics.LONGITUDE) > 0 && (result.lon2 = lon2) + (outmask & Geodesics.AZIMUTH) > 0 && (result.azi2 = azi2) + (outmask & Geodesics.REDUCEDLENGTH) > 0 && (result.m12 = m12) + if (outmask & Geodesics.GEODESICSCALE) > 0 + result.M12 = M12 + result.M21 = M21 + end + (outmask & Geodesics.AREA) > 0 && (result.S12 = S12) + result +end + +""" + ArcDirect(geodesic, lat1, lon1, azi1, a12, outmask=Geodesics.STANDARD) -> result::Geodesics.Result + +Solve the direct geodesic problem and return a `Geodesics.Result` +containing the parameters of interest + +Input parameters: +`lat1`: latitude of the first point in degrees +`lon1`: longitude of the first point in degrees +`azi1`: azimuth at the first point in degrees +`a12`: the angular distance from the first point to the second in ° +`outmask`: a mask setting which output values are computed (see note below) + +Compute geodesic starting at (`lat1`, `lon1`)° with azimuth `azi1`° +and arc length `a12`°. The default value of `outmask` is STANDARD, i.e., +the `lat1`, `lon1`, `azi1`, `lat2`, `lon2`, `azi2`, `s12`, `a12` +entries are returned. + +### Output mask + +May be any combination of: +`Geodesics.EMPTY`, `Geodesics.LATITUDE`, `Geodesics.LONGITUDE`, +`Geodesics.AZIMUTH`, `Geodesics.DISTANCE`, `Geodesics.STANDARD`, +`Geodesics.DISTANCE_IN`, `Geodesics.REDUCEDLENGTH`, `Geodesics.GEODESICSCALE`, +`Geodesics.AREA`, `Geodesics.ALL` or `Geodesics.LONG_UNROLL`. +See the docstring for each for more information. + +Flags are combined by bitwise or-ing values together, e.g. +`Geodesics.AZIMUTH | Geodesics.DISTANCE`. +""" +function ArcDirect(self::Geodesics.Geodesic, lat1::T1, lon1::T2, azi1::T3, a12::T4, + outmask = GeodesicCapability.STANDARD) where {T1,T2,T3,T4} + + T = float(promote_type(Float64, T1, T2, T3, T4)) + lat1, lon1, azi1, a12 = T.((lat1, lon1, azi1, a12)) + outmask = Int(outmask) + + a12, lat2, lon2, azi2, s12, m12, M12, M21, S12 = _GenDirect(self, + lat1, lon1, azi1, true, a12, outmask) + outmask &= Geodesics.OUT_MASK + result = Geodesics.Result() + result.lat1 = Math.LatFix(lat1) + result.lon1 = (outmask & Geodesics.LONG_UNROLL) > 0 ? + lon1 : + Math.AngNormalize(lon1) + result.azi1 = Math.AngNormalize(azi1) + result.a12 = a12 + (outmask & Geodesics.DISTANCE) > 0 && (result.s12 = s12) + (outmask & Geodesics.LATITUDE) > 0 && (result.lat2 = lat2) + (outmask & Geodesics.LONGITUDE) > 0 && (result.lon2 = lon2) + (outmask & Geodesics.AZIMUTH) > 0 && (result.azi2 = azi2) + (outmask & Geodesics.REDUCEDLENGTH) > 0 && (result.m12 = m12) + if (outmask & Geodesics.GEODESICSCALE) > 0 + result.M12 = M12 + result.M21 = M21 + end + (outmask & Geodesics.AREA) > 0 && (result.S12 = S12) + result +end + +"""Private: General version of direct problem""" +function _GenDirect(self::Geodesics.Geodesic, lat1, lon1, azi1, arcmode, s12_a12, outmask) + outmask = Int(outmask) + # Automatically supply DISTANCE_IN if necessary + if ! arcmode + outmask |= Geodesics.DISTANCE_IN + end + line = GeodesicLines.GeodesicLine(self, lat1, lon1, azi1, caps=outmask) + GeodesicLines._GenPosition(line, arcmode, s12_a12, outmask) +end diff --git a/src/GeographicLib/inverse_line.jl b/src/GeographicLib/inverse_line.jl new file mode 100644 index 0000000..c92fecb --- /dev/null +++ b/src/GeographicLib/inverse_line.jl @@ -0,0 +1,84 @@ +"""Define a GeodesicLine object in terms of the direct geodesic +problem specified in terms of spherical arc length + +:param lat1: latitude of the first point in degrees +:param lon1: longitude of the first point in degrees +:param azi1: azimuth at the first point in degrees +:param s12: the distance from the first point to the second in +meters +:param caps: the :ref:`capabilities ` +:return: a :class:`~geographiclib.geodesicline.GeodesicLine` + +This function sets point 3 of the GeodesicLine to correspond to +point 2 of the direct geodesic problem. The default value of *caps* +is STANDARD | DISTANCE_IN, allowing direct geodesic problem to be +solved. +""" +function DirectLine(geod::Geodesic, lat1, lon1, azi1, s12, + caps = GeodesicCapability.STANDARD | GeodesicCapability.DISTANCE_IN) + _GenDirectLine(geod, lat1, lon1, azi1, false, s12, caps) +end + +"""Define a GeodesicLine object in terms of the direct geodesic +problem specified in terms of spherical arc length + +:param lat1: latitude of the first point in degrees +:param lon1: longitude of the first point in degrees +:param azi1: azimuth at the first point in degrees +:param a12: spherical arc length from the first point to the second +in degrees +:param caps: the :ref:`capabilities ` +:return: a :class:`~geographiclib.geodesicline.GeodesicLine` + +This function sets point 3 of the GeodesicLine to correspond to +point 2 of the direct geodesic problem. The default value of *caps* +is STANDARD | DISTANCE_IN, allowing direct geodesic problem to be +solved. +""" +function ArcDirectLine(geod::Geodesic, lat1, lon1, azi1, a12, + caps = GeodesicCapability.STANDARD | GeodesicCapability.DISTANCE_IN) + _GenDirectLine(geod, lat1, lon1, azi1, true, a12, caps) +end + +"""Define a GeodesicLine object in terms of the invese geodesic problem + +:param lat1: latitude of the first point in degrees +:param lon1: longitude of the first point in degrees +:param lat2: latitude of the second point in degrees +:param lon2: longitude of the second point in degrees +:param caps: the :ref:`capabilities ` +:return: a :class:`~geographiclib.geodesicline.GeodesicLine` + +This function sets point 3 of the GeodesicLine to correspond to +point 2 of the inverse geodesic problem. The default value of *caps* +is STANDARD | DISTANCE_IN, allowing direct geodesic problem to be +solved. +""" +function InverseLine(geod::Geodesic, lat1, lon1, lat2, lon2, + caps=GeodesicCapability.STANDARD | GeodesicCapability.DISTANCE_IN) + a12, _, salp1, calp1, _, _, _, _, _, _ = Geodesics._GenInverse(geod, lat1, lon1, lat2, lon2, 0) + azi1 = Math.atan2d(salp1, calp1) + if (caps & (Geodesics.OUT_MASK & Geodesics.DISTANCE_IN)) > 0 + caps |= Geodesics.DISTANCE + end + line = GeodesicLine(geod, lat1, lon1, azi1, caps=caps, salp1=salp1, calp1=calp1) + line = SetArc(line, a12) + line +end + +"""Private: general form of DirectLine""" +function _GenDirectLine(geod::Geodesic, lat1, lon1, azi1, arcmode::Bool, s12_a12, + caps=GeodesicCapability.STANDARD | GeodesicCapability.DISTANCE_IN) + # Automatically supply DISTANCE_IN if necessary + if !arcmode + caps |= Geodesics.DISTANCE_IN + end + line = GeodesicLine(geod, lat1, lon1, azi1, caps=caps) + line = if arcmode + SetArc(line, s12_a12) + else + SetDistance(line, s12_a12) + end + line +end + diff --git a/src/GeographicLib/results.jl b/src/GeographicLib/results.jl new file mode 100644 index 0000000..6ea3df5 --- /dev/null +++ b/src/GeographicLib/results.jl @@ -0,0 +1,56 @@ +# Struct holding results of the GeographicLib interface + +""" + Result + +Type holding the result of a forward or inverse calculation using a geodesic. + +If fields have not been calculated because of the presence/absence of any +flags when performing the calculation, the field will be `nothing`. + +### Fields + +These are user-accesible. + +- `lat1`: Latitude of starting point (°) +- `lon1`: Longitude of starting point (°) +- `lat2`: Latitude of end point (°) +- `lon2`: Longitude of end point (°) +- `a12`: Angular distace between start and end points (°) +- `s12`: Distance between start and end points (°) +- `azi1`: Forward azimuth from start point to end point (°) +- `azi2`: Forward azimuth at end point from start point along great cricle +- `m12`: Reduced length of the geodesic +- `M12`: First geodesic scale +- `M21`: Second geodesic scale +- `S12`: Area between the path from the first to the second point, and the equator (m²) +""" +mutable struct Result + "Latitude of starting point (°)" + lat1::Union{Float64,Nothing} + "Longitude of starting point (°)" + lon1::Union{Float64,Nothing} + "Latitude of end point (°)" + lat2::Union{Float64,Nothing} + "Longitude of end point (°)" + lon2::Union{Float64,Nothing} + "Angular distace between start and end points (°)" + a12::Union{Float64,Nothing} + "Distance between start and end points (°)" + s12::Union{Float64,Nothing} + "Forward azimuth from start point to end point (°)" + azi1::Union{Float64,Nothing} + "Forward azimuth at end point from start point along great cricle" + azi2::Union{Float64,Nothing} + "Reduced length of the geodesic" + m12::Union{Float64,Nothing} + "First geodesic scale" + M12::Union{Float64,Nothing} + "Second geodesic scale" + M21::Union{Float64,Nothing} + "Area between the path from the first to the second point, and the equator (m²)" + S12::Union{Float64,Nothing} +end + +Result() = Result(nothing, nothing, nothing, nothing, nothing, nothing, nothing, + nothing, nothing, nothing, nothing, nothing) diff --git a/src/GeographicLib/show.jl b/src/GeographicLib/show.jl new file mode 100644 index 0000000..a7b601e --- /dev/null +++ b/src/GeographicLib/show.jl @@ -0,0 +1,13 @@ +# Displaying custom types + +const SHOW_MAXLEN = maximum(x -> first(x) != '_' ? length(x) : 0, + String(f) for T in (Geodesic, GeodesicLine, Result, Polygon) for f in fieldnames(T)) + +function Base.show(io::IO, x::Union{Geodesic, GeodesicLine, Result, Polygon}) + println(io, typeof(x), ":") + for f in fieldnames(typeof(x)) + first(String(f)) == '_' && continue + v = getfield(x, f) + println(io, lpad(String(f), SHOW_MAXLEN), ": ", v === nothing ? "<>" : v) + end +end diff --git a/src/distances.jl b/src/distances.jl index 2f9a65e..1a7154c 100644 --- a/src/distances.jl +++ b/src/distances.jl @@ -1,3 +1,5 @@ +# Cartesian distances + using LinearAlgebra: norm """ @@ -32,5 +34,3 @@ distance(a::UTM, b::UTM, zone::Integer, isnorth::Bool, datum = wgs84) = distance distance(a, b::UTM, zone::Integer, isnorth::Bool, datum = wgs84) = distance(a, ECEFfromUTM(zone, isnorth, datum)(b), datum) distance(a::UTM, b, zone::Integer, isnorth::Bool, datum = wgs84) = distance(ECEFfromUTM(zone, isnorth, datum)(a), b, datum) - -# Also add geodesic distances here diff --git a/src/geodesics.jl b/src/geodesics.jl new file mode 100644 index 0000000..18645a6 --- /dev/null +++ b/src/geodesics.jl @@ -0,0 +1,173 @@ +# Great circle distances + +include("GeographicLib/GeographicLib.jl") +import .GeographicLib + +""" + GreatCircle + +A `GreatCircle` object stores a cache of values needed to quickly compute +great circles. +""" +struct GreatCircle + geod::GeographicLib.Geodesic + datum::Datum +end + +""" + GreatCircle(datum) -> gc + +Construct a `GreatCircle` `gc` from a datum. + +The object `gc` can then be used as a function to calculate great circle properties. +""" +GreatCircle(datum::Datum) = + (el = ellipsoid(datum); GreatCircle(GeographicLib.Geodesic(el.a, el.f), datum)) + +""" + (::GreatCircle)(a, b) -> (azi, baz, dist, angle) + +Compute the great circle between points `a` and `b`. The function returns a named +tuple giving the forward azimuth (`azi`, °) from `a` to `b`, +the backazimuth (`baz`, °) from `b` to `a`, the great circle distance (`dist`, m) +and the spherical equivalent angular distance (`angle`, °). + +#### Example + +Compute the distance `dist` (m) from Nelson's Column, Trafalgar Square, London to +the Empire State Building, New York, USA, as well as the forward azimuth (`azi`°) +from London, the backazimuth (`baz`°) from New York and the spherical equivalent angular +distance (`angle`°). + +``` +julia> using Geodesy + +julia> gc = GreatCircle(wgs84); + +julia> nelson = LatLon(51.5077, -0.1277); esb = LatLon(40.7484, -73.9857); + +julia> gc(nelson, esb) +(azi = -71.60766384249631, baz = 51.2691499649614, dist = 5.581417932416286e6, angle = 50.20783561306145) +``` + +Distances are typically accurate to the nanometre. + +--- + + (::GreatCircle)(a, b, zone, isnorth) + (::GreatCircle)(a, zone_a, isnorth_a, b, zone_b, isnorth_b) + +If one of `a` or `b` are `UTM`s, then provide the UTM `zone` integer and whether +the zone is in the northern hemisphere (`isnorth=true`) or not. If both `a` and +`b` are `UTM`s, then both will be taken to be in the same zone and hemisphere. + +If both `a` and `b` are `UTM`s but are not in the same zone and hemisphere, provide +the zone numbers `zone_a` and `zone_b`, respectively for points `a` and `b`, and +likewise `isnorth_a` and `isnorth_b` for the hemisphere. + +--- + + (::GreatCircle)(a, b, origin) + (::GreatCircle)(a, origin_a, b, origin_b) + +If one or both of `a` and `b` are `ENU`s, then provide the `origin` of the local +coordinate system. If both `a` and `b` are in `ENU` format but are related to +`LLA` with different origins, provide `origin_a` and `origin_b` respectively for +points `a` and `b`. + +--- + + (::GreatCircle)(utm, zone, isnorth, enu, origin) + (::GreatCircle)(enu, origin, utm, zone, isnorth) + +If one of the points is a `UTM` and the other an `ENU`, supply the `zone` +and `isnorth` of the former, and the `origin` of the latter. +""" +(g::GreatCircle)(a::Union{LatLon,LLA}, b::Union{LatLon,LLA}) = + GeographicLib.inverse(g.geod, a.lon, a.lat, b.lon, b.lat) + +# UTM conversions +(g::GreatCircle)(a::UTM, b, zone::Integer, isnorth::Bool) = + g(LLA(a, zone, isnorth, g.datum), b) +(g::GreatCircle)(a, b::UTM, zone::Integer, isnorth::Bool) = + g(a, LLA(b, zone, isnorth, g.datum)) +(g::GreatCircle)(a::UTM, b::UTM, zone::Integer, isnorth::Bool) = + g(a, zone, isnorth, b, zone, isnorth) +(g::GreatCircle)(a::UTM, zone_a::Integer, isnorth_a::Bool, + b::UTM, zone_b::Integer, isnorth_b::Bool) = + g(LLA(a, zone_a, isnorth_a, g.datum), LLA(b, zone_b, isnorth_b, g.datum)) + +# ENU conversions +(g::GreatCircle)(a::ENU, origin_a, b::ENU, origin_b) = + g(LLA(a, origin_a, g.datum), LLA(b, origin_b, g.datum)) +(g::GreatCircle)(a::ENU, b::ENU, origin) = g(a, origin, b, origin) +(g::GreatCircle)(a::ENU, b, origin) = g(LLA(a, origin, g.datum), b) +(g::GreatCircle)(a, b::ENU, origin) = g(a, LLA(b, origin, g.datum)) + +# ENU and UTM conversions combined +(g::GreatCircle)(a::ENU, origin, b::UTM, zone::Integer, isnorth::Bool) = + g(LLA(a, origin, g.datum), LLA(b, zone, isnorth, g.datum)) +(g::GreatCircle)(a::UTM, zone::Integer, isnorth::Bool, b::ENU, origin) = + g(LLA(a, zone, isnorth, g.datum), LLA(b, origin, g.datum)) + +# Other conversions +(g::GreatCircle)(a::Union{LatLon,LLA}, b) = g(a, LLA(b, g.datum)) +(g::GreatCircle)(a, b::Union{LatLon,LLA}) = g(LLA(a, g.datum), b) +(g::GreatCircle)(a, b) = g(LLA(a, g.datum), LLA(b, g.datum)) + +""" + (::GreatCircle)(p; azi, dist, angle) -> (lon, lat, baz, dist, angle) + +Compute the endpoint when travelling either `dist` m or an equivalent spherical angular +distance `angle`° along an azimuth `azi`° from starting point `p`. + +One and only one of `dist` and `angle` should be provided. `azi` must always be given. + +The function returns a named tuple giving the longitude (`lon`, °), latitude (`lat`, °), +backazimuth from the end point (`baz`, °), distance (`dist`, m) and spherical equivalent +angular distance (`angle`, °). + +#### Example + +Calculate the end point from travelling northeast from Nelson's Column, London for +1 km: + +``` +julia> using Geodesy + +julia> gc = GreatCircle(wgs84); + +julia> nelson = LatLon(51.5077, -0.1277); + +juila> gc(nelson, azi=45, dist=1000) +(lon = -0.11751395294238295, lat = 51.51405511443968, baz = -134.99202711278312, dist = 1000, angle = 0.008994870339334382) +``` + +--- + + (::GreatCircle)(p, zone, isnorth; azi, dist, angle) + +If `p` is a `UTM`, then provide the UTM `zone` number and whether the zone is in the +northern hemisphere (`isnorth=true`) or not. + +--- + + (::GreatCircle)(p, origin; azi, dist, angle) + +If `p` is a `ENU`, then provide the `origin` of the local ENU coordinate system. +""" +function (g::GreatCircle)(p::Union{LatLon,LLA}; azi, dist=nothing, angle=nothing) + sum(x -> x === nothing, (dist, angle)) == 1 || + throw(ArgumentError("must specify one and only one of dist or angle")) + if dist !== nothing + GeographicLib.forward(g.geod, p.lon, p.lat, azi, dist) + elseif angle !== nothing + GeographicLib.forward_deg(g.geod, p.lon, p.lat, azi, angle) + end +end + +# Conversions +(g::GreatCircle)(p::UTM, zone::Integer, isnorth::Bool; kwargs...) = + g(LLA(p, zone, isnorth, g.datum); kwargs...) +(g::GreatCircle)(p::ENU, origin; kwargs...) = g(LLA(p, origin, g.datum); kwargs...) +(g::GreatCircle)(p; kwargs...) = g(LLA(p, g.datum); kwargs...) diff --git a/test/geodesics.jl b/test/geodesics.jl new file mode 100644 index 0000000..9e73f3f --- /dev/null +++ b/test/geodesics.jl @@ -0,0 +1,89 @@ +@testset "Geodesics" begin + gc = GreatCircle(wgs84) + datum = wgs84 + ellps = ellipsoid(datum) + + # Construction + @test gc isa GreatCircle + @test gc.datum == datum + + # Invocation + @test_throws UndefKeywordError gc(LLA(0, 0, 0), dist=0) + @test_throws ArgumentError gc(LLA(0, 0, 0), azi=0) + @test_throws ArgumentError gc(LLA(0, 0, 0), azi=0, dist=0, angle=0) + + # Output + @test propertynames(gc(LLA(0, 0, 0), LLA(1, 1, 1))) == (:azi, :baz, :dist, :angle) + @test propertynames(gc(LLA(0, 0, 0), azi=0, dist=1)) == (:lon, :lat, :baz, :dist, :angle) + + # Conversion between LLA and other types + point1 = LLA(lat = 180rand() - 90, lon = 360rand(), alt=100_000rand()) + point2 = LLA(lat = 180rand() - 90, lon = 360rand(), alt=100_000rand()) + origin1 = point1 + zone1 = UTMZ(point1, datum).zone + isnorth1 = point1.lat >= 0 + origin2 = point2 + zone2 = UTMZ(point2, datum).zone + isnorth2 = point2.lat >= 0 + point_types1 = (identity, + x -> ECEF(x, datum), + LatLon, + x -> ENU(x, origin1, datum), + x -> UTMZ(x, datum), + x -> UTM(x, zone1, isnorth1, datum)) + point_types2 = (identity, + x -> ECEF(x, datum), + LatLon, + x -> ENU(x, origin2, datum), + x -> UTMZ(x, datum), + x -> UTM(x, zone2, isnorth2, datum)) + azi = 720rand() - 360 + dist = 100_000rand() + arc = 720rand() + + for type1 in point_types1 + tpoint1 = type1(point1) + # Approximate equality because of numerical error in conversion + # Forward + if tpoint1 isa ENU + @test gc(tpoint1, origin1, azi=azi, dist=dist) ≈ gc(point1, azi=azi, dist=dist) + @test gc(tpoint1, origin1, azi=azi, angle=arc) ≈ gc(point1, azi=azi, angle=arc) + elseif tpoint1 isa UTM + @test gc(tpoint1, zone1, isnorth1, azi=azi, angle=dist) ≈ gc(point1, azi=azi, angle=dist) + @test gc(tpoint1, zone1, isnorth1, azi=azi, angle=arc) ≈ gc(point1, azi=azi, angle=arc) + else + @test gc(tpoint1, azi=azi, dist=dist) ≈ gc(point1, azi=azi, dist=dist) + @test gc(tpoint1, azi=azi, angle=arc) ≈ gc(point1, azi=azi, angle=arc) + end + # Inverse + for type2 in point_types2 + result = gc(point1, point2) + tpoint2 = type2(point2) + if tpoint1 isa ENU + if tpoint2 isa ENU + @test gc(tpoint1, origin1, tpoint2, origin2) ≈ result + elseif tpoint2 isa UTM + @test gc(tpoint1, origin1, tpoint2, zone2, isnorth2) ≈ result + else + @test gc(tpoint1, tpoint2, origin1) ≈ result + end + elseif tpoint1 isa UTM + if tpoint2 isa ENU + @test gc(tpoint1, zone1, isnorth1, tpoint2, origin2) ≈ result + elseif tpoint2 isa UTM + @test gc(tpoint1, zone1, isnorth1, tpoint2, zone2, isnorth2) ≈ result + else + @test gc(tpoint1, tpoint2, zone1, isnorth1) ≈ result + end + else + if tpoint2 isa ENU + @test gc(tpoint1, tpoint2, origin2) ≈ result + elseif tpoint2 isa UTM + @test gc(tpoint1, tpoint2, zone2, isnorth2) ≈ result + else + @test gc(tpoint1, tpoint2) ≈ result + end + end + end + end +end diff --git a/test/runtests.jl b/test/runtests.jl index b42a14f..5b9ed45 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,7 +6,7 @@ using Test ################################################ # Interesting that this isn't in Base... -Base.isapprox(a::T, b::T; kwargs...) where {T<:Tuple} = all(ntuple(i->isapprox(a[i],b[i]), length(a)); kwargs...) +Base.isapprox(a::T, b::T; kwargs...) where {T<:Union{Tuple,NamedTuple}} = all(ntuple(i->isapprox(a[i],b[i]), length(a)); kwargs...) @testset "Geodesy" begin include("points.jl") @@ -17,4 +17,5 @@ Base.isapprox(a::T, b::T; kwargs...) where {T<:Tuple} = all(ntuple(i->isapprox(a include("transformations.jl") include("conversion.jl") include("datum_transformations.jl") + include("geodesics.jl") end # @testset "Geodesy"