Orbit.hxx

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
 * Copyright (c) 2025, Davide Stocco and Enrico Bertolazzi.                  *
 *                                                                           *
 * The Orbit project is distributed under the GNU GPLv3.                     *
 *                                                                           *
 * Davide Stocco                                           Enrico Bertolazzi *
 * University of Trento                                 University of Trento *
 * e-mail: davide.stocco@unitn.it         e-mail: enrico.bertolazzi@unitn.it *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#pragma once
 
#ifndef ASTRO_ORBIT_HXX
#define ASTRO_ORBIT_HXX
 
namespace Astro
{
 
  /*\
   |    ___       _     _ _
   |   / _ \ _ __| |__ (_) |_
   |  | | | | '__| '_ \| | __|
   |  | |_| | |  | |_) | | |_
   |   \___/|_|  |_.__/|_|\__|
   |
  \*/

  class Orbit
  {
    using Cartesian    = OrbitalElements::Cartesian;
    using Keplerian    = OrbitalElements::Keplerian;
    using Equinoctial  = OrbitalElements::Equinoctial;
    using Quaternionic = OrbitalElements::Quaternionic;
    using Anomaly      = OrbitalElements::Anomaly;
 
  protected:
    Cartesian    m_cart;                      
    Keplerian    m_kepl;                      
    Equinoctial  m_equi;                      
    Quaternionic m_quat;                      
    Anomaly      m_anom;                      
    Type         m_type{Type::UNDEFINED};     
    Factor       m_factor{Factor::UNDEFINED}; 
    Real         m_mu{QUIET_NAN};             
 
  public:
    Orbit() {}

    Orbit(Orbit const &) = default;

    Orbit(Orbit &&) = default;

    Orbit & operator=(const Orbit &) = default;

    Orbit & operator=(Orbit &&) = default;

    Cartesian const & cartesian() const {return this->m_cart;}

    void set_cartesian(Real r_x, Real r_y, Real r_z, Real v_x, Real v_y, Real v_z)
    {
      this->m_cart.r << r_x, r_y, r_z;
      this->m_cart.v << v_x, v_y, v_z;
      OrbitalElements::cartesian_to_keplerian(this->m_cart, this->m_mu, this->m_kepl);
      OrbitalElements::keplerian_to_equinoctial(this->m_kepl, this->m_factor, this->m_equi);
      // TODO: OrbitalElements::cartesian_to_quaternionic(this->m_cart, this->m_quat);
      this->type(this->m_kepl.e);
    }

    void cartesian(Vector3 const & t_r, Vector3 const & t_v)
    {
      this->m_cart.r = t_r;
      this->m_cart.v = t_v;
      OrbitalElements::cartesian_to_keplerian(this->m_cart, this->m_mu, this->m_kepl);
      OrbitalElements::keplerian_to_equinoctial(this->m_kepl, this->m_factor, this->m_equi);
      // TODO: OrbitalElements::cartesian_to_quaternionic(this->m_cart, this->m_quat);
      this->type(this->m_kepl.e);
    }

    void cartesian(Cartesian const & t_cart)
    {
      this->m_cart = t_cart;
      OrbitalElements::cartesian_to_keplerian(this->m_cart, this->m_mu, this->m_kepl);
      OrbitalElements::keplerian_to_equinoctial(this->m_kepl, this->m_factor, this->m_equi);
      // TODO: OrbitalElements::cartesian_to_quaternionic(this->m_cart, this->m_quat);
      this->type(this->m_kepl.e);
    }

    Keplerian const & keplerian() const {return this->m_kepl;}

    void keplerian(Real t_a, Real t_e, Real t_i, Real t_Omega, Real t_omega)
    {
      this->m_kepl.a     = t_a;
      this->m_kepl.e     = t_e;
      this->m_kepl.i     = t_i;
      this->m_kepl.Omega = t_Omega;
      this->m_kepl.omega = t_omega;
      OrbitalElements::keplerian_to_cartesian(this->m_kepl, this->m_anom, this->m_mu, this->m_cart);
      OrbitalElements::keplerian_to_equinoctial(this->m_kepl, this->m_factor, this->m_equi);
      // TODO: OrbitalElements::keplerian_to_quaternionic(this->m_kepl, this->m_quat);
      this->type(this->m_kepl.e);
    }

    void keplerian(Keplerian const & t_kepl)
    {
      this->m_kepl = t_kepl;
      OrbitalElements::keplerian_to_cartesian(this->m_kepl, this->m_anom, this->m_mu, this->m_cart);
      OrbitalElements::keplerian_to_equinoctial(this->m_kepl, this->m_factor, this->m_equi);
      // TODO: OrbitalElements::keplerian_to_quaternionic(this->m_kepl, this->m_quat);
      this->type(this->m_kepl.e);
    }

    Equinoctial const & equinoctial() const {return this->m_equi;}

    void set_equinoctial(Real t_p, Real t_f, Real t_g, Real t_h, Real t_k)
    {
      this->m_equi.p = t_p;
      this->m_equi.f = t_f;
      this->m_equi.g = t_g;
      this->m_equi.h = t_h;
      this->m_equi.k = t_k;
      OrbitalElements::equinoctial_to_keplerian(this->m_equi, this->m_kepl);
      OrbitalElements::equinoctial_to_cartesian(this->m_equi, this->m_anom, this->m_mu, this->m_cart);
      // TODO: OrbitalElements::equinoctial_to_quaternionic(this->m_equi, this->m_quat);
      this->type(this->m_kepl.e);
    }

    void set_equinoctial(Equinoctial const & t_equi)
    {
      this->m_equi = t_equi;
      OrbitalElements::equinoctial_to_keplerian(this->m_equi, this->m_kepl);
      OrbitalElements::equinoctial_to_cartesian(this->m_equi, this->m_anom, this->m_mu, this->m_cart);
      // TODO: OrbitalElements::equinoctial_to_quaternionic(this->m_equi, this->m_quat);
      this->type(this->m_kepl.e);
    }

    Quaternionic const & quaternionic() const {return this->m_quat;}

    void set_quaternionic(Real t_q_1, Real t_q_2, Real t_q_3, Real t_q_4)
    {
      this->m_quat.q.x() = t_q_1;
      this->m_quat.q.y() = t_q_2;
      this->m_quat.q.z() = t_q_3;
      this->m_quat.q.w() = t_q_4;
      // TODO: OrbitalElements::quaternionic_to_keplerian(this->m_quat, this->m_kepl);
      // TODO: OrbitalElements::quaternionic_to_equinoctial(this->m_quat, this->m_factor, this->m_equi);
      // TODO: OrbitalElements::quaternionic_to_cartesian(this->m_quat, this->m_anom, this->m_mu, this->m_cart);
      this->type(this->m_kepl.e);
    }

    void set_quaternionic(Quaternion const & t_quat)
    {
      this->m_quat.q = t_quat;
      // TODO: OrbitalElements::quaternionic_to_keplerian(this->m_quat, this->m_kepl);
      // TODO: OrbitalElements::quaternionic_to_equinoctial(this->m_quat, this->m_factor, this->m_equi);
      // TODO: OrbitalElements::quaternionic_to_cartesian(this->m_quat, this->m_anom, this->m_mu, this->m_cart);
      this->type(this->m_kepl.e);
    }

    Anomaly const & anomaly() const {return this->m_anom;}

    void anomaly(Anomaly const & t_anom)
    {
      ASTRO_ASSERT(t_anom.sanity_check(),
        "Astro::Orbit::anomaly(...): invalid orbital anomalies detected.");
      this->m_anom = t_anom;
    }

    Type type() const {return this->m_type;}

    void type(Real e)
    {
      if (e > 0.0 && e < 1.0) {
        this->m_type = Type::ELLIPTIC;
      } else if (e == 1.0) {
        this->m_type = Type::PARABOLIC;
      } else if (e > 1.0) {
        this->m_type = Type::HYPERBOLIC;
      } else {
        ASTRO_ERROR("Astro::Orbit::type(...): invalid eccentricity detected.");
      }
    }

    Factor factor() const {return this->m_factor;}

    void factor(Factor t_factor) {this->m_factor = t_factor;}

    Real mu() const {return this->m_mu;}

    void mu(Real t_mu) {this->m_mu = t_mu;}

    std::string info() const {
      std::ostringstream os;
      os <<
        "Cartesian orbit parameters:" << std::endl <<
        this->m_cart.info() << std::endl <<
        "Keplerian orbit parameters:" << std::endl <<
        this->m_kepl.info() << std::endl <<
        "Equinoctial orbit parameters:" << std::endl <<
        this->m_equi.info() << std::endl <<
        "Quaternionic orbit parameters:" << std::endl <<
        this->m_quat.info() << std::endl <<
        "Orbital anomalies:" << std::endl <<
        this->m_anom.info() << std::endl <<
        "Orbit:" << std::endl <<
        "µ: grav. const. = " << this->m_mu << " (UA³/day²)" << std::endl <<
        "type            = " << (this->m_type == Type::UNDEFINED ? "undefined" :
          (this->m_type == Type::ELLIPTIC ? "elliptic" :
          (this->m_type == Type::PARABOLIC ? "parabolic" :
          (this->m_type == Type::HYPERBOLIC ? "hyperbolic" : "undefined")))) << std::endl <<
        "factor          = " << (this->m_factor == Factor::UNDEFINED ? "undefined" :
          (this->m_factor == Factor::POSIGRADE ? "posigrade" :
          (this->m_factor == Factor::RETROGRADE ? "retrograde" : "undefined"))) << std::endl;
        return os.str();
    }

    void info(std::ostream & os) {os << this->info();}

    void reset() {
      this->m_cart.reset();
      this->m_kepl.reset();
      this->m_equi.reset();
      this->m_quat.reset();
      this->m_anom.reset();
      this->m_type   = Type::UNDEFINED;
      this->m_factor = Factor::UNDEFINED;
    }

    Rotation keplerian_to_reference() const
    {
      Real s_Omega{std::sin(this->m_kepl.Omega)};
      Real c_Omega{std::cos(this->m_kepl.Omega)};
      Real s_omega{std::sin(this->m_kepl.omega)};
      Real c_omega{std::cos(this->m_kepl.omega)};
      Real s_i{std::sin(this->m_kepl.i)};
      Real c_i{std::cos(this->m_kepl.i)};
 
      Rotation R;
      R <<
        c_Omega*c_omega - s_Omega*s_omega*c_i, -c_Omega*s_omega - s_Omega*c_omega*c_i, s_Omega*s_i,
        s_Omega*c_omega + c_Omega*s_omega*c_i, -s_Omega*s_omega + c_Omega*c_omega*c_i, -c_Omega*s_i,
        s_i*s_omega,                           s_i*c_omega,                            c_i;
      return R;
    }

    Rotation equinoctial_to_reference() const
    {
      Real h{this->m_equi.h};
      Real k{this->m_equi.k};
      Real h2{h*h};
      Real k2{k*k};
      Real hk{h*k};
      Real I{static_cast<Real>(this->m_factor)};
 
      Rotation R;
      R <<
        1.0-h2+k2, 2.0*I*hk,      2.0*h,
        2.0*hk,    I*(1.0+h2-k2), -2.0*k,
        -2.0*I*h,  2.0*k,         I*(1.0-h2-k2);
      return R / (1.0+h2+k2);
    }

    Rotation cartesian_to_frenet_rtn() const
    {
      // Initialize the Frenet-Serret frame
      Rotation rtn;
 
      // Compute the radial vector
      rtn.col(0) = this->m_cart.r;
      rtn.col(0).normalize();
 
      // Compute the tangential vector
      rtn.col(1) = rtn.col(0).cross(this->m_cart.v);
      rtn.col(1).normalize();
 
      // Compute the binormal vector
      rtn.col(2) = rtn.col(0).cross(rtn.col(1));
 
      // Return the Frenet-Serret frame
      return rtn;
    }

    Rotation equinoctial_to_frenet_rtn() const
    {
      Real h{this->m_equi.h};
      Real k{this->m_equi.k};
      Real c_L{std::cos(this->m_anom.L)};
      Real s_L{std::sin(this->m_anom.L)};
      Real h2{h*h};
      Real k2{k*k};
      Real hk{h*k};
      Real bf{1.0+h2+k2};
      Real I{static_cast<Real>(this->m_factor)};
 
      Rotation R;
      R <<
        ((h2-k2+1)*c_L+2.0*I*hk*s_L)/bf,   ((k2-h2-1.0)*s_L+2.0*hk*I*c_L)/bf, 2.0*k/bf,
        (2.0*hk*c_L-I*(h2-k2-1.0)*s_L)/bf, (I*(1.0+k2-h2)*c_L-2.0*hk*s_L)/bf, -2.0*h/bf,
        (2.0*(h*s_L-c_L*k*I))/bf,          (2.0*k*I*s_L+2.0*h*c_L)/bf,        I*(1-h2-k2)/bf;
      return R;
    }

    Vector3 cartesian_rtn_to_xyz(Vector3 const & vec) const
    {
      return this->cartesian_to_frenet_rtn() * vec;
    }

    Vector3 equinoctial_rtn_to_xyz(Vector3 const & vec) const
    {
      return this->equinoctial_to_frenet_rtn() * vec;
    }

    Vector3 cartesian_xyz_to_rtn(Vector3 const & vec) const
    {
      return this->cartesian_to_frenet_rtn().transpose() * vec;
    }

    Vector3 equinoctial_xyz_to_rtn(Vector3 const & vec) const
    {
      return this->equinoctial_to_frenet_rtn().transpose() * vec;
    }
 
  }; // class Orbit
 
} // namespace Astro
 
#endif // ASTRO_ORBIT_HXX