Utilities.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
 * Copyright (c) 2025, Davide Stocco and Enrico Bertolazzi.                                      *
 *                                                                                               *
 * The Astro project is distributed under the BSD 2-Clause License.                              *
 *                                                                                               *
 * 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_UTILITIES_HH
#define ASTRO_UTILITIES_HH
 
#include "Astro.hh"
 
namespace Astro {
 
  /*\
   |   ____
   |  |  _ \ _____      _____ _ __
   |  | |_) / _ \ \ /\ / / _ \ '__|
   |  |  __/ (_) \ V  V /  __/ |
   |  |_|   \___/ \_/\_/ \___|_|
   |
  \*/

  Real Power2(Real x) {return x*x;}

  template<typename Real>
  Real Power3(Real x) {return x*x*x;}

  Real Power4(Real x) {return x*x*x*x;}
 
  /*\
   |      _                _
   |     / \   _ __   __ _| | ___
   |    / _ \ | '_ \ / _` | |/ _ \
   |   / ___ \| | | | (_| | |  __/
   |  /_/   \_\_| |_|\__, |_|\___|
   |                 |___/
  \*/
 
  static Real const PI         = Real(3.141592653589793238462643383279502884197); 
  static Real const PIMUL2     = Real(6.283185307179586476925286766559005768394); 
  static Real const PIDIV2     = Real(1.570796326794896619231321691639751442098); 
  static Real const DEG_TO_RAD = Real(0.017453292519943295769236907684886127134); 
  static Real const RAD_TO_DEG = Real(57.29577951308232087679815481410517033240); 

  Real Deg_To_Rad(Real x) {return DEG_TO_RAD*x;}

  Real Rad_To_Deg(Real x) {return RAD_TO_DEG*x;}

  Real AngleInRange(Real x)
  {
    x = std::fmod(x, PIMUL2);
    while (x < Real(0.0)) {x += PIMUL2;}
    while (x > PIMUL2) {x -= PIMUL2;}
    return x;
  }

  Real AngleInRangeSym(Real x)
  {
    x = std::fmod(x, PIMUL2);
    while (x < -PI) {x += PIMUL2;}
    while (x > PI) {x -= PIMUL2;}
    return x;
  }
 
  /*\
   |   ____  _     _
   |  |  _ \‍(_)___| |_ __ _ _ __   ___ ___
   |  | | | | / __| __/ _` | '_ \ / __/ _ \
   |  | |_| | \__ \ || (_| | | | | (_|  __/
   |  |____/|_|___/\__\__,_|_| |_|\___\___|
   |
  \*/
 
  static Real const AU_TO_KM{1.49597870707e+08}; 
  static Real const KM_TO_AU{1.0/AU_TO_KM};      
  static Real const KM_TO_M{1000.0};             
  static Real const M_TO_KM{1.0/KM_TO_M};        
  static Real const AU_TO_M{AU_TO_KM*KM_TO_M};   
  static Real const M_TO_AU{AU_TO_KM*KM_TO_M};   
 
  // Time units
  static Real const DAY_TO_SEC{86400.0};        
  static Real const SEC_TO_DAY{1.0/DAY_TO_SEC}; 
 
  static Real const KG_M_SEC2_TO_KG_AU_DAY2{M_TO_AU/(SEC_TO_DAY*SEC_TO_DAY)}; // Kg * m/s^2 => Kg * UA / day^2
 
  static Real const gravity_kg_m_s2{9.80665};
  static Real const gravity_kg_AU_DAY2{gravity_kg_m_s2*KG_M_SEC2_TO_KG_AU_DAY2};
 
 
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
  Real KM_To_M(Real x) {return x * KM_TO_M;}
  Real M_To_KM(Real x) {return x * M_TO_KM;}
 
  Real KM_To_AU(Real x) {return x * KM_TO_AU;}
  Real AU_To_KM(Real x) {return x * AU_TO_KM;}
 
  Real M_To_AU(Real x) {return x * M_TO_AU;}
  Real AU_To_M(Real x) {return x * AU_TO_M;}
 
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
  Real KM_S_To_M_S(Real x) {return x * KM_TO_M;}
  Real M_S_To_KM_S(Real x) {return x * M_TO_KM;}
 
  Real KM_S_To_AU_S(Real x) {return x * (KM_TO_AU/SEC_TO_DAY);}
  Real AU_S_To_KM_S(Real x) {return x * (AU_TO_KM*DAY_TO_SEC);}
 
  Real M_S_To_AU_S(Real x) {return x * (M_TO_AU/SEC_TO_DAY);}
  Real AU_S_To_M_S(Real x) {return x * (AU_TO_M*DAY_TO_SEC);}
 
  Real M_S_To_AU_DAY(Real x) {return x * (M_TO_AU/DAY_TO_SEC);}
  Real AU_DAY_To_M_S(Real x) {return x * (AU_TO_M*DAY_TO_SEC);}
 
  Real KM_S_To_AU_DAY(Real x) {return x * (KM_TO_AU/DAY_TO_SEC);}
  Real AU_DAY_To_KM_S(Real x) {return x * (AU_TO_KM*DAY_TO_SEC);}
 
  Real M_S_To_KM_DAY(Real x) {return x * (M_TO_KM/DAY_TO_SEC);}
  Real KM_DAY_To_M_S(Real x) {return x * (KM_TO_M*DAY_TO_SEC);}
 
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
  Real AU_DAY2_To_KM_S2(Real x) {return x * (AU_TO_KM/(DAY_TO_SEC*DAY_TO_SEC));}
 
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
  Real AU_DAY3_To_KM_S3(Real x) {return x * (AU_TO_KM/(DAY_TO_SEC*DAY_TO_SEC*DAY_TO_SEC));}
 
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
  Real KM3_S2_To_KM3_DAY2(Real x) {return x * (1.0/(SEC_TO_DAY*SEC_TO_DAY));}
 
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
  Real KM3_S2_To_AU3_DAY2(Real x) {return x * (1.0/(AU_TO_KM*AU_TO_KM*AU_TO_KM)/(SEC_TO_DAY*SEC_TO_DAY));}
 
  /*\
   |   _____ _
   |  |_   _(_)_ __ ___   ___
   |    | | | | '_ ` _ \ / _ \
   |    | | | | | | | | |  __/
   |    |_| |_|_| |_| |_|\___|
   |
  \*/

  Real JulianDate(Integer year, Integer month, Integer day)
  {
    if (month <= 2)
    {
      year -= 1;
      month += 12;
    }
    Integer A{year / 100};
    Integer B{2 - A + (A / 4)};
    Real JD{std::floor(365.25 * (year + 4716))
            + std::floor(30.6001 * (month + 1))
            + day + B - 1524.5};
    return JD;
  }
 
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  Real ModifiedJulianDate(Integer year, Integer month, Integer day)
  {
    return JulianDate(year, month, day) - Real{2400000.5};
  }
 
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  void JD_To_Date(Real JD, Integer & year, Integer & month, Integer & day)
  {
    Integer Z{Integer(std::floor(JD + 0.5))};
    Real F{(JD + 0.5) - Z};
    Integer A{Z};
    if (Z >= 2299161)
    {
      Integer alpha{Integer((Z - 1867216.25) / 36524.25)};
      A = Z + 1 + alpha - alpha / 4;
    }
    Integer B{A + 1524};
    Integer C{Integer((B - 122.1) / 365.25)};
    Integer D{Integer(365.25 * C)};
    Integer E{Integer((B - D) / 30.6001)};
    Real dayD{B - D - std::floor(30.6001 * E) + F};
    day = Integer(dayD);
    if (E < 14)
      month = E - 1;
    else
      month = E - 13;
    if (month > 2)
      year = C - 4716;
    else
      year = C - 4715;
  }
 
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  void MJD_To_Date(Real MJD, Integer & year, Integer & month, Integer & day)
  {
    JD_To_Date(MJD + Real{2400000.5}, year, month, day);
  }
 
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
} // namespace Astro
 
#endif // ASTRO_UTILITIES_HH