Angle.h

Go to the documentation of this file.

/*
 * Developed for the LSST Data Management System.
 * This product includes software developed by the LSST Project
 * (https://www.lsst.org).
 * See the COPYRIGHT file at the top-level directory of this distribution
 * for details of code ownership.
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <https://www.gnu.org/licenses/>.
 */
#ifndef LSST_GEOM_ANGLE_H
#define LSST_GEOM_ANGLE_H
 
#include <cmath>
#include <iostream>
#include <type_traits>
 
#include "lsst/sphgeom/Angle.h"
#include "boost/math/constants/constants.hpp"
 
namespace lsst {
namespace geom {
 
/*
 * None of C99, C++98, and C++11 define M_PI, so we'll do it ourselves
 */
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wunused-variable"
double constexpr PI = boost::math::constants::pi<double>();
double constexpr TWOPI = boost::math::constants::pi<double>() * 2.0;
double constexpr HALFPI = boost::math::constants::pi<double>() * 0.5;
double constexpr ONE_OVER_PI = 1.0 / boost::math::constants::pi<double>();
// sqrt is not a constexpr on OS X
double const SQRTPI = sqrt(boost::math::constants::pi<double>());
double const INVSQRTPI = 1.0 / sqrt(boost::math::constants::pi<double>());
double constexpr ROOT2 = boost::math::constants::root_two<double>();  // sqrt(2)
#pragma clang diagnostic pop
 
// These shouldn't be necessary if the Angle class is used, but sometimes you just need
// them.  Better to define them once here than have *180/PI throughout the code...
inline constexpr double degToRad(double x) noexcept { return x * PI / 180.; }
inline constexpr double radToDeg(double x) noexcept { return x * 180. / PI; }
inline constexpr double radToArcsec(double x) noexcept { return x * 3600. * 180. / PI; }
inline constexpr double radToMas(double x) noexcept { return x * 1000. * 3600. * 180. / PI; }
inline constexpr double arcsecToRad(double x) noexcept { return (x / 3600.) * PI / 180.; }
inline constexpr double masToRad(double x) noexcept { return (x / (1000. * 3600.)) * PI / 180.; }
 
class Angle;
class AngleUnit final {
    friend class Angle;
    template <typename T>
    friend constexpr Angle operator*(T lhs, AngleUnit rhs) noexcept;
 
public:
    // Current implementation does not throw exceptions, but somebody may want
    // to add input validation later
    explicit constexpr AngleUnit(double val) : _val(val) {}
 
    constexpr bool operator==(AngleUnit const& rhs) const noexcept;
 
    std::size_t hash_value() const noexcept { return std::hash<double>()(_val); }
 
private:
    double _val;
};
 
inline constexpr bool AngleUnit::operator==(AngleUnit const& rhs) const noexcept {
    return (_val == rhs._val);
}
 
AngleUnit constexpr radians = AngleUnit(1.0);                     
AngleUnit constexpr degrees = AngleUnit(PI / 180.0);              
AngleUnit constexpr hours = AngleUnit(PI * 15.0 / 180.0);         
AngleUnit constexpr arcminutes = AngleUnit(PI / 60 / 180.0);      
AngleUnit constexpr arcseconds = AngleUnit(PI / 180.0 / 3600.0);  
// Note: if we use PI / 180.0 / 3.6e6 then 60*60*180*1000*lsst.geom.milliarcseconds
// does not test exactly equal to 180*lsst.geom.degrees
AngleUnit constexpr milliarcseconds =
        AngleUnit(PI / (180.0 * 3.6e6));  
 
class Angle final {
    friend class AngleUnit;
 
public:
    explicit constexpr Angle(double val, AngleUnit units = radians) noexcept : _val(val * units._val) {}
 
    constexpr Angle() noexcept : _val(0) {}
 
    explicit Angle(sphgeom::Angle const & sgAngle) noexcept : Angle(sgAngle.asRadians()) {}
 
    constexpr Angle(Angle const& other) noexcept = default;
 
    constexpr Angle(Angle&& other) noexcept = default;
 
    Angle& operator=(Angle const& other) noexcept = default;
 
    Angle& operator=(Angle&& other) noexcept = default;
 
    ~Angle() = default;
 
    constexpr operator double() const noexcept { return _val; }
 
    operator sphgeom::Angle() const noexcept { return sphgeom::Angle::fromRadians(_val); }
 
    constexpr double asAngularUnits(AngleUnit const& units) const noexcept { return _val / units._val; }
 
    constexpr double asRadians() const noexcept { return asAngularUnits(radians); }
 
    constexpr double asDegrees() const noexcept { return asAngularUnits(degrees); }
 
    constexpr double asHours() const noexcept { return asAngularUnits(hours); }
 
    constexpr double asArcminutes() const noexcept { return asAngularUnits(arcminutes); }
 
    constexpr double asArcseconds() const noexcept { return asAngularUnits(arcseconds); }
 
    constexpr double asMilliarcseconds() const noexcept { return asAngularUnits(milliarcseconds); }
 
    Angle wrap() const noexcept;
 
    Angle wrapCtr() const noexcept;
 
    Angle wrapNear(Angle const& refAng) const noexcept;
 
    Angle separation(Angle const& other) const noexcept;
 
#define ANGLE_OPUP_TYPE(OP, TYPE)                \
    Angle& operator OP(TYPE const& d) noexcept { \
        _val OP d;                               \
        return *this;                            \
    }
 
 
    ANGLE_OPUP_TYPE(*=, double)
    ANGLE_OPUP_TYPE(*=, int)
 
 
    ANGLE_OPUP_TYPE(+=, double)
    ANGLE_OPUP_TYPE(+=, int)
 
 
    ANGLE_OPUP_TYPE(-=, double)
    ANGLE_OPUP_TYPE(-=, int)
 
#undef ANGLE_OPUP_TYPE
 
#define ANGLE_COMP(OP) \
    constexpr bool operator OP(const Angle& rhs) const noexcept { return _val OP rhs._val; }
 
 
    ANGLE_COMP(==)
    ANGLE_COMP(!=)
 
 
    ANGLE_COMP(<=)
    ANGLE_COMP(>=)
    ANGLE_COMP(<)
    ANGLE_COMP(>)
 
#undef ANGLE_COMP
 
    std::size_t hash_value() const noexcept { return std::hash<double>()(_val); }
 
private:
    double _val;
};
 
/*
 * Operators for Angles.
 */
#define ANGLE_OP(OP)                                                    \
    inline constexpr Angle operator OP(Angle a, Angle d) noexcept {     \
        return Angle(static_cast<double>(a) OP static_cast<double>(d)); \
    }
 
// We need both int and double versions to avoid ambiguous overloading due to
// implicit conversion of Angle to double
#define ANGLE_OP_TYPE(OP, TYPE)                                    \
    inline constexpr Angle operator OP(Angle a, TYPE d) noexcept { \
        return Angle(static_cast<double>(a) OP d);                 \
    }                                                              \
                                                                   \
    inline constexpr Angle operator OP(TYPE d, Angle a) noexcept { \
        return Angle(d OP static_cast<double>(a));                 \
    }
 
ANGLE_OP(+)
ANGLE_OP(-)
ANGLE_OP(*)
 
ANGLE_OP_TYPE(*, double)
ANGLE_OP_TYPE(*, int)
 
#undef ANGLE_OP
#undef ANGLE_OP_TYPE
 
inline constexpr Angle operator-(Angle angle) { return Angle(-static_cast<double>(angle)); }
 
// Apparently @relatesalso doesn't work with grouping
inline constexpr Angle operator/(Angle a, int d) noexcept { return Angle(static_cast<double>(a) / d); }
 
inline constexpr Angle operator/(Angle a, double d) noexcept { return Angle(static_cast<double>(a) / d); }
 
// Division is different.  Don't allow division by an Angle
template <typename T>
constexpr double operator/(T const lhs, Angle rhs) noexcept = delete;
 
std::ostream& operator<<(std::ostream& s, Angle a);
 
template <typename T>
inline constexpr bool isAngle(T) noexcept {
    return std::is_base_of<Angle, T>::value;
};
 
template <typename T>
inline constexpr Angle operator*(T lhs, AngleUnit rhs) noexcept {
    static_assert(std::is_arithmetic<T>::value,
                  "Only numeric types may be multiplied by an AngleUnit to create an Angle!");
    return Angle(lhs * rhs._val);
}
 
// Inline method definitions, placed last in order to benefit from Angle's full API
 
inline Angle Angle::wrap() const noexcept {
    double wrapped = std::fmod(_val, TWOPI);
    // wrapped is in the range (-TWOPI, TWOPI)
    if (wrapped < 0.0) wrapped += TWOPI;
    // if wrapped is small enough, adding 2 pi gives 2 pi
    if (wrapped >= TWOPI) wrapped = 0.0;
    return wrapped * radians;
}
 
inline Angle Angle::wrapCtr() const noexcept {
    double wrapped = std::fmod(_val, TWOPI);
    // wrapped is in the range [-TWOPI, TWOPI]
    if (wrapped < -PI) {
        wrapped += TWOPI;
        if (wrapped >= PI) {
            // handle roundoff error, however unlikely
            wrapped = -PI;
        }
    } else if (wrapped >= PI) {
        wrapped -= TWOPI;
        if (wrapped < -PI) {
            // handle roundoff error, however unlikely
            wrapped = -PI;
        }
    }
    return wrapped * radians;
}
 
inline Angle Angle::wrapNear(Angle const& refAng) const noexcept {
    // compute (this - refAng).wrapCtr() + refAng
    // which is correct except for roundoff error at the edges
    double const refAngRad = refAng.asRadians();
    double wrapped = (*this - refAng).wrapCtr().asRadians() + refAngRad;
 
    // roundoff can cause slightly out-of-range values; fix those as bast we can;
    // note that both conditionals are wanted, since either or both may execute
    // (though if/else could be used if the lower limit was squishy for radians)
    if (wrapped - refAngRad >= PI) {
        wrapped -= TWOPI;
    }
    if (wrapped - refAngRad < -PI) {
        wrapped += TWOPI;
    }
    return wrapped * radians;
}
 
inline Angle Angle::separation(Angle const& other) const noexcept { return (*this - other).wrapCtr(); }
 
}  // namespace geom
}  // namespace lsst
 
namespace std {
template <>
struct hash<lsst::geom::AngleUnit> {
    using argument_type = lsst::geom::AngleUnit;
    using result_type = size_t;
    size_t operator()(argument_type const& x) const noexcept { return x.hash_value(); }
};
 
template <>
struct hash<lsst::geom::Angle> {
    using argument_type = lsst::geom::Angle;
    using result_type = size_t;
    size_t operator()(argument_type const& x) const noexcept { return x.hash_value(); }
};
}  // namespace std
 
#endif  // if !defined(LSST_GEOM_ANGLE_H)