cartesian.hpp

Go to the documentation of this file.

#pragma once
 
#include "../astro_context.hpp"
#include "../centers.hpp"
#include "../ffi_core.hpp"
#include "../frames.hpp"
#include "../time.hpp"
 
#include <qtty/qtty.hpp>
 
#include "detail/stream.hpp"
 
#include <cmath>
#include <ostream>
#include <type_traits>
 
// Forward-declare spherical Position to avoid circular include.
namespace siderust {
namespace spherical {
template <typename C, typename F, typename U> struct Position;
}
} // namespace siderust
 
namespace siderust {
namespace cartesian {
 
template <typename F> struct Direction {
  static_assert(frames::is_frame_v<F>, "F must be a valid frame tag");
 
  double x; 
  double y; 
  double z; 
 
  Direction() : x(0), y(0), z(0) {}
  Direction(double x_, double y_, double z_) : x(x_), y(y_), z(z_) {}
 
  static constexpr siderust_frame_t frame_id() { return frames::FrameTraits<F>::ffi_id; }
 
  double dot(const Direction &other) const { return x * other.x + y * other.y + z * other.z; }
 
  double angle_to(const Direction &other) const {
    double d = dot(other);
    if (d > 1.0)
      d = 1.0;
    if (d < -1.0)
      d = -1.0;
    return std::acos(d);
  }
 
  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>(x, y, z);
    } else {
      siderust_cartesian_pos_t out{};
      check_status(siderust_cartesian_dir_transform_frame(x, y, z, frames::FrameTraits<F>::ffi_id,
                                                          frames::FrameTraits<Target>::ffi_id,
                                                          jd.value(), &out),
                   "cartesian::Direction::to_frame");
      return Direction<Target>(out.x, out.y, out.z);
    }
  }
 
  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>(x, y, z);
    } else {
      siderust_cartesian_pos_t out{};
      detail::OwnedFfiContext fctx(ctx);
      check_status(siderust_cartesian_dir_transform_frame_with_context(
                       x, y, z, frames::FrameTraits<F>::ffi_id, frames::FrameTraits<Target>::ffi_id,
                       jd.value(), fctx.get(), &out),
                   "cartesian::Direction::to_frame_with");
      return Direction<Target>(out.x, out.y, out.z);
    }
  }
 
  template <typename Target>
  auto to(const Time<TT, JD> &jd) const -> decltype(this->template to_frame<Target>(jd)) {
    return to_frame<Target>(jd);
  }
};
 
template <typename F> inline std::ostream &operator<<(std::ostream &os, const Direction<F> &dir) {
  coordinates::detail::write_frame<F>(os);
  return os << " direction (x=" << dir.x << ", y=" << dir.y << ", z=" << dir.z << ')';
}
 
// Forward-declare Position for Displacement operators.
template <typename C, typename F, typename U> struct Position;
 
template <typename F, typename U> struct Displacement {
  static_assert(frames::is_frame_v<F>, "F must be a valid frame tag");
 
  U comp_x; 
  U comp_y; 
  U comp_z; 
 
  Displacement() : comp_x(U(0)), comp_y(U(0)), comp_z(U(0)) {}
 
  Displacement(U x_, U y_, U z_) : comp_x(x_), comp_y(y_), comp_z(z_) {}
 
  Displacement(double x_, double y_, double z_) : comp_x(U(x_)), comp_y(U(y_)), comp_z(U(z_)) {}
 
  U x() const { return comp_x; }
  U y() const { return comp_y; }
  U z() const { return comp_z; }
 
  U magnitude() const {
    using std::sqrt;
    const double vx = comp_x.value();
    const double vy = comp_y.value();
    const double vz = comp_z.value();
    return U(sqrt(vx * vx + vy * vy + vz * vz));
  }
 
  static constexpr siderust_frame_t frame_id() { return frames::FrameTraits<F>::ffi_id; }
 
  Displacement operator+(const Displacement &other) const {
    return Displacement(comp_x + other.comp_x, comp_y + other.comp_y, comp_z + other.comp_z);
  }
 
  Displacement operator-(const Displacement &other) const {
    return Displacement(comp_x - other.comp_x, comp_y - other.comp_y, comp_z - other.comp_z);
  }
 
  Displacement operator-() const { return Displacement(-comp_x, -comp_y, -comp_z); }
 
  Displacement operator*(double scalar) const {
    return Displacement(comp_x * scalar, comp_y * scalar, comp_z * scalar);
  }
 
  template <typename Target>
  std::enable_if_t<frames::has_frame_transform_v<F, Target>, Displacement<Target, U>>
  to_frame(const Time<TT, JD> &jd) const {
    if constexpr (std::is_same_v<F, Target>) {
      return Displacement<Target, U>(comp_x, comp_y, comp_z);
    } else {
      siderust_cartesian_pos_t out{};
      check_status(siderust_cartesian_dir_transform_frame(
                       comp_x.value(), comp_y.value(), comp_z.value(),
                       frames::FrameTraits<F>::ffi_id, frames::FrameTraits<Target>::ffi_id,
                       jd.value(), &out),
                   "cartesian::Displacement::to_frame");
      return Displacement<Target, U>(out.x, out.y, out.z);
    }
  }
 
  template <typename Target>
  std::enable_if_t<frames::has_frame_transform_v<F, Target>, Displacement<Target, U>>
  to_frame_with(const Time<TT, JD> &jd, const AstroContext &ctx) const {
    if constexpr (std::is_same_v<F, Target>) {
      return Displacement<Target, U>(comp_x, comp_y, comp_z);
    } else {
      siderust_cartesian_pos_t out{};
      detail::OwnedFfiContext fctx(ctx);
      check_status(siderust_cartesian_dir_transform_frame_with_context(
                       comp_x.value(), comp_y.value(), comp_z.value(),
                       frames::FrameTraits<F>::ffi_id, frames::FrameTraits<Target>::ffi_id,
                       jd.value(), fctx.get(), &out),
                   "cartesian::Displacement::to_frame_with");
      return Displacement<Target, U>(out.x, out.y, out.z);
    }
  }
};
 
template <typename F, typename U>
inline std::ostream &operator<<(std::ostream &os, const Displacement<F, U> &d) {
  coordinates::detail::write_frame<F>(os);
  return os << " displacement (dx=" << d.x() << ", dy=" << d.y() << ", dz=" << d.z() << ')';
}
 
template <typename F, typename U>
inline Displacement<F, U> operator*(double scalar, const Displacement<F, U> &d) {
  return d * scalar;
}
 
template <typename C, typename F, typename U> 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");
 
  U comp_x; 
  U comp_y; 
  U comp_z; 
 
  Position() : comp_x(U(0)), comp_y(U(0)), comp_z(U(0)) {}
 
  Position(U x_, U y_, U z_) : comp_x(x_), comp_y(y_), comp_z(z_) {}
 
  Position(double x_, double y_, double z_) : comp_x(U(x_)), comp_y(U(y_)), comp_z(U(z_)) {}
 
  U x() const { return comp_x; }
  U y() const { return comp_y; }
  U z() const { return comp_z; }
 
  U distance() const {
    using std::sqrt;
    const double vx = comp_x.value();
    const double vy = comp_y.value();
    const double vz = comp_z.value();
    return U(sqrt(vx * vx + vy * vy + vz * vz));
  }
 
  spherical::Position<C, F, U> to_spherical() const;
 
  U distance_to(const Position &other) const {
    using std::sqrt;
    const double dx = comp_x.value() - other.comp_x.value();
    const double dy = comp_y.value() - other.comp_y.value();
    const double dz = comp_z.value() - other.comp_z.value();
    return U(sqrt(dx * dx + dy * dy + dz * dz));
  }
 
  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; }
 
  siderust_cartesian_pos_t to_c() const {
    return {comp_x.value(), comp_y.value(), comp_z.value(), frame_id(), center_id()};
  }
 
  static Position from_c(const siderust_cartesian_pos_t &c) { return Position(c.x, c.y, c.z); }
 
  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 {
    if constexpr (std::is_same_v<F, Target>) {
      return *this;
    } else {
      siderust_cartesian_pos_t out{};
      check_status(siderust_cartesian_pos_transform_frame(
                       to_c(), frames::FrameTraits<Target>::ffi_id, jd.value(), &out),
                   "cartesian::Position::to_frame");
      return Position<C, Target, U>(out.x, out.y, out.z);
    }
  }
 
  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 {
    if constexpr (std::is_same_v<F, Target>) {
      return *this;
    } else {
      siderust_cartesian_pos_t out{};
      detail::OwnedFfiContext fctx(ctx);
      check_status(siderust_cartesian_pos_transform_frame_with_context(
                       to_c(), frames::FrameTraits<Target>::ffi_id, jd.value(), fctx.get(), &out),
                   "cartesian::Position::to_frame_with");
      return Position<C, Target, U>(out.x, out.y, out.z);
    }
  }
 
  template <typename Target>
  auto to(const Time<TT, JD> &jd) const -> decltype(this->template to_frame<Target>(jd)) {
    return to_frame<Target>(jd);
  }
 
  template <typename TargetC>
  std::enable_if_t<centers::has_center_transform_v<C, TargetC>, Position<TargetC, F, U>>
  to_center(const Time<TT, JD> &jd) const {
    if constexpr (std::is_same_v<C, TargetC>) {
      return *this;
    } else if constexpr (std::is_same_v<F, frames::EclipticMeanJ2000>) {
      // Direct FFI call — shift vectors and position are both in ecliptic.
      siderust_cartesian_pos_t out{};
      check_status(siderust_cartesian_pos_transform_center(
                       to_c(), centers::CenterTraits<TargetC>::ffi_id, jd.value(), &out),
                   "cartesian::Position::to_center");
      return Position<TargetC, F, U>(out.x, out.y, out.z);
    } else {
      // Route through ecliptic so the shift vectors match the frame.
      auto ecl = to_frame<frames::EclipticMeanJ2000>(jd);
      auto shifted = ecl.template to_center<TargetC>(jd);
      return shifted.template to_frame<F>(jd);
    }
  }
 
  template <typename TargetC, typename TargetF>
  std::enable_if_t<frames::has_frame_transform_v<F, TargetF> &&
                       centers::has_center_transform_v<C, TargetC>,
                   Position<TargetC, TargetF, U>>
  transform(const Time<TT, JD> &jd) const {
    auto ecl = to_frame<frames::EclipticMeanJ2000>(jd);
    auto shifted = ecl.template to_center<TargetC>(jd);
    return shifted.template to_frame<TargetF>(jd);
  }
 
  Displacement<F, U> operator-(const Position &other) const {
    return Displacement<F, U>(U(comp_x.value() - other.comp_x.value()),
                              U(comp_y.value() - other.comp_y.value()),
                              U(comp_z.value() - other.comp_z.value()));
  }
 
  Position operator+(const Displacement<F, U> &displacement) const {
    return Position(comp_x + displacement.comp_x, comp_y + displacement.comp_y,
                    comp_z + displacement.comp_z);
  }
 
  Position operator-(const Displacement<F, U> &displacement) const {
    return Position(comp_x - displacement.comp_x, comp_y - displacement.comp_y,
                    comp_z - displacement.comp_z);
  }
 
  U magnitude() const { return distance(); }
};
 
// ============================================================================
// Stream operators
// ============================================================================
 
template <typename C, typename F, typename U>
inline std::ostream &operator<<(std::ostream &os, const Position<C, F, U> &pos) {
  coordinates::detail::write_center_frame<C, F>(os);
  return os << " (x=" << pos.x() << ", y=" << pos.y() << ", z=" << pos.z() << ')';
}
 
} // namespace cartesian
} // namespace siderust