Allow SpherePoint.bearingTo and offset at poles

For a SpherePoint constructed from a 3-vector, set longitude to 0
if at a pole and document the same. This allows bearingTo and offset
to be well defined for all cases.

Also change operator== to compare longitude even at the poles,
since different longitudes result in different behavior for bearingTo
and offset even at the poles.

include/lsst/afw/geom/SpherePoint.h

@@ -63,6 +63,8 @@ class SpherePoint final {
      * Construct a SpherePoint from a longitude and latitude.
      *
      * @param longitude The longitude of the point.
+     *                  It will be wrapped into the range [0, 2 π) radians,
+     *                  and `+/-inf` is wrapped to `nan`
      * @param latitude The latitude of the point. Must be in the
      *                 interval [-π/2, π/2] radians.
      *
@@ -77,6 +79,8 @@ class SpherePoint final {
      * Construct a SpherePoint from double longitude and latitude.
      *
      * @param longitude The longitude of the point.
+     *                  It will be wrapped into the range [0, 2 π) radians,
+     *                  and `+/-inf` is wrapped to `nan`
      * @param latitude The latitude of the point. Must be in the
      *                 interval [-π/2, π/2] radians.
      * @param units The units of longitude and latitude
@@ -85,6 +89,9 @@ class SpherePoint final {
      *         Thrown if `latitude` isout of range.
      *
      * @exceptsafe Provides strong exception guarantee.
+     *
+     * @warning If longitude is +/- inf then the resulting longitude is `NaN`
+     * since infinity cannot be wrapped into the range
      */
     SpherePoint(double longitude, double latitude, AngleUnit units);
 
@@ -98,6 +105,9 @@ class SpherePoint final {
      * @throws pex::exceptions::InvalidParameterError
      *         Thrown if `vector` is the zero vector.
      *
+     * @note If the SpherePoint is at a pole then longitude is set to 0.
+     * That provides predictable behavior for @ref bearingTo and @ref offset.
+     *
      * @exceptsafe Provides strong exception guarantee.
      */
     explicit SpherePoint(Point3D const& vector);
@@ -148,10 +158,8 @@ class SpherePoint final {
     /**
      * The longitude of this point.
      *
-     * If this point is at a coordinate pole, the longitude is undefined, and
-     * this method may return any value. If the SpherePoint implementation
-     * allows multiple values of longitude from a pole, they shall all be
-     * treated as valid representations of the same point.
+     * If this point is at a pole, the longitude is the value provided to the constructor,
+     * if using a longitude, latitude constructor, or 0 if using the Point3D constructor.
      *
      * @returns the longitude, in the interval [0, 2π) radians.
      *
@@ -202,10 +210,8 @@ class SpherePoint final {
      * @exceptsafe Shall not throw exceptions.
      */
     bool atPole() const noexcept {
-        // Unit tests indicate we don't need to worry about epsilon-errors
-        // Objects constructed from lat=90*degrees, <0,0,1>, etc. all have
-        // atPole() = true. More complex operations like bearingTo have also
-        // been tested near the poles with no ill effects
+        // Unit tests indicate we don't need to worry about epsilon-errors, in that
+        // Objects constructed from lat=90*degrees, <0,0,1>, etc. all have atPole() = true.
         return fabs(_latitude) >= HALFPI;
     }
 
@@ -225,12 +231,12 @@ class SpherePoint final {
      */
 
     /**
-     * `true` if two points represent the same position.
+     * `true` if two points have the same longitude and latitude.
      *
-     * Points are considered equal if and only if they represent the same
-     * location, regardless of how they were constructed. In particular,
-     * representations of the north or south pole with different longitudes
-     * are considered equal.
+     * @note Two points at the same pole will compare unequal if they have different longitudes,
+     * despite representing the same point on the unit sphere.
+     * This is important because the behavior of @ref offset and @ref bearingTo
+     * depend on longitude, even at the pole.
      *
      * @param other the point to test for equality
      * @returns true if this point matches `other` exactly, false otherwise
@@ -254,20 +260,19 @@ class SpherePoint final {
     bool operator!=(SpherePoint const& other) const noexcept;
 
     /**
-     * Direction from one point to another.
+     * Orientation at this point of the great circle arc from this point to another point.
+     *
+     * The orientation is measured in a plane tangent to the sphere at this point, with 0 degrees along
+     * the direction of increasing longitude and 90 degrees along the direction of increasing latitude.
+     * Thus for any point except the pole, the arc is due east at 0 degrees and due north at 90 degrees.
+     * To understand the behavior at the poles it may help to imagine the point has been shifted slightly
+     * by decreasing the absolute value of its latitude.
      *
-     * This method finds the shortest (great-circle) arc between two
-     * points, and characterizes its direction by the angle between
-     * it and a line of latitude passing through this point. 0 represents
-     * due east, &pi;/2 represents due north. If the points are on
-     * opposite sides of the sphere, the bearing may be any value.
+     * If the points are on opposite sides of the sphere then the returned bearing may be any value.
      *
      * @param other the point to which to measure the bearing
      * @returns the direction, as defined above, in the interval [0, 2&pi;).
      *
-     * @throws pex::exceptions::DomainError
-     *         Thrown if `this.atPole()`.
-     *
      * @exceptsafe Provides strong exception guarantee.
      *
      * @note For two points `A` and `B`, `A.bearingTo(B)` will in
@@ -296,7 +301,8 @@ class SpherePoint final {
      * @param axis a point defining the north pole of the rotation axis.
      * @param amount the amount of rotation, where positive values
      *               represent right-handed rotations around `axis`.
-     * @returns a new point created by rotating this point
+     * @returns a new point created by rotating this point. If the point is at the pole
+     * then its longitude will be 0.
      *
      * @exceptsafe Shall not throw exceptions.
      */
@@ -305,20 +311,17 @@ class SpherePoint final {
     /**
      * Return a point offset from this one along a great circle.
      *
-     * For any point `A` not at a coordinate pole, and any two angles `b`
-     * and `delta`, `A.bearingTo(A.offset(b, delta))` = `b` and
-     * `A.separationTo(A.offset(b, delta))` = `delta`.
+     * For any point `A` not at a coordinate pole, any angle `bearing` and any non-negative angle `amount`,
+     * if `B` = A.offset(bearing, amount)` then `A.bearingTo(B) = amount` and `A.separationTo(B) = amount`.
+     * Negative values of `amount` are supported in the obvious manner:
+     * `A.offset(bearing, delta) = A.offset(bearing + 180*degrees, -delta)`
      *
      * @param bearing the direction in which to move this point, following
      *                the conventions described in @ref bearingTo.
      * @param amount the distance by which to move along the great
      *               circle defined by `bearing`
-     * @returns a new point created by shifting this point
-     *
-     * @throws pex::exceptions::DomainError
-     *         Thrown if `this.atPole()`.
-     * @throws pex::exceptions::InvalidParameterError
-     *         Thrown if `amount` is negative.
+     * @returns a new point created by rotating this point. If the point is at the pole
+     * then its longitude will be 0.
      *
      * @exceptsafe Provides strong exception guarantee.
      */

src/geom/SpherePoint.cc

@@ -90,10 +90,14 @@ SpherePoint::SpherePoint(Point3D const& vector) {
     double const y = vector.getY() / norm;
     double const z = vector.getZ() / norm;
 
-    // Need to convert to Angle, Angle::wrap, and convert back to radians
-    //     to handle _longitude = -1e-16 without code duplication
-    _longitude = (atan2(y, x) * radians).wrap().asRadians();
     _latitude = asin(z);
+    if (!atPole()) {
+        // Need to convert to Angle, Angle::wrap, and convert back to radians
+        //     to handle _longitude = -1e-16 without code duplication
+        _longitude = (atan2(y, x) * radians).wrap().asRadians();
+    } else {
+        _longitude = 0;
+    }
 }
 
 SpherePoint::SpherePoint(SpherePoint const& other) noexcept = default;
@@ -124,22 +128,12 @@ bool SpherePoint::isFinite() const noexcept { return isfinite(_longitude) && isf
 bool SpherePoint::operator==(SpherePoint const& other) const noexcept {
     // Deliberate override of Style Guide 5-12
     // Approximate FP comparison would make object equality intransitive
-    if (atPole()) {
-        return _latitude == other._latitude;
-    } else {
-        return _longitude == other._longitude && _latitude == other._latitude;
-    }
+    return _longitude == other._longitude && _latitude == other._latitude;
 }
 
 bool SpherePoint::operator!=(SpherePoint const& other) const noexcept { return !(*this == other); }
 
 Angle SpherePoint::bearingTo(SpherePoint const& other) const {
-    if (atPole()) {
-        stringstream buffer;
-        buffer << "Cannot calculate offset from pole " << *this << ".";
-        throw pexExcept::DomainError(buffer.str());
-    }
-
     Angle const deltaLon = other.getLongitude() - this->getLongitude();
 
     double const sinDelta1 = sin(getLatitude().asRadians());
@@ -167,16 +161,6 @@ SpherePoint SpherePoint::rotated(SpherePoint const& axis, Angle const& amount) c
 
 SpherePoint SpherePoint::offset(Angle const& bearing, Angle const& amount) const {
     double const phi = bearing.asRadians();
-    if (atPole()) {
-        stringstream buffer;
-        buffer << "Cannot define offset direction from pole " << *this << ".";
-        throw pexExcept::DomainError(buffer.str());
-    }
-    if (amount < 0.0) {
-        stringstream buffer;
-        buffer << "Negative offset of " << amount.asDegrees() << " degrees is not allowed.";
-        throw pexExcept::InvalidParameterError(buffer.str());
-    }
 
     // let v = vector in the direction bearing points (tangent to surface of sphere)
     // To do the rotation, use rotate() method.