Celestial.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <cmath>
#include <iostream>
#include <string>
 
#include "NumCpp/Coordinates/Cartesian.hpp"
#include "NumCpp/Core/Internal/Error.hpp"
#include "NumCpp/Core/Types.hpp"
#include "NumCpp/Functions/deg2rad.hpp"
#include "NumCpp/Functions/dot.hpp"
#include "NumCpp/Functions/rad2deg.hpp"
#include "NumCpp/NdArray.hpp"
#include "NumCpp/Utils/essentiallyEqual.hpp"
#include "NumCpp/Utils/num2str.hpp"
#include "NumCpp/Utils/sqr.hpp"
#include "NumCpp/Vector/Vec3.hpp"
 
namespace nc::coordinates::reference_frames
{
    //================================================================================
    class RA
    {
    public:
        //============================================================================
        RA() = default;
 
        //============================================================================
        explicit RA(double inDegrees) :
            degrees_(inDegrees),
            radians_(deg2rad(inDegrees))
        {
            if (inDegrees < 0 || inDegrees >= 360)
            {
                THROW_INVALID_ARGUMENT_ERROR("input degrees must be of the range [0, 360)");
            }
 
            hours_                  = static_cast<uint8>(std::floor(degrees_ / 15.));
            const double decMinutes = (degrees_ - static_cast<double>(hours_) * 15.) * 4.;
            minutes_                = static_cast<uint8>(std::floor(decMinutes));
            seconds_                = static_cast<double>((decMinutes - static_cast<double>(minutes_)) * 60.);
        }
 
        //============================================================================
        RA(uint8 inHours, uint8 inMinutes, double inSeconds) noexcept :
            hours_(inHours),
            minutes_(inMinutes),
            seconds_(inSeconds)
        {
            degrees_ = static_cast<double>(hours_) * 15. + static_cast<double>(minutes_) / 4. + seconds_ / 240.;
            radians_ = deg2rad(degrees_);
        }
 
        //============================================================================
        [[nodiscard]] double radians() const noexcept
        {
            return radians_;
        }
 
        //============================================================================
        [[nodiscard]] double degrees() const noexcept
        {
            return degrees_;
        }
 
        //============================================================================
        [[nodiscard]] uint8 hours() const noexcept
        {
            return hours_;
        }
 
        //============================================================================
        [[nodiscard]] uint8 minutes() const noexcept
        {
            return minutes_;
        }
 
        //============================================================================
        [[nodiscard]] double seconds() const noexcept
        {
            return seconds_;
        }
 
        //============================================================================
        [[nodiscard]] std::string str() const
        {
            std::string out =
                "RA hms: " + utils::num2str(hours_) + " hours, " + utils::num2str(minutes_) + " minutes, ";
            out += utils::num2str(seconds_) + " seconds\nRA degrees: " + utils::num2str(degrees_) + '\n';
            return out;
        }
 
        //============================================================================
        void print() const
        {
            std::cout << *this;
        }
 
        //============================================================================
        bool operator==(const RA& inRhs) const noexcept
        {
            return utils::essentiallyEqual(degrees_, inRhs.degrees_);
        }
 
        //============================================================================
        bool operator!=(const RA& inRhs) const noexcept
        {
            return !(*this == inRhs);
        }
 
        //============================================================================
        friend std::ostream& operator<<(std::ostream& inStream, const RA& inRa)
        {
            inStream << inRa.str();
            return inStream;
        }
 
    private:
        //====================================Attributes==============================
        uint8  hours_{ 0 };
        uint8  minutes_{ 0 };
        double seconds_{ 0. };
        double degrees_{ 0. };
        double radians_{ 0. };
    };
 
    //================================================================================
    class Dec
    {
    public:
        //================================================================================
        enum class Sign
        {
            NEGATIVE = 0,
            POSITIVE
        };
 
        //============================================================================
        Dec() = default;
 
        //============================================================================
        explicit Dec(double inDegrees) :
            degrees_(inDegrees),
            radians_(deg2rad(inDegrees))
        {
            if (inDegrees < -90 || inDegrees > 90)
            {
                THROW_INVALID_ARGUMENT_ERROR("input degrees must be of the range [-90, 90]");
            }
 
            sign_                   = degrees_ < 0 ? Sign::NEGATIVE : Sign::POSITIVE;
            const double absDegrees = std::abs(degrees_);
            degreesWhole_           = static_cast<uint8>(std::floor(absDegrees));
 
            const double decMinutes = (absDegrees - static_cast<double>(degreesWhole_)) * 60.;
            minutes_                = static_cast<uint8>(std::floor(decMinutes));
            seconds_                = (decMinutes - static_cast<double>(minutes_)) * 60.;
        }
 
        //============================================================================
        Dec(Sign inSign, uint8 inDegrees, uint8 inMinutes, double inSeconds) noexcept :
            sign_(inSign),
            degreesWhole_(inDegrees),
            minutes_(inMinutes),
            seconds_(inSeconds)
        {
            degrees_ = static_cast<double>(degreesWhole_) + static_cast<double>(minutes_) / 60. + seconds_ / 3600.;
            degrees_ *= sign_ == Sign::NEGATIVE ? -1 : 1;
 
            radians_ = deg2rad(degrees_);
        }
 
        //============================================================================
        [[nodiscard]] Sign sign() const noexcept
        {
            return sign_;
        }
 
        //============================================================================
        [[nodiscard]] double degrees() const noexcept
        {
            return degrees_;
        }
 
        //============================================================================
        [[nodiscard]] double radians() const noexcept
        {
            return radians_;
        }
 
        //============================================================================
        [[nodiscard]] uint8 degreesWhole() const noexcept
        {
            return degreesWhole_;
        }
 
        //============================================================================
        [[nodiscard]] uint8 minutes() const noexcept
        {
            return minutes_;
        }
 
        //============================================================================
        [[nodiscard]] double seconds() const noexcept
        {
            return seconds_;
        }
 
        //============================================================================
        [[nodiscard]] std::string str() const
        {
            std::string strSign = sign_ == Sign::NEGATIVE ? "-" : "+";
            std::string out     = "Dec dms: " + strSign + utils::num2str(degreesWhole_) + " degrees, " +
                              utils::num2str(minutes_) + " minutes, ";
            out += utils::num2str(seconds_) + " seconds\nDec degrees = " + utils::num2str(degrees_) + '\n';
            return out;
        }
 
        //============================================================================
        void print() const
        {
            std::cout << *this;
        }
 
        //============================================================================
        bool operator==(const Dec& inRhs) const noexcept
        {
            return utils::essentiallyEqual(degrees_, inRhs.degrees_);
        }
 
        //============================================================================
        bool operator!=(const Dec& inRhs) const noexcept
        {
            return !(*this == inRhs);
        }
 
        //============================================================================
        friend std::ostream& operator<<(std::ostream& inStream, const Dec& inDec)
        {
            inStream << inDec.str();
            return inStream;
        }
 
    private:
        //====================================Attributes==============================
        Sign   sign_{ Sign::POSITIVE };
        uint8  degreesWhole_{ 0 };
        uint8  minutes_{ 0 };
        double seconds_{ 0. };
        double degrees_{ 0. };
        double radians_{ 0. };
    };
 
    //================================================================================
    class Celestial
    {
    public:
        //============================================================================
        Celestial() = default;
 
        //============================================================================
        Celestial(double inRaDegrees, double inDecDegrees) :
            ra_(inRaDegrees),
            dec_(inDecDegrees)
        {
            polarToCartesian();
        }
 
        //============================================================================
        Celestial(uint8     inRaHours,
                  uint8     inRaMinutes,
                  double    inRaSeconds,
                  Dec::Sign inSign,
                  uint8     inDecDegreesWhole,
                  uint8     inDecMinutes,
                  double    inDecSeconds) :
            ra_(inRaHours, inRaMinutes, inRaSeconds),
            dec_(inSign, inDecDegreesWhole, inDecMinutes, inDecSeconds)
        {
            polarToCartesian();
        }
 
        //============================================================================
        Celestial(const RA& inRA, const Dec& inDec) noexcept :
            ra_(inRA),
            dec_(inDec)
        {
            polarToCartesian();
        }
 
        //============================================================================
        Celestial(double inX, double inY, double inZ) :
            x_(inX),
            y_(inY),
            z_(inZ)
        {
            cartesianToPolar();
        }
 
        //============================================================================
        Celestial(const Cartesian& inCartesianVector) :
            x_(inCartesianVector.x),
            y_(inCartesianVector.y),
            z_(inCartesianVector.z)
        {
            cartesianToPolar();
        }
 
        //============================================================================
        Celestial(const Vec3& inCartesianVector) :
            x_(inCartesianVector.x),
            y_(inCartesianVector.y),
            z_(inCartesianVector.z)
        {
            cartesianToPolar();
        }
 
        //============================================================================
        Celestial(const NdArray<double>& inCartesianVector)
        {
            if (inCartesianVector.size() != 3)
            {
                THROW_INVALID_ARGUMENT_ERROR("NdArray input must be of length 3.");
            }
 
            x_ = inCartesianVector[0];
            y_ = inCartesianVector[1];
            z_ = inCartesianVector[2];
 
            cartesianToPolar();
        }
 
        //============================================================================
        [[nodiscard]] const Dec& dec() const noexcept
        {
            return dec_;
        }
 
        //============================================================================
        [[nodiscard]] const RA& ra() const noexcept
        {
            return ra_;
        }
 
        //============================================================================
        [[nodiscard]] double x() const noexcept
        {
            return x_;
        }
 
        //============================================================================
        [[nodiscard]] double y() const noexcept
        {
            return y_;
        }
 
        //============================================================================
        [[nodiscard]] double z() const noexcept
        {
            return z_;
        }
 
        //============================================================================
        [[nodiscard]] NdArray<double> xyz() const
        {
            NdArray<double> out = { x_, y_, z_ };
            return out;
        }
 
        //============================================================================
        [[nodiscard]] double degreeSeperation(const Celestial& inOtherCelestial) const
        {
            return rad2deg(radianSeperation(inOtherCelestial));
        }
 
        //============================================================================
        [[nodiscard]] double degreeSeperation(const NdArray<double>& inVector) const
        {
            return rad2deg(radianSeperation(inVector));
        }
 
        //============================================================================
        [[nodiscard]] double radianSeperation(const Celestial& inOtherCelestial) const
        {
            return std::acos(dot(xyz(), inOtherCelestial.xyz()).item());
        }
 
        //============================================================================
        [[nodiscard]] double radianSeperation(const NdArray<double>& inVector) const
        {
            if (inVector.size() != 3)
            {
                THROW_INVALID_ARGUMENT_ERROR("input vector must be of length 3.");
            }
 
            return std::acos(dot(xyz(), inVector.flatten()).item());
        }
 
        //============================================================================
        [[nodiscard]] std::string str() const
        {
            std::string returnStr;
            returnStr = ra_.str();
            returnStr += dec_.str();
            returnStr += "Cartesian = " + xyz().str();
            return returnStr;
        }
 
        //============================================================================
        void print() const
        {
            std::cout << *this;
        }
 
        //============================================================================
        bool operator==(const Celestial& inRhs) const noexcept
        {
            return ra_ == inRhs.ra_ && dec_ == inRhs.dec_;
        }
 
        //============================================================================
        bool operator!=(const Celestial& inRhs) const noexcept
        {
            return !(*this == inRhs);
        }
 
        //============================================================================
        friend std::ostream& operator<<(std::ostream& inStream, const Celestial& inCoord)
        {
            inStream << inCoord.str();
            return inStream;
        }
 
    private:
        //====================================Attributes==============================
        RA     ra_{};
        Dec    dec_{};
        double x_{ 1. };
        double y_{ 0. };
        double z_{ 0. };
 
        //============================================================================
        void cartesianToPolar()
        {
            double degreesRa = rad2deg(std::atan2(y_, x_));
            if (degreesRa < 0)
            {
                degreesRa += 360;
            }
            ra_ = RA(degreesRa);
 
            const double r          = std::sqrt(utils::sqr(x_) + utils::sqr(y_) + utils::sqr(z_));
            const double degreesDec = rad2deg(std::asin(z_ / r));
            dec_                    = Dec(degreesDec);
        }
 
        //============================================================================
        void polarToCartesian() noexcept
        {
            const double raRadians  = deg2rad(ra_.degrees());
            const double decRadians = deg2rad(dec_.degrees());
 
            x_ = std::cos(raRadians) * std::cos(decRadians);
            y_ = std::sin(raRadians) * std::cos(decRadians);
            z_ = std::sin(decRadians);
        }
    };
} // namespace nc::coordinates::reference_frames