Body.hxx

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
 * Copyright (c) 2025, Davide Stocco and Enrico Bertolazzi.                  *
 *                                                                           *
 * The Astro 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_BODY_HXX
#define ASTRO_BODY_HXX
 
namespace Astro
{
 
  /*\
   |   ____            _
   |  | __ )  ___   __| |_   _
   |  |  _ \ / _ \ / _` | | | |
   |  | |_) | (_) | (_| | |_| |
   |  |____/ \___/ \__,_|\__, |
   |                     |___/
  \*/

  class Body : public Orbit
  {
    using MatrixA = Eigen::Matrix<Real, 6, 3>; 
    using VectorB = Eigen::Vector<Real, 6>;    
 
    Real m_mass{QUIET_NAN}; 
 
  public:
    Body() {}

    Body(Body const &) = default;

    Body(Body &&) = default;

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

    Body & operator=(Body &&) = default;

    Real mass() const {return this->m_mass;}

    void mass(Real t_mass) {this->m_mass = t_mass;}

    MatrixA equinoctial_eom_A() const
    {
      Real const & p{this->m_equi.p};
      Real const & f{this->m_equi.f};
      Real const & g{this->m_equi.g};
      Real const & h{this->m_equi.h};
      Real const & k{this->m_equi.k};
      Real s_L{std::sin(this->m_anom.L)};
      Real c_L{std::cos(this->m_anom.L)};
 
      Real w{1.0 + (f*c_L+g*s_L)};
      Real s2{1.0 + (h*h+k*k)};
      Real bf{std::sqrt(p/this->m_mu)};
      Real bf1{bf/w};
      Real bf2{(h*s_L-k*c_L)*bf1};
 
      MatrixA A;
      A <<
        0.0,     bf1*2.0*p,           0.0,
        bf*s_L,  bf1*((w+1.0)*c_L+f), -bf2*g,
        -bf*c_L, bf1*((w+1.0)*s_L+g), bf2*f,
        0.0,     0.0,                 bf1*s2*c_L/2.0,
        0.0,     0.0,                 bf1*s2*s_L/2.0,
        0.0,     0.0,                 bf2;
      return A;
    }

    VectorB equinoctial_eom_b() const
    {
      Real const & p{this->m_equi.p};
      Real const & f{this->m_equi.f};
      Real const & g{this->m_equi.g};
      Real w{1.0 + (f*std::cos(this->m_anom.L) + g*std::sin(this->m_anom.L))};
      return VectorB(0.0, 0.0, 0.0, 0.0, 0.0, std::sqrt(p*this->m_mu)*power2(w/p));
    }

    Vector6 cartesian_eom(Vector3 const & thrust_rtn) const
    {
      // Compute the cartesian perturbation
      Vector3 thrust_xyz(this->cartesian_rtn_to_xyz(thrust_rtn));
      thrust_xyz /= this->m_mass;
 
      // Compute the equations of motion
      Real r_norm{this->m_cart.r.norm()};
      Vector3 r_mur3(this->m_mu/power3(r_norm) * this->m_cart.r);
      return Vector6(
        this->m_cart.v.x(),
        this->m_cart.v.y(),
        this->m_cart.v.z(),
        thrust_xyz.x() - r_mur3.x(),
        thrust_xyz.y() - r_mur3.y(),
        thrust_xyz.z() - r_mur3.z()
      );
    }

    Vector6 cartesian_eom(Real thrust_rad, Real thrust_tan, Real thrust_nor) const
    {
      return this->cartesian_eom(Vector3(thrust_rad, thrust_tan, thrust_nor));
    }

    Vector6 equinoctial_eom(Vector3 const & thrust) const
    {
      Real const & p{this->m_equi.p};
      Real const & f{this->m_equi.f};
      Real const & g{this->m_equi.g};
      Real const & h{this->m_equi.h};
      Real const & k{this->m_equi.k};
      Real const & L{this->m_anom.L};
 
      Real s_L{std::sin(L)};
      Real c_L{std::cos(L)};
 
      Real w{1.0 + (f*c_L+g*s_L)};
      Real s2{1.0 + (h*h+k*k)};
      Real bf{std::sqrt(p/this->m_mu)};
 
      Real C_rad{thrust(0)*bf/this->m_mass};
      Real C_tan{thrust(1)*bf/this->m_mass/w};
      Real C_nor{thrust(2)*bf/this->m_mass/w};
 
      Vector6 res;
      res.setZero();
 
      res(1) += C_rad*s_L;
      res(2) -= C_rad*c_L;
 
      res(0) += C_tan*2.0*p;
      res(1) += C_tan*((w+1)*c_L+f);
      res(2) += C_tan*((w+1)*s_L+g);
 
      res(1) -= C_nor*(h*s_L-k*c_L)*g;
      res(2) += C_nor*(h*s_L-k*c_L)*f;
      res(3) += C_nor*s2*c_L/2.0;
      res(4) += C_nor*s2*s_L/2.0;
      res(5) += C_nor*(h*s_L-k*c_L);
 
      res(5) += std::sqrt(p*this->m_mu)*power2(w/p);
 
      return res;
    }

    Vector6 equinoctial_eom(Real thrust_rad, Real thrust_tan, Real thrust_nor) const
    {
      return this->equinoctial_eom(Vector3(thrust_rad, thrust_tan, thrust_nor));
    }
 
  }; // class Orbit
 
} // namespace Astro
 
#endif // ASTRO_ORBIT_HXX