Planets.hh

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 *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef ASTRO_PLANETS_HH
#define ASTRO_PLANETS_HH
 
#include "Astro/Body.hh"
#include "Astro/Utilities.hh"
 
namespace Astro
{
  namespace Planets
  {
 
    static const Real KM3_S2_TO_AU3_DAY2{(KM_TO_AU*KM_TO_AU*KM_TO_AU)/(SEC_TO_DAY*SEC_TO_DAY)}; 
    static const Real AU3_DAY2_TO_KM3_S2{1.0/KM3_S2_TO_AU3_DAY2}; 
 
    /*\
     |   ____
     |  / ___| _   _ _ __
     |  \___ \| | | | '_ \
     |   ___) | |_| | | | |
     |  |____/ \__,_|_| |_|
     |
    \*/
 
    static Real const Sun_mass_KG{1.9885E30}; 
    static Real const Sun_radius_KM{695700.0}; 
    static Real const Sun_radius_AU{KM_To_AU(Sun_radius_KM)}; 
    static Real const Sun_mu_KM3_S2{1.32712440018E11}; 
    static Real const Sun_mu_AU3_DAY2{KM3_S2_To_AU3_DAY2(Sun_mu_KM3_S2)}; 
 
    /*\
     |   _____           _   _
     |  | ____|__ _ _ __| |_| |__
     |  |  _| / _` | '__| __| '_ \
     |  | |__| (_| | |  | |_| | | |
     |  |_____\__,_|_|   \__|_| |_|
     |
    \*/
 
    static Real const Earth_mass_KG{5.97219E24}; 
    static Real const Earth_radius_KM{6378.1370}; 
    static Real const Earth_radius_AU{KM_To_AU(Earth_radius_KM)}; 
    static Real const Earth_mu_KM3_S2{398600.4418}; 
    static Real const Earth_mu_KM3_DAY2{KM3_S2_To_KM3_DAY2(Earth_mu_KM3_S2)}; 
    static Real const Earth_mu_AU3_DAY2{KM3_S2_To_AU3_DAY2(Earth_mu_KM3_S2)}; 

    struct KeplerianEarth : public OrbitalElements::Keplerian
    {
      KeplerianEarth()
      {
        this->a     = 1.00000011; // Semi-major axis (AU)
        this->e     = 0.01671022; // Eccentricity (-)
        this->i     = Deg_To_Rad(0.00005); // Inclination (rad)
        this->Omega = Deg_To_Rad(-11.26064); // Longitude of the ascending node (rad)
        this->omega = Deg_To_Rad(102.93768193); // Argument of periapsis (rad)
      }
    };

    inline Body Earth()
    {
      Body earth("Earth", Earth_mass_KG, Earth_radius_AU);
      earth.set_orbit().set_factor(Factor::POSIGRADE);
      earth.set_orbit().set_mu(Sun_mu_AU3_DAY2);
      earth.set_orbit().set_keplerian(KeplerianEarth(), 0.0);
      earth.set_epoch(51544.0); // J2000 epoch (MJD 2451545.0 = JD 51544.5)
      earth.set_epoch_anomaly().set_M(Deg_To_Rad(100.46435), earth.orbit().keplerian(), Factor::POSIGRADE);
      return earth;
    }
 
    static constexpr Real EARTH_MAG_MOMENT{7.96e15}; 

    inline Vector3 EarthMagneticFieldDipole(Vector3 const & position)
    {
      // Compute radial distance
      Real r{position.norm()};
      if (r < EPSILON_LOW) {return ZEROS_VEC3;} // avoid division by zero
 
      // Unit vector
      Vector3 r_hat{position/r};
 
      // Magnetic dipole moment aligned with z-axis (simplified)
      Vector3 m(0.0, 0.0, EARTH_MAG_MOMENT);
 
      Vector3 B{(3.0 * r_hat * m.dot(r_hat) - m) / Power3(r)};
 
      return B;
    }
 
    /*\
     |   __  __
     |  |  \/  | ___   ___  _ __
     |  | |\/| |/ _ \ / _ \| '_ \
     |  | |  | | (_) | (_) | | | |
     |  |_|  |_|\___/ \___/|_| |_|
     |
    \*/
 
    static Real const Moon_mass_KG{7.34767309E22}; 
    static Real const Moon_radius_KM{1737.400}; 
    static Real const Moon_radius_AU{Moon_radius_KM*KM_TO_AU}; 
    static Real const Moon_mu_M3S2{4.9048695E12}; 
    static Real const Moon_mu_KM3_S2{Moon_mu_M3S2/1.0E9}; 
    static Real const Moon_mu_AU3_DAY2{Moon_mu_KM3_S2*KM3_S2_TO_AU3_DAY2}; 

    struct KeplerianMoon : public OrbitalElements::Keplerian
    {
      KeplerianMoon()
      {
        this->a     = 0.002569555; // Semi-major axis (AU)
        this->e     = 0.0549006; // Eccentricity (-)
        this->i     = Deg_To_Rad(5.145396); // Inclination (rad)
        this->Omega = Deg_To_Rad(125.1228); // Longitude of the ascending node (rad)
        this->omega = Deg_To_Rad(318.0634); // Argument of periapsis (rad)
      }
    };

    inline Body Moon()
    {
      Body moon("Moon", Moon_mass_KG, Moon_radius_AU);
      moon.set_orbit().set_factor(Factor::POSIGRADE);
      moon.set_orbit().set_mu(Earth_mu_AU3_DAY2);
      moon.set_orbit().set_keplerian(KeplerianMoon(), 0.0);
      moon.set_epoch(51544.0); // J2000 epoch (JD 2451545.0 = MJD 51544.5)
      moon.set_epoch_anomaly().set_M(Deg_To_Rad(134.96340251), moon.orbit().keplerian(), Factor::POSIGRADE);
      return moon;
    }
 
  } // namespace Planets
 
} // namespace Astro
 
#endif // ASTRO_PLANETS_HH