spherical.hpp

Go to the documentation of this file.

#pragma once

 
#include "../astro_context.hpp"
#include "../centers.hpp"
#include "../constants.hpp"
#include "../frames.hpp"
#include "../time.hpp"
#include "geodetic.hpp"
 
#include <qtty/qtty.hpp>
 
#include "detail/stream.hpp"
 
#include <cmath>
#include <ostream>
#include <type_traits>
 
// Forward-declare cartesian Position to avoid circular include.
namespace siderust {
namespace cartesian {
template <typename C, typename F, typename U> struct Position;
}
} // namespace siderust
 
namespace siderust {
namespace spherical {

template <typename F> struct Direction {
  static_assert(frames::is_frame_v<F>, "F must be a valid frame tag");
 
private:
  qtty::Degree azimuth_; 
  qtty::Degree polar_;   
 
public:
  Direction() : azimuth_(qtty::Degree(0)), polar_(qtty::Degree(0)) {}
 
  Direction(qtty::Degree azimuth, qtty::Degree polar) : azimuth_(azimuth), polar_(polar) {}

  static constexpr siderust_frame_t frame_id() { return frames::FrameTraits<F>::ffi_id; }
  static constexpr const char *frame_name() { return frames::FrameTraits<F>::name(); }

  template <typename F_ = F, std::enable_if_t<frames::has_ra_dec_v<F_>, int> = 0>
  qtty::Degree ra() const {
    return azimuth_;
  }
 
  template <typename F_ = F, std::enable_if_t<frames::has_ra_dec_v<F_>, int> = 0>
  qtty::Degree dec() const {
    return polar_;
  }


  template <typename F_ = F, std::enable_if_t<frames::has_az_alt_v<F_>, int> = 0>
  qtty::Degree az() const {
    return azimuth_;
  }
 
  template <typename F_ = F, std::enable_if_t<frames::has_az_alt_v<F_>, int> = 0>
  qtty::Degree al() const {
    return polar_;
  }
 
  template <typename F_ = F, std::enable_if_t<frames::has_az_alt_v<F_>, int> = 0>
  qtty::Degree alt() const {
    return polar_;
  }
 
  template <typename F_ = F, std::enable_if_t<frames::has_az_alt_v<F_>, int> = 0>
  qtty::Degree altitude() const {
    return polar_;
  }


  template <typename F_ = F, std::enable_if_t<frames::has_lon_lat_v<F_>, int> = 0>
  qtty::Degree lon() const {
    return azimuth_;
  }
 
  template <typename F_ = F, std::enable_if_t<frames::has_lon_lat_v<F_>, int> = 0>
  qtty::Degree lat() const {
    return polar_;
  }
 
  template <typename F_ = F, std::enable_if_t<frames::has_lon_lat_v<F_>, int> = 0>
  qtty::Degree longitude() const {
    return azimuth_;
  }
 
  template <typename F_ = F, std::enable_if_t<frames::has_lon_lat_v<F_>, int> = 0>
  qtty::Degree latitude() const {
    return polar_;
  }


  siderust_spherical_dir_t to_c() const { return {polar_.value(), azimuth_.value(), frame_id()}; }
 
  static Direction from_c(const siderust_spherical_dir_t &c) {
    return Direction(qtty::Degree(c.azimuth_deg), qtty::Degree(c.polar_deg));
  }


  qtty::Degree angular_separation(const Direction &other) const {
    constexpr double DEG2RAD = constants::pi / 180.0;
    constexpr double RAD2DEG = 180.0 / constants::pi;
    const double az1 = azimuth_.value() * DEG2RAD;
    const double po1 = polar_.value() * DEG2RAD;
    const double az2 = other.azimuth_.value() * DEG2RAD;
    const double po2 = other.polar_.value() * DEG2RAD;
    const double x =
        std::cos(po1) * std::sin(po2) - std::sin(po1) * std::cos(po2) * std::cos(az2 - az1);
    const double y = std::cos(po2) * std::sin(az2 - az1);
    const double z =
        std::sin(po1) * std::sin(po2) + std::cos(po1) * std::cos(po2) * std::cos(az2 - az1);
    return qtty::Degree(std::atan2(std::sqrt(x * x + y * y), z) * RAD2DEG);
  }

  template <typename Target>
  std::enable_if_t<frames::has_frame_transform_v<F, Target>, Direction<Target>>
  to_frame(const Time<TT, JD> &jd) const {
    if constexpr (std::is_same_v<F, Target>) {
      return Direction<Target>(azimuth_, polar_);
    } else {
      siderust_spherical_dir_t out;
      check_status(siderust_spherical_dir_transform_frame(
                       polar_.value(), azimuth_.value(), frames::FrameTraits<F>::ffi_id,
                       frames::FrameTraits<Target>::ffi_id, jd.value(), &out),
                   "Direction::to_frame");
      return Direction<Target>::from_c(out);
    }
  }

  template <typename Target>
  std::enable_if_t<frames::has_frame_transform_v<F, Target>, Direction<Target>>
  to_frame_with(const Time<TT, JD> &jd, const AstroContext &ctx) const {
    if constexpr (std::is_same_v<F, Target>) {
      return Direction<Target>(azimuth_, polar_);
    } else {
      siderust_spherical_dir_t out;
      detail::OwnedFfiContext fctx(ctx);
      check_status(siderust_spherical_dir_transform_frame_with_context(
                       polar_.value(), azimuth_.value(), frames::FrameTraits<F>::ffi_id,
                       frames::FrameTraits<Target>::ffi_id, jd.value(), fctx.get(), &out),
                   "Direction::to_frame_with");
      return Direction<Target>::from_c(out);
    }
  }

  template <typename Target>
  auto to(const Time<TT, JD> &jd) const -> decltype(this->template to_frame<Target>(jd)) {
    return to_frame<Target>(jd);
  }

  cartesian::Direction<F> to_cartesian() const {
    constexpr double DEG2RAD = constants::pi / 180.0;
    const double az = azimuth_.value() * DEG2RAD;
    const double po = polar_.value() * DEG2RAD;
    const double cp = std::cos(po);
    return cartesian::Direction<F>(std::cos(az) * cp, std::sin(az) * cp, std::sin(po));
  }

  template <typename F_ = F>
  std::enable_if_t<frames::has_horizontal_transform_v<F_>, Direction<frames::Horizontal>>
  to_horizontal(const Time<TT, JD> &jd, const Geodetic &observer) const {
    siderust_spherical_dir_t out;
    check_status(siderust_spherical_dir_to_horizontal(polar_.value(), azimuth_.value(),
                                                      frames::FrameTraits<F>::ffi_id, jd.value(),
                                                      observer.to_c(), &out),
                 "Direction::to_horizontal");
    return Direction<frames::Horizontal>::from_c(out);
  }

  template <typename F_ = F>
  std::enable_if_t<frames::has_horizontal_transform_v<F_>, Direction<frames::Horizontal>>
  to_horizontal_with(const Time<TT, JD> &jd, const Geodetic &observer,
                     const AstroContext &ctx) const {
    siderust_spherical_dir_t out;
    detail::OwnedFfiContext fctx(ctx);
    check_status(siderust_spherical_dir_to_horizontal_precise_with_context(
                     polar_.value(), azimuth_.value(), frames::FrameTraits<F>::ffi_id, jd.value(),
                     jd.value(), observer.to_c(), fctx.get(), &out),
                 "Direction::to_horizontal_with");
    return Direction<frames::Horizontal>::from_c(out);
  }

  template <typename F_ = F>
  std::enable_if_t<frames::has_horizontal_transform_v<F_>, Direction<frames::Horizontal>>
  to_horizontal_precise(const Time<TT, JD> &jd_tt, const Time<UT1, JD> &jd_ut1,
                        const Geodetic &observer) const {
    siderust_spherical_dir_t out;
    check_status(siderust_spherical_dir_to_horizontal_precise(
                     polar_.value(), azimuth_.value(), frames::FrameTraits<F>::ffi_id,
                     jd_tt.value(), jd_ut1.value(), observer.to_c(), &out),
                 "Direction::to_horizontal_precise");
    return Direction<frames::Horizontal>::from_c(out);
  }
};

template <typename C, typename F, typename U = qtty::Meter> struct Position {
  static_assert(frames::is_frame_v<F>, "F must be a valid frame tag");
  static_assert(centers::is_center_v<C>, "C must be a valid center tag");
 
private:
  qtty::Degree azimuth_;
  qtty::Degree polar_;
  U dist_;
 
public:
  Position(qtty::Degree azimuth, qtty::Degree polar, U distance)
      : azimuth_(azimuth), polar_(polar), dist_(distance) {}
 
  Position(const Direction<F> &dir, U distance)
      : azimuth_(dir.azimuth()), polar_(dir.polar()), dist_(distance) {}

  Direction<F> direction() const { return Direction<F>(azimuth_, polar_); }

  template <typename F_ = F, std::enable_if_t<frames::has_ra_dec_v<F_>, int> = 0>
  qtty::Degree ra() const {
    return azimuth_;
  }
 
  template <typename F_ = F, std::enable_if_t<frames::has_ra_dec_v<F_>, int> = 0>
  qtty::Degree dec() const {
    return polar_;
  }
 
  template <typename F_ = F, std::enable_if_t<frames::has_az_alt_v<F_>, int> = 0>
  qtty::Degree az() const {
    return azimuth_;
  }
 
  template <typename F_ = F, std::enable_if_t<frames::has_az_alt_v<F_>, int> = 0>
  qtty::Degree al() const {
    return polar_;
  }
 
  template <typename F_ = F, std::enable_if_t<frames::has_az_alt_v<F_>, int> = 0>
  qtty::Degree alt() const {
    return polar_;
  }
 
  template <typename F_ = F, std::enable_if_t<frames::has_lon_lat_v<F_>, int> = 0>
  qtty::Degree lon() const {
    return azimuth_;
  }
 
  template <typename F_ = F, std::enable_if_t<frames::has_lon_lat_v<F_>, int> = 0>
  qtty::Degree lat() const {
    return polar_;
  }

 
  static constexpr siderust_frame_t frame_id() { return frames::FrameTraits<F>::ffi_id; }
  static constexpr siderust_center_t center_id() { return centers::CenterTraits<C>::ffi_id; }
 
  U distance() const { return dist_; }
  cartesian::Position<C, F, U> to_cartesian() const;

  template <typename Target>
  std::enable_if_t<frames::has_frame_transform_v<F, Target>, Position<C, Target, U>>
  to_frame(const Time<TT, JD> &jd) const;
 
  template <typename Target>
  std::enable_if_t<frames::has_frame_transform_v<F, Target>, Position<C, Target, U>>
  to_frame_with(const Time<TT, JD> &jd, const AstroContext &ctx) const;

  template <typename Target>
  auto to(const Time<TT, JD> &jd) const -> decltype(this->template to_frame<Target>(jd)) {
    return to_frame<Target>(jd);
  }

  qtty::Degree angular_separation(const Position &other) const {
    return direction().angular_separation(other.direction());
  }
 
  U distance_to(const Position &other) const {
    using std::sqrt;
    // Values in underlying unit (e.g. meters)
    const double r = dist_.value();
    const double s = other.dist_.value();
 
    // convert degrees to radians
    constexpr double DEG2RAD = constants::pi / 180.0;
    const double a1 = azimuth_.value() * DEG2RAD;
    const double p1 = polar_.value() * DEG2RAD;
    const double a2 = other.azimuth_.value() * DEG2RAD;
    const double p2 = other.polar_.value() * DEG2RAD;
 
    // dot product of unit direction vectors (spherical -> cartesian)
    double cos_p1 = std::cos(p1);
    double cos_p2 = std::cos(p2);
    double dot = cos_p1 * cos_p2 * std::cos(a1 - a2) + std::sin(p1) * std::sin(p2);
    if (dot > 1.0)
      dot = 1.0;
    if (dot < -1.0)
      dot = -1.0;
 
    double d = std::sqrt(r * r + s * s - 2.0 * r * s * dot);
    return U(d);
  }
};
 
// ============================================================================
// Stream operators
// ============================================================================

template <typename F, std::enable_if_t<frames::has_ra_dec_v<F>, int> = 0>
inline std::ostream &operator<<(std::ostream &os, const Direction<F> &dir) {
  coordinates::detail::write_frame<F>(os);
  return os << " (ra=" << dir.ra() << ", dec=" << dir.dec() << ')';
}

template <typename F, std::enable_if_t<frames::has_az_alt_v<F>, int> = 0>
inline std::ostream &operator<<(std::ostream &os, const Direction<F> &dir) {
  coordinates::detail::write_frame<F>(os);
  return os << " (az=" << dir.az() << ", alt=" << dir.alt() << ')';
}

template <typename F, std::enable_if_t<frames::has_lon_lat_v<F>, int> = 0>
inline std::ostream &operator<<(std::ostream &os, const Direction<F> &dir) {
  coordinates::detail::write_frame<F>(os);
  return os << " (lon=" << dir.lon() << ", lat=" << dir.lat() << ')';
}

template <typename C, typename F, typename U, std::enable_if_t<frames::has_ra_dec_v<F>, int> = 0>
inline std::ostream &operator<<(std::ostream &os, const Position<C, F, U> &pos) {
  coordinates::detail::write_center_frame<C, F>(os);
  return os << " (ra=" << pos.ra() << ", dec=" << pos.dec() << ", r=" << pos.distance() << ')';
}

template <typename C, typename F, typename U, std::enable_if_t<frames::has_az_alt_v<F>, int> = 0>
inline std::ostream &operator<<(std::ostream &os, const Position<C, F, U> &pos) {
  coordinates::detail::write_center_frame<C, F>(os);
  return os << " (az=" << pos.az() << ", alt=" << pos.alt() << ", r=" << pos.distance() << ')';
}

template <typename C, typename F, typename U, std::enable_if_t<frames::has_lon_lat_v<F>, int> = 0>
inline std::ostream &operator<<(std::ostream &os, const Position<C, F, U> &pos) {
  coordinates::detail::write_center_frame<C, F>(os);
  return os << " (lon=" << pos.lon() << ", lat=" << pos.lat() << ", r=" << pos.distance() << ')';
}
 
} // namespace spherical
} // namespace siderust