Program Listing for File constants.h¶

↰ Return to documentation for file (core/constants.h)

#pragma once

#include <Eigen/Core>
#include <Eigen/Dense>
#include <autodiff/forward/real.hpp>
#include <autodiff/forward/real/eigen.hpp>
#include <filesystem>
#include <iostream>
#include <limits>
#include <memory>
#include <sstream>
#include <stdexcept>
#include <string>
#include <string_view>
#include <vector>

#include "lupnt/core/definitions.h"
#include "lupnt/core/file.h"

namespace lupnt {

  // Math constants
  static constexpr double PI = 3.14159265358979323846264338327950288419716939937511;
  static constexpr double TWO_PI = 2.0 * PI;
  static constexpr double PI_OVER_TWO = PI / 2.0;
  static constexpr double E = 2.71828182845904523536028747135266249775724709369996;
  static constexpr double EPS = 1.0e-16;

  // Unit system constants ******************************************************
  //
  // LuPNT's historical constants remain SI-valued for backward compatibility.
  // Use UnitSystem and GetPhysicalConstants(...) when a simulation state is
  // expressed in another coherent unit set, e.g. kilometers, seconds, kilograms.
  static constexpr double METER = 1.0;      // [m/m]
  static constexpr double KILOMETER = 1e3;  // [m/km]
  static constexpr double SECOND = 1.0;     // [s/s]
  static constexpr double KILOGRAM = 1.0;   // [kg/kg]

  constexpr double PowInt(double base, int exponent) {
    double result = 1.0;
    int n = exponent < 0 ? -exponent : exponent;
    for (int i = 0; i < n; i++) result *= base;
    return exponent < 0 ? 1.0 / result : result;
  }

  struct UnitSystem {
    double length = METER;   // SI meters per simulation length unit
    double time = SECOND;    // SI seconds per simulation time unit
    double mass = KILOGRAM;  // SI kilograms per simulation mass unit

    constexpr double FromSI(double value, int length_power, int time_power = 0,
                            int mass_power = 0) const {
      return value
             / (PowInt(length, length_power) * PowInt(time, time_power) * PowInt(mass, mass_power));
    }

    constexpr double ToSI(double value, int length_power, int time_power = 0,
                          int mass_power = 0) const {
      return value * PowInt(length, length_power) * PowInt(time, time_power)
             * PowInt(mass, mass_power);
    }

    constexpr double Length(double value_m) const { return FromSI(value_m, 1); }
    constexpr double Area(double value_m2) const { return FromSI(value_m2, 2); }
    constexpr double Volume(double value_m3) const { return FromSI(value_m3, 3); }
    constexpr double Duration(double value_s) const { return FromSI(value_s, 0, 1); }
    constexpr double Mass(double value_kg) const { return FromSI(value_kg, 0, 0, 1); }
    constexpr double Velocity(double value_m_s) const { return FromSI(value_m_s, 1, -1); }
    constexpr double Acceleration(double value_m_s2) const { return FromSI(value_m_s2, 1, -2); }
    constexpr double GravitationalParameter(double value_m3_s2) const {
      return FromSI(value_m3_s2, 3, -2);
    }
    constexpr double AngularVelocity(double value_rad_s) const {
      return FromSI(value_rad_s, 0, -1);
    }
    constexpr double Frequency(double value_hz) const { return FromSI(value_hz, 0, -1); }
    constexpr double Force(double value_n) const { return FromSI(value_n, 1, -2, 1); }
    constexpr double Pressure(double value_pa) const { return FromSI(value_pa, -1, -2, 1); }
    constexpr double AreaPerMass(double value_m2_kg) const { return FromSI(value_m2_kg, 2, 0, -1); }

    constexpr bool operator==(const UnitSystem& other) const {
      return length == other.length && time == other.time && mass == other.mass;
    }
    constexpr bool operator!=(const UnitSystem& other) const { return !(*this == other); }
  };

  static constexpr UnitSystem SI_UNITS{METER, SECOND, KILOGRAM};
  static constexpr UnitSystem M_S_KG_UNITS = SI_UNITS;
  static constexpr UnitSystem KM_S_KG_UNITS{KILOMETER, SECOND, KILOGRAM};

  // Angle conversion
  static constexpr double RAD = PI / 180.0;               // [rad/deg]
  static constexpr double DEG = 180.0 / PI;               // [deg/rad]
  static constexpr double ARCSEC_DEG = 3600.0;            // [arcsec/deg]
  static constexpr double DEG_ARCSEC = 1.0 / ARCSEC_DEG;  // [deg/arcsec]
  static constexpr double RAD_ARCSEC = DEG_ARCSEC * RAD;  // [rad/arcsec]
  static constexpr double ARCSEC_RAD = 1.0 / RAD_ARCSEC;  // [arcsec/rad]

  // Length
  static constexpr double INCH_M = 0.0254;    // [in/m]
  static constexpr double FOOT_M = 0.3048;    // [ft/m]
  static constexpr double MILE_M = 1609.344;  // [mile/m]
  static constexpr double KM_M = 1e-3;        // [km/m]
  static constexpr double M_KM = 1e3;         // [m/km]
  static constexpr double MM_KM = 1.0e6;      // [mm/km]
  static constexpr double KM_MM = 1.0e-6;     // [km/mm]
  static constexpr double MM_M = 1e3;         // [mm/m]
  static constexpr double M_MM = 1e-3;        // [m/mm]
  static constexpr double M_CM = 1e-2;        // [m/cm]
  static constexpr double CM_M = 1e2;         // [cm/m]

  // Time system constants ******************************************************
  static constexpr double SECS_DAY = 86400.0;         // [s/day]
  static constexpr double SECS_HOUR = 3600.0;         // [s/hour]
  static constexpr double SECS_MINUTE = 60.0;         // [s/minute]
  static constexpr double MINS_HOUR = 60.0;           // [min/hour]
  static constexpr double MINS_DAY = 1440.0;          // [min/day]
  static constexpr double HOURS_DAY = 24.0;           // [hour/day]
  static constexpr double DAYS_SEC = 1.0 / SECS_DAY;  // [day/s]
  static constexpr double DAYS_WEEK = 7.0;            // [days/week]
  static constexpr double DAYS_YEAR = 365.25;         // [days/year]
  static constexpr double DAYS_CENTURY = 36525.00;    // [days/century]

  static constexpr double JD_MJD_OFFSET = 2400000.5;  // [days]
  static constexpr double TT_TAI_OFFSET = 32.184;     // [s]
  static constexpr double A1_TAI_OFFSET = 0.0343817;  // [s]

  static constexpr double JD_CCSDS_TAI = 2436203.5;  // [days]
  static constexpr double JD_J2000_TT = 2451545.0;   // [days]
  static constexpr double MJD_CCSDS_TAI = 36203.0;   // [days]
  static constexpr double MJD_J2000_TT = 51544.5;    // [days]

  // Vallado page 194
  static constexpr double MJD_COORDINATE_TAI = 2443144.5 - JD_MJD_OFFSET;  // [days]
  static constexpr double MJD_COORDINATE_TT_TCG_TCB
      = MJD_COORDINATE_TAI + TT_TAI_OFFSET / SECS_DAY;  // [days]

  // IAU/Turyshev (2026, ApJ 997:97) Table 1 relativistic rate constants.
  // L_B, L_G, L_L define: TDB = TCB − L_B·(TCB−T_0) + TDB_0
  //                        TT  = TCG − L_G·(TCG−T_0)
  //                        TL  = TCL − L_L·(TCL−T_0)
  static constexpr double L_B = 1.550519768e-8;   // TCB→TDB rate  (IAU 2006 B3)
  static constexpr double L_G = 6.969290134e-10;  // TCG→TT  rate  (IAU 1997 B1.9, defining)
  static constexpr double L_L = 3.13905e-11;      // TCL→TL  rate  (Turyshev 2026, selenoid)
  static constexpr double L_H = 1.48253624e-8;    // TCL−TCB mean rate
  static constexpr double L_M = 1.485675290e-8;   // TL−TCB  mean rate
  static constexpr double L_EM = 1.7093906e-11;   // TL−TT   mean rate [s/s TDB]
  // DE405 TDB offset: TDB_0 = −65.5 μs  (Turyshev 2026, Table 1)
  static constexpr double TDB_0 = -65.5e-6;  // [s]

  enum class CoordinateScale {
    TCB,
    TDB,
    TCG,
    TT,
    TCL,
    TL,
  };

  constexpr bool AreCoordinateScalesConvertible(CoordinateScale from, CoordinateScale to) {
    if (from == to) return true;
    bool barycentric_pair = (from == CoordinateScale::TCB || from == CoordinateScale::TDB)
                            && (to == CoordinateScale::TCB || to == CoordinateScale::TDB);
    bool geocentric_pair = (from == CoordinateScale::TCG || from == CoordinateScale::TT)
                           && (to == CoordinateScale::TCG || to == CoordinateScale::TT);
    bool lunicentric_pair = (from == CoordinateScale::TCL || from == CoordinateScale::TL)
                            && (to == CoordinateScale::TCL || to == CoordinateScale::TL);
    return barycentric_pair || geocentric_pair || lunicentric_pair;
  }

  constexpr double CoordinateScaleFactor(CoordinateScale scale) {
    switch (scale) {
      case CoordinateScale::TCB:
      case CoordinateScale::TCG:
      case CoordinateScale::TCL: return 1.0;
      case CoordinateScale::TDB: return 1.0 - L_B;
      case CoordinateScale::TT: return 1.0 - L_G;
      case CoordinateScale::TL: return 1.0 - L_L;
    }
    return 1.0;
  }

  constexpr double CoordinateScaleRatio(CoordinateScale from, CoordinateScale to) {
    return AreCoordinateScalesConvertible(from, to)
               ? CoordinateScaleFactor(to) / CoordinateScaleFactor(from)
               : std::numeric_limits<double>::quiet_NaN();
  }

  constexpr double ScaleLengthForCoordinateScale(double value, CoordinateScale from,
                                                 CoordinateScale to) {
    return value * CoordinateScaleRatio(from, to);
  }

  constexpr double ScaleGravitationalParameterForCoordinateScale(double value, CoordinateScale from,
                                                                 CoordinateScale to) {
    return value * CoordinateScaleRatio(from, to);
  }

  inline double CheckedCoordinateScaleRatio(CoordinateScale from, CoordinateScale to) {
    if (!AreCoordinateScalesConvertible(from, to)) {
      throw std::invalid_argument(
          "Coordinate scales are not related by a constant IAU scale factor");
    }
    return CoordinateScaleFactor(to) / CoordinateScaleFactor(from);
  }

  inline double ScaleLengthForCoordinateScaleChecked(double value, CoordinateScale from,
                                                     CoordinateScale to) {
    return value * CheckedCoordinateScaleRatio(from, to);
  }

  inline double ScaleGravitationalParameterForCoordinateScaleChecked(double value,
                                                                     CoordinateScale from,
                                                                     CoordinateScale to) {
    return value * CheckedCoordinateScaleRatio(from, to);
  }

  inline Vec3 ScalePositionForCoordinateScale(const Vec3& r, CoordinateScale from,
                                              CoordinateScale to) {
    return r * CheckedCoordinateScaleRatio(from, to);
  }

  inline Vec6 ScaleStateForCoordinateScale(const Vec6& rv, CoordinateScale from,
                                           CoordinateScale to) {
    Vec6 out = rv;
    out.head(3) = ScalePositionForCoordinateScale(rv.head(3), from, to);
    return out;
  }

  // Coordinate system constants DE440 *******************************************
  static constexpr double GM_SUN = 132712440041.279419e9;          // [m^3/s^2]
  static constexpr double GM_MERCURY = 22031.868551e9;             // [m^3/s^2]
  static constexpr double GM_VENUS = 324858.592000e9;              // [m^3/s^2]
  static constexpr double GM_EARTH = 398600.435507e9;              // [m^3/s^2]
  static constexpr double GM_MOON = 4902.800118e9;                 // [m^3/s^2]
  static constexpr double GM_MARS_SYSTEM = 42828.375816e9;         // [m^3/s^2]
  static constexpr double GM_JUPITER_SYSTEM = 126712764.100000e9;  // [m^3/s^2]
  static constexpr double GM_SATURN_SYSTEM = 37940584.841800e9;    // [m^3/s^2]
  static constexpr double GM_URANUS_SYSTEM = 5794556.400000e9;     // [m^3/s^2]
  static constexpr double GM_NEPTUNE_SYSTEM = 6836527.100580e9;    // [m^3/s^2]
  static constexpr double GM_PLUTO_SYSTEM = 977.000000e9;          // [m^3/s^2]
  static constexpr double GM_MARS = 0.4282837566395650e14;         // [m^3/s^2]
  static constexpr double GM_JUPITER = 0.1267127646799999e17;      // [m^3/s^2]
  static constexpr double GM_SATURN = 0.3794058480000000e16;       // [m^3/s^2]
  static constexpr double GM_URANUS = 0.5794556400000000e15;       // [m^3/s^2]
  static constexpr double GM_NEPTUNE = 0.6836527100580000e15;      // [m^3/s^2]

  static constexpr double GM_CERES = 62.62890e9;   // [m^3/s^2]
  static constexpr double GM_VESTA = 17.288245e9;  // [m^3/s^2]

  static constexpr double D_EARTH_MOON = 384400.0e3;  // [m]
  static constexpr double D_EARTH_EMB = 4671.0e3;     // [m]
  static constexpr double R_EARTH = 6378.137e3;       // [m]
  static constexpr double R_MOON = 1737.4e3;          // [m]
  static constexpr double R_SUN = 696342.0e3;         // [m]
  static constexpr double R_MERCURY = 2439.7e3;       // [m]
  static constexpr double R_VENUS = 6051.8e3;         // [m]
  static constexpr double R_MARS = 3396.0e3;          // [m]
  static constexpr double R_JUPITER = 71492.0e3;      // [m]
  static constexpr double R_SATURN = 60268.0e3;       // [m]
  static constexpr double R_URANUS = 25559.0e3;       // [m]
  static constexpr double R_NEPTUNE = 24764.0e3;      // [m]
  static constexpr double R_PLUTO = 1188.3e3;         // [m]

  static constexpr double OMEGA_EARTH_MOON = 2.6617e-6;             // [rad/s]
  static constexpr double D_MOON_EMB = D_EARTH_MOON - D_EARTH_EMB;  // [m]

  // Flattening
  static constexpr double WGS84_A = 6378.137e3;           // [m]
  static constexpr double WGS84_F = 1.0 / 298.257223563;  // [-]

  static constexpr double SUN_F = 0.0000;        // [-]
  static constexpr double EARTH_F = WGS84_F;     // [-]
  static constexpr double MOON_F = 0.0012;       // [-]
  static constexpr double MERCURY_F = 0.0009;    // [-]
  static constexpr double VENUS_F = 0.0000;      // [-]
  static constexpr double MARS_F = 1.0 / 169.8;  // [-]
  static constexpr double JUPITER_F = 0.06487;   // [-]
  static constexpr double SATURN_F = 0.09796;    // [-]
  static constexpr double URANUS_F = 0.02293;    // [-]
  static constexpr double NEPUTUNE_F = 0.01708;  // [-]

  // Sideral Rotation rate
  static constexpr double OMEGA_SUN = 2.903e-6;           // [rad/s]
  static constexpr double OMEGA_MERCURY = 1.244e-5;       // [rad/s]
  static constexpr double OMEGA_VENUS = -1.9521515e-7;    // [rad/s]
  static constexpr double OMEGA_EARTH = 7.2921151467e-5;  // [rad/s]
  static constexpr double OMEGA_MOON = 2.6617e-6;         // [rad/s]
  static constexpr double OMEGA_MARS = 7.0882185e-5;      // [rad/s]
  static constexpr double OMEGA_JUPITER = 1.758e-4;       // [rad/s]
  static constexpr double OMEGA_SATURN = 1.624e-4;        // [rad/s]
  static constexpr double OMEGA_URANUS = -1.036e-4;       // [rad/s]
  static constexpr double OMEGA_NEPTUNE = 1.083e-4;       // [rad/s]

  // Spherical harmonics
  static constexpr double J2_EARTH = 1.08262668e-3;
  // static constexpr double J2_MOON = 9.08901807506000e-5;
  static constexpr double J2_MOON
      = 9.094278450270e-5;  // Zonal value adjusted for permanent tide - Rigid J2
  static constexpr double C22_MOON
      = 3.470983013194e-5;  // Sectorial value adjusted for perm. tide - Rigid C22
  static constexpr double J2_MARS = 1.96045e-3;  // J2 value for Mars

  // Transformations Between GCRF and Mean Equator and Equinox at J2000
  static constexpr double FRAME_BIAS_XI0 = -16.6170e-3 * RAD_ARCSEC;   // [rad]
  static constexpr double FRAME_BIAS_ETA0 = -6.8192e-3 * RAD_ARCSEC;   // [rad]
  static constexpr double FRAME_BIAS_DALPHA0 = -14.6e-3 * RAD_ARCSEC;  // [rad]

  // Solar Radiation Pressure Constants
  static constexpr double AU = 149597970e3;      // AU [m]
  static constexpr double SOLAR_FLUX_AU = 1367;  // Mean Solar Flux at 1 AU [W/m^2]
  static constexpr double C = 299792458;         // Light speed [m/s]
  static constexpr double P_SUN
      = SOLAR_FLUX_AU / C;  // Solar radiation pressure at 1 AU [N/m^2] = 4.56e-6 N/m^2

  struct PhysicalConstants {
    double GM_SUN;
    double GM_MERCURY;
    double GM_VENUS;
    double GM_EARTH;
    double GM_MOON;
    double GM_MARS_SYSTEM;
    double GM_JUPITER_SYSTEM;
    double GM_SATURN_SYSTEM;
    double GM_URANUS_SYSTEM;
    double GM_NEPTUNE_SYSTEM;
    double GM_PLUTO_SYSTEM;
    double GM_MARS;
    double GM_JUPITER;
    double GM_SATURN;
    double GM_URANUS;
    double GM_NEPTUNE;
    double GM_CERES;
    double GM_VESTA;

    double D_EARTH_MOON;
    double D_EARTH_EMB;
    double D_MOON_EMB;
    double R_EARTH;
    double R_MOON;
    double R_SUN;
    double R_MERCURY;
    double R_VENUS;
    double R_MARS;
    double R_JUPITER;
    double R_SATURN;
    double R_URANUS;
    double R_NEPTUNE;
    double R_PLUTO;
    double WGS84_A;

    double OMEGA_EARTH_MOON;
    double OMEGA_SUN;
    double OMEGA_MERCURY;
    double OMEGA_VENUS;
    double OMEGA_EARTH;
    double OMEGA_MOON;
    double OMEGA_MARS;
    double OMEGA_JUPITER;
    double OMEGA_SATURN;
    double OMEGA_URANUS;
    double OMEGA_NEPTUNE;

    double AU;
    double C;
    double P_SUN;

    CoordinateScale coordinate_scale = CoordinateScale::TDB;
  };

  constexpr PhysicalConstants GetPhysicalConstants(const UnitSystem& units = SI_UNITS) {
    return {
        .GM_SUN = units.GravitationalParameter(GM_SUN),
        .GM_MERCURY = units.GravitationalParameter(GM_MERCURY),
        .GM_VENUS = units.GravitationalParameter(GM_VENUS),
        .GM_EARTH = units.GravitationalParameter(GM_EARTH),
        .GM_MOON = units.GravitationalParameter(GM_MOON),
        .GM_MARS_SYSTEM = units.GravitationalParameter(GM_MARS_SYSTEM),
        .GM_JUPITER_SYSTEM = units.GravitationalParameter(GM_JUPITER_SYSTEM),
        .GM_SATURN_SYSTEM = units.GravitationalParameter(GM_SATURN_SYSTEM),
        .GM_URANUS_SYSTEM = units.GravitationalParameter(GM_URANUS_SYSTEM),
        .GM_NEPTUNE_SYSTEM = units.GravitationalParameter(GM_NEPTUNE_SYSTEM),
        .GM_PLUTO_SYSTEM = units.GravitationalParameter(GM_PLUTO_SYSTEM),
        .GM_MARS = units.GravitationalParameter(GM_MARS),
        .GM_JUPITER = units.GravitationalParameter(GM_JUPITER),
        .GM_SATURN = units.GravitationalParameter(GM_SATURN),
        .GM_URANUS = units.GravitationalParameter(GM_URANUS),
        .GM_NEPTUNE = units.GravitationalParameter(GM_NEPTUNE),
        .GM_CERES = units.GravitationalParameter(GM_CERES),
        .GM_VESTA = units.GravitationalParameter(GM_VESTA),
        .D_EARTH_MOON = units.Length(D_EARTH_MOON),
        .D_EARTH_EMB = units.Length(D_EARTH_EMB),
        .D_MOON_EMB = units.Length(D_MOON_EMB),
        .R_EARTH = units.Length(R_EARTH),
        .R_MOON = units.Length(R_MOON),
        .R_SUN = units.Length(R_SUN),
        .R_MERCURY = units.Length(R_MERCURY),
        .R_VENUS = units.Length(R_VENUS),
        .R_MARS = units.Length(R_MARS),
        .R_JUPITER = units.Length(R_JUPITER),
        .R_SATURN = units.Length(R_SATURN),
        .R_URANUS = units.Length(R_URANUS),
        .R_NEPTUNE = units.Length(R_NEPTUNE),
        .R_PLUTO = units.Length(R_PLUTO),
        .WGS84_A = units.Length(WGS84_A),
        .OMEGA_EARTH_MOON = units.AngularVelocity(OMEGA_EARTH_MOON),
        .OMEGA_SUN = units.AngularVelocity(OMEGA_SUN),
        .OMEGA_MERCURY = units.AngularVelocity(OMEGA_MERCURY),
        .OMEGA_VENUS = units.AngularVelocity(OMEGA_VENUS),
        .OMEGA_EARTH = units.AngularVelocity(OMEGA_EARTH),
        .OMEGA_MOON = units.AngularVelocity(OMEGA_MOON),
        .OMEGA_MARS = units.AngularVelocity(OMEGA_MARS),
        .OMEGA_JUPITER = units.AngularVelocity(OMEGA_JUPITER),
        .OMEGA_SATURN = units.AngularVelocity(OMEGA_SATURN),
        .OMEGA_URANUS = units.AngularVelocity(OMEGA_URANUS),
        .OMEGA_NEPTUNE = units.AngularVelocity(OMEGA_NEPTUNE),
        .AU = units.Length(AU),
        .C = units.Velocity(C),
        .P_SUN = units.Pressure(P_SUN),
        .coordinate_scale = CoordinateScale::TDB,
    };
  }

  inline PhysicalConstants ScalePhysicalConstantsForCoordinateScale(PhysicalConstants constants,
                                                                    CoordinateScale from,
                                                                    CoordinateScale to) {
    double ratio = CheckedCoordinateScaleRatio(from, to);
    double inverse_ratio = 1.0 / ratio;

    constants.GM_SUN *= ratio;
    constants.GM_MERCURY *= ratio;
    constants.GM_VENUS *= ratio;
    constants.GM_EARTH *= ratio;
    constants.GM_MOON *= ratio;
    constants.GM_MARS_SYSTEM *= ratio;
    constants.GM_JUPITER_SYSTEM *= ratio;
    constants.GM_SATURN_SYSTEM *= ratio;
    constants.GM_URANUS_SYSTEM *= ratio;
    constants.GM_NEPTUNE_SYSTEM *= ratio;
    constants.GM_PLUTO_SYSTEM *= ratio;
    constants.GM_MARS *= ratio;
    constants.GM_JUPITER *= ratio;
    constants.GM_SATURN *= ratio;
    constants.GM_URANUS *= ratio;
    constants.GM_NEPTUNE *= ratio;
    constants.GM_CERES *= ratio;
    constants.GM_VESTA *= ratio;

    constants.D_EARTH_MOON *= ratio;
    constants.D_EARTH_EMB *= ratio;
    constants.D_MOON_EMB *= ratio;
    constants.R_EARTH *= ratio;
    constants.R_MOON *= ratio;
    constants.R_SUN *= ratio;
    constants.R_MERCURY *= ratio;
    constants.R_VENUS *= ratio;
    constants.R_MARS *= ratio;
    constants.R_JUPITER *= ratio;
    constants.R_SATURN *= ratio;
    constants.R_URANUS *= ratio;
    constants.R_NEPTUNE *= ratio;
    constants.R_PLUTO *= ratio;
    constants.WGS84_A *= ratio;
    constants.AU *= ratio;

    constants.OMEGA_EARTH_MOON *= inverse_ratio;
    constants.OMEGA_SUN *= inverse_ratio;
    constants.OMEGA_MERCURY *= inverse_ratio;
    constants.OMEGA_VENUS *= inverse_ratio;
    constants.OMEGA_EARTH *= inverse_ratio;
    constants.OMEGA_MOON *= inverse_ratio;
    constants.OMEGA_MARS *= inverse_ratio;
    constants.OMEGA_JUPITER *= inverse_ratio;
    constants.OMEGA_SATURN *= inverse_ratio;
    constants.OMEGA_URANUS *= inverse_ratio;
    constants.OMEGA_NEPTUNE *= inverse_ratio;

    constants.coordinate_scale = to;
    return constants;
  }

  inline PhysicalConstants GetPhysicalConstants(const UnitSystem& units, CoordinateScale scale) {
    return ScalePhysicalConstantsForCoordinateScale(GetPhysicalConstants(units),
                                                    CoordinateScale::TDB, scale);
  }

  inline PhysicalConstants GetPhysicalConstants(CoordinateScale scale) {
    return GetPhysicalConstants(SI_UNITS, scale);
  }

  // Ionosphere
  static const double IONOSPHERIC_CONSTANT = 40.3e16;  // Adjusted constant for units
  static const double L1_FREQ = 1.57542e9;             // L1 frequency [Hz]

  // File Path *******************************************************************
  static constexpr std::string_view TAI_UTC_FILENAME = "tai-utc.dat";
  static constexpr std::string_view EOP_FILENAME = "eopc04_08.62-now";
  static constexpr std::string_view IAU_SOFA_FILENAME = "IAU_SOFA.DAT";

  static constexpr std::string_view UNDEFINED = "Unknown";  // Used for unknown values

  // NAIF Intefer ID codes
  // Reference:
  // https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/FORTRAN/req/naif_ids.html
  enum class BodyId {
    SSB = 0,
    SOLAR_SYSTEM_BARYCENTER = SSB,
    MERCURY_BARYCENTER = 1,
    VENUS_BARYCENTER = 2,
    EMB = 3,
    EARTH_MOON_BARYCENTER = EMB,
    MARS_BARYCENTER = 4,
    JUPITER_BARYCENTER = 5,
    SATURN_BARYCENTER = 6,
    URANUS_BARYCENTER = 7,
    NEPTUNE_BARYCENTER = 8,
    PLUTO_BARYCENTER = 9,
    SUN = 10,
    MERCURY = 199,
    VENUS = 299,
    EARTH = 399,
    MOON = 301,
    MARS = 499,
    PHOBOS = 401,
    DEIMOS = 402,
    JUPITER = 599,
    SATURN = 699,
    URANUS = 799,
    NEPTUNE = 899,
  };

  double GetBodyRadius(BodyId body);

  enum class Time {
    UT1,     // Universal Time 1
    UTC,     // Coordinated Universal Time
    TAI,     // International Atomic Time
    TDB,     // Barycentric Dynamical Time
    TT,      // Terrestrial Time
    TCG,     // Geocentric Coordinate Time
    TCB,     // Barycentric Coordinate Time
    GPS,     // Global Positioning System Time
    JD_TT,   // Julian Date relative to TT
    JD_TDB,  // Julian Date relative to TDB
    TCL,     // Lunar Coordinate Time
    LT,      // Lunar Time
  };

}  // namespace lupnt