orbit.h

Go to the documentation of this file.

/* Conquer Space
 * Copyright (C) 2021-2025 Conquer Space
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <https://www.gnu.org/licenses/>.
 */
#pragma once
 
#include <math.h>
 
#include <ostream>
 
#include <entt/entt.hpp>
#include <glm/vec3.hpp>
 
#include "core/components/bodies.h"
#include "core/components/coordinates.h"
#include "core/components/units.h"
 
namespace cqsp::core::components::types {
typedef glm::dvec3 Vec3AU;
 
double GetOrbitingRadius(const double e, const kilometer a, const radian v);
 
struct Orbit {
    double eccentricity = 0;
 
    kilometer semi_major_axis = 0;
 
    radian inclination = 0;
 
    radian LAN = 0;
 
    radian w = 0;
 
    radian M0 = 0;
 
    double epoch = 0;
 
    radian v = 0;
 
    double GM = SunMu;
 
    // So we can prepare for moons and stuff
    entt::entity reference_body = entt::null;
 
    Orbit() = default;
 
    // TODO(EhWhoAmI): Mean anomaly is not the true anomaly, we need to convert
    // mean anomaly to true anomaly
    Orbit(kilometer semi_major_axis, double eccentricity, radian inclination, radian LAN, radian w, radian M0)
        : eccentricity(eccentricity),
          semi_major_axis(semi_major_axis),
          inclination(std::fmod(inclination, PI)),
          LAN(std::fmod(LAN, TWOPI)),
          w(w),
          M0(M0),
          v(M0),
          epoch(0),
          reference_body(entt::null) {}
 
    Orbit(kilometer semi_major_axis, double eccentricity, radian inclination, radian LAN, radian w, radian M0,
          entt::entity reference)
        : eccentricity(eccentricity),
          semi_major_axis(semi_major_axis),
          inclination(std::fmod(inclination, PI)),
          LAN(std::fmod(LAN, TWOPI)),
          w(w),
          M0(M0),
          v(M0),
          epoch(0),
          reference_body(reference) {}
 
    Orbit(const Orbit& orbit)
        : eccentricity(orbit.eccentricity),
          semi_major_axis(orbit.semi_major_axis),
          inclination(orbit.inclination),
          LAN(orbit.LAN),
          w(orbit.w),
          M0(orbit.M0),
          v(orbit.M0),
          epoch(orbit.epoch),
          GM(orbit.GM),
          reference_body(orbit.reference_body) {}
 
    double GetMtElliptic(double time) const { return normalize_radian(M0 + (time - epoch) * nu()); }
 
    double GetOrbitingRadius() const;
 
    double GetOrbitingRadius(double v) const;
 
    // Orbital period
    double T() const { return 2 * PI * sqrt(semi_major_axis * semi_major_axis * semi_major_axis / GM); }
 
    // Mean motion
    double nu() const { return sqrt(GM / abs(semi_major_axis * semi_major_axis * semi_major_axis)); }
 
    double GetApoapsis() const { return semi_major_axis * (1 + eccentricity); }
 
    double GetPeriapsis() const { return semi_major_axis * (1 - eccentricity); }
 
    double TimeToTrueAnomaly(double v2) const;
 
    double OrbitalVelocity() const;
 
    double OrbitalVelocityAtTrueAnomaly(double true_anomaly) const;
 
    std::string ToHumanString() const;
};
 
inline std::ostream& operator<<(std::ostream& outs, const Orbit& orb) {
    return outs << "(a=" << orb.semi_major_axis << ", e=" << orb.eccentricity << ", i=" << orb.inclination
                << ", LAN=" << orb.LAN << ", w=" << orb.w << ", GM=" << orb.GM << ", v=" << orb.v << ", M0=" << orb.M0
                << ", t=" << orb.epoch << ", ref=" << (uint32_t)orb.reference_body << ")";
}
 
constexpr double apoapsis = PI;
constexpr double periapsis = 0;
 
struct SetTrueAnomaly {
    double true_anomaly;
};
 
glm::dvec3 ConvertOrbParams(const double LAN, const double i, const double w, const glm::dvec3& vec);
 
glm::dvec3 ConvertToOrbitalVector(const double LAN, const double i, const double w, const double v,
                                  const glm::dvec3& vec);
 
glm::dvec3 MatrixConvertOrbParams(const double LAN, const double i, const double w, const glm::dvec3& vec);
 
void UpdateOrbit(Orbit& orb, const second& time);
 
double GetTrueAnomaly(const Orbit& orb, const second& epoch);
 
glm::dvec3 CalculateVelocity(const double& E, const kilometer& r, const double& GM, const kilometer& a,
                             const double& e);
 
glm::dvec3 CalculateVelocityElliptic(const double& E, const kilometer& r, const double& GM, const kilometer& a,
                                     const double& e);
 
glm::dvec3 CalculateVelocityHyperbolic(const double& E, const double& r, const double& GM, const double& a,
                                       const double& e);
 
double GetCircularOrbitingVelocity(const double& GM, const double& radius);
 
double GetTrueAnomalyToAsymptope(const Orbit& orbit);
 
Orbit Vec3ToOrbit(const glm::dvec3& position, const glm::dvec3& velocity, const double& GM, const double& time);
 
glm::dvec3 OrbitToVec3(const double& a, const double& e, const radian& i, const radian& LAN, const radian& w,
                       const radian& v);
 
double OrbitVelocity(const double v, const double e, const double a, const double GM);
double OrbitVelocityAtR(const double GM, const double a, const double r);
double AvgOrbitalVelocity(const Orbit& orb);
 
glm::dvec3 OrbitVelocityToVec3(const Orbit& orb, double v);
glm::dvec3 OrbitVelocityToVec3(const Orbit& orb);
 
double SolveKeplerElliptic(const double& mean_anomaly, const double& ecc, const int steps = 200);
 
double SolveKeplerHyperbolic(const double& mean_anomaly, const double& ecc, const int steps = 200);
 
double EccentricAnomalyToTrueAnomaly(const double& ecc, const double& E);
 
double HyperbolicAnomalyToTrueAnomaly(const double& ecc, const double& H);
 
double GetMtElliptic(const double& M0, const double& nu, const double& time, const double& epoch);
 
double GetMtHyperbolic(const double& M0, const double& nu, const double& time, const double& epoch);
 
radian TrueAnomalyElliptic(const Orbit& orbit, const second& time);
radian TrueAnomalyElliptic(const Orbit& orbit, const second& time, double& E_out);
 
radian TrueAnomalyHyperbolic(const Orbit& orbit, const second& time);
 
radian EccentricAnomaly(double v, double e);
radian HyperbolicAnomaly(double v, double e);
inline Vec3AU toVec3AU(const Orbit& orb, radian theta) {
    glm::dvec3 vec = OrbitToVec3(orb.semi_major_axis, orb.eccentricity, orb.inclination, orb.LAN, orb.w, theta);
    return vec / KmInAu;
}
 
Orbit ApplyImpulse(const Orbit& orbit, const glm::dvec3& impulse, double time);
 
inline glm::dvec3 toVec3(const Orbit& orb, radian theta) {
    return OrbitToVec3(orb.semi_major_axis, orb.eccentricity, orb.inclination, orb.LAN, orb.w, theta);
}
 
glm::dvec3 OrbitTimeToVec3(const Orbit& orb, const second time = 0);
glm::dvec3 OrbitTimeToVelocityVec3(const Orbit& orb, const second time = 0);
 
inline glm::dvec3 toVec3(const Orbit& orb) { return toVec3(orb, orb.v); }
inline Vec3AU toVec3AU(const Orbit& orb) { return toVec3AU(orb, orb.v); }
 
inline void UpdatePos(Kinematics& kin, const Orbit& orb) {
    // Calculate time
    kin.position = toVec3AU(orb);
}
 
double CalculateTransferTime(const Orbit& orb1, const Orbit& orb2);
double CalculateTransferAngle(const Orbit& start_orbit, const Orbit& end_orbit);
double CalculatePhaseAngle(const Orbit& start_orbit, const Orbit& end_orbit, double epoch);
 
// https://orbital-mechanics.space/the-orbit-equation/hyperbolic-trajectories.html
// True anomaly of the asymptope for a hyperbolic orbit
// The hyperbolic asymtope is defined by -GetHyperbolicAsymtope < v <GetHyperbolicAsymtope
double GetHyperbolicAsymptopeAnomaly(double eccentricity);
 
double FlightPathAngle(double eccentricity, double v);
 
glm::dvec3 GetOrbitNormal(const Orbit& orbit);
glm::dvec3 GetRadialVector(const Orbit& orbit);
glm::dvec3 GetRadialVector(const Orbit& orbit, double true_anomaly);
glm::dvec3 InvertOrbitalVector(const double LAN, const double i, const double w, const double v, const glm::dvec3& vec);
double AngleWith(const Orbit& orbit, const Orbit& second_orbit);
double TrueAnomalyFromVector(const Orbit& orbit, const glm::dvec3& vec);
 
double AscendingTrueAnomaly(const Orbit& start, const Orbit& dest);
double GetEccentricity(double apoapsis, double periapsis);
glm::dvec3 GetEccentricityVector(const glm::dvec3& h, const glm::dvec3& r, const glm::dvec3& v, const double GM);
glm::dvec3 GetEccentricityVector(const glm::dvec3& r, const glm::dvec3& v, const double GM);
}  // namespace cqsp::core::components::types