diff --git a/src/game/astronomy/celestial-mechanics.cpp b/src/game/astronomy/celestial-mechanics.cpp index 74fd563..a18a7fc 100644 --- a/src/game/astronomy/celestial-mechanics.cpp +++ b/src/game/astronomy/celestial-mechanics.cpp @@ -129,72 +129,53 @@ double3x3 approx_moon_ecliptic_rotation(double jd) return rz0 * rx * rz1; } -double3 solve_kepler(const kepler_orbit& orbit, double t) +double solve_kepler(double ec, double ma, double tolerance, std::size_t iterations) { - // Orbital period (Kepler's third law) - double pr = std::sqrt(orbit.a * orbit.a * orbit.a); + // Approximate eccentric anomaly, E + double ea0 = ma + ec * std::sin(ma) * (1.0 + ec * std::cos(ma)); - // Mean anomaly - double epoch = 0.0; - double ma = (math::two_pi * (t - epoch)) / pr; - - // Semi-minor axis - const double b = std::sqrt(1.0 - orbit.ec * orbit.ec); - - // Trig of orbital elements (can be precalculated?) - const double cos_ma = std::cos(ma); - const double sin_ma = std::sin(ma); - const double cos_om = std::cos(orbit.om); - const double sin_om = std::sin(orbit.om); - const double cos_i = std::cos(orbit.i); - const double sin_i = std::sin(orbit.i); - - // Eccentric anomaly - double ea = ma + orbit.ec * sin_ma * (1.0 + orbit.ec * cos_ma); - - // Find radial distance (r) and true anomaly (v) - double x = orbit.a * (std::cos(ea) - orbit.ec); - double y = b * std::sin(ea); - double r = std::sqrt(x * x + y * y); - double v = std::atan2(y, x); - - // Convert (r, v) to ecliptic rectangular coordinates - double cos_wv = std::cos(orbit.w + v); - double sin_wv = std::sin(orbit.w + v); - return double3 + // Iteratively converge ea0 and ea1 + for (std::size_t i = 0; i < iterations; ++i) { - r * (cos_om * cos_wv - sin_om * sin_wv * cos_i), - r * (sin_om * cos_wv + cos_om * sin_wv * cos_i), - r * sin_wv * sin_i - }; + double ea1 = ea0; + ea0 = ea1 - (ea1 - ec * std::sin(ea1) - ma) / (1.0 - ec * std::cos(ea1)); + if (std::abs(ea1 - ea0) < tolerance) + break; + } + + return ea0; } -double3 solve_kepler(double a, double ec, double w, double ma, double i, double om) +orbital_state orbital_elements_to_state(const orbital_elements& elements, double ke_tolerance, std::size_t ke_iterations) { - // Semi-minor axis - double b = std::sqrt(1.0 - ec * ec); + orbital_state state; + + // Calculate semi-minor axis, b + double b = std::sqrt(1.0 - elements.ec * elements.ec); - // Eccentric anomaly - double ea = ma + ec * std::sin(ma) * (1.0 + ec * std::cos(ma)); + // Solve Kepler's equation for eccentric anomaly, E + double ea = solve_kepler(elements.ec, elements.ma, ke_tolerance, ke_iterations); // Radial distance (r) and true anomaly (v) - double x = a * (std::cos(ea) - ec); + double x = elements.a * (std::cos(ea) - elements.ec); double y = b * std::sin(ea); double r = std::sqrt(x * x + y * y); double v = std::atan2(y, x); // Convert (r, v) to ecliptic rectangular coordinates - double cos_om = std::cos(om); - double sin_om = std::sin(om); - double cos_i = std::cos(i); - double cos_wv = std::cos(w + v); - double sin_wv = std::sin(w + v); - return double3 + double cos_om = std::cos(elements.om); + double sin_om = std::sin(elements.om); + double cos_i = std::cos(elements.i); + double cos_wv = std::cos(elements.w + v); + double sin_wv = std::sin(elements.w + v); + state.r = { r * (cos_om * cos_wv - sin_om * sin_wv * cos_i), r * (sin_om * cos_wv + cos_om * sin_wv * cos_i), - r * sin_wv * std::sin(i) + r * sin_wv * std::sin(elements.i) }; + + return state; } } // namespace ast diff --git a/src/game/astronomy/celestial-mechanics.hpp b/src/game/astronomy/celestial-mechanics.hpp index 60443a8..9ef4617 100644 --- a/src/game/astronomy/celestial-mechanics.hpp +++ b/src/game/astronomy/celestial-mechanics.hpp @@ -25,6 +25,28 @@ namespace ast { +/** + * Contains six orbital elements which describe a Keplerian orbit. + */ +struct orbital_elements +{ + double ec; ///< Eccentricity, e + double a; ///< Semi-major axis, a + double i; ///< Inclination, i (radians) + double om; ///< Longitude of the ascending node, OMEGA (radians) + double w; ///< Argument of periapsis, w (radians) + double ma; ///< Mean anomaly, M (radians) +}; + +/** + * Orbital state vectors. + */ +struct orbital_state +{ + double3 r; ///< Cartesian position, r. + double3 v; ///< Cartesian velocity, v. +}; + /** * Approximates the obliquity of the ecliptic. * @@ -57,28 +79,26 @@ double3 approx_moon_ecliptic(double jd); */ double3x3 approx_moon_ecliptic_rotation(double jd); -struct kepler_orbit -{ - double a; ///< Semi-major axis, a - double ec; ///< Eccentricity, e - double w; ///< Argument of periapsis, w (radians) - double ma; ///< Mean anomaly, M (radians) - double i; ///< Inclination, i (radians) - double om; ///< Longitude of the ascending node, OMEGA (radians) -}; - -double3 solve_kepler(const kepler_orbit& orbit, double t); - /** + * Iteratively solves Kepler's equation for eccentric anomaly, E. * - * @param a Semi-major axis, a. * @param ec Eccentricity, e. - * @param w Argument of periapsis, w (radians). * @param ma Mean anomaly, M (radians). - * @param i Inclination, i (radians). - * @param om Longitude of the ascending node, OMEGA (radians). + * @param tolerance Tolerance of solution. + * @param iterations Maximum number of iterations. + * @return Eccentric anomaly. + */ +double solve_kepler(double ec, double ma, double tolerance, std::size_t iterations); + +/** + * Calculates orbital state vectors from Keplerian orbital elements. + * + * @param elements Orbital elements. + * @param ke_tolerance Kepler's equation tolerance. + * @param ke_iterations Kepler's equation iterations. + * @return Orbital state. */ -double3 solve_kepler(double a, double ec, double w, double ma, double i, double om); +orbital_state orbital_elements_to_state(const orbital_elements& elements, double ke_tolerance, std::size_t ke_iterations); } // namespace ast diff --git a/src/game/systems/solar-system.cpp b/src/game/systems/solar-system.cpp index 59e13a3..691f862 100644 --- a/src/game/systems/solar-system.cpp +++ b/src/game/systems/solar-system.cpp @@ -23,6 +23,7 @@ #include "game/astronomy/celestial-time.hpp" #include "game/components/orbit-component.hpp" #include "game/components/transform-component.hpp" +#include using namespace ecs; @@ -42,10 +43,39 @@ void solar_system::update(double t, double dt) // Add scaled timestep to Julian date set_julian_date(julian_date + (dt * time_scale) / seconds_per_day); + + const double ke_tolerance = 1e-6; + const std::size_t ke_iterations = 10; + + ast::orbital_elements sun_elements; + sun_elements.a = 1.0; + sun_elements.ec = 0.016709; + sun_elements.w = math::radians(282.9404); + sun_elements.ma = math::radians(356.0470); + sun_elements.i = math::radians(0.0); + sun_elements.om = math::radians(0.0); + + const double j2k_day = julian_date - 2451545.0; + const double j2k_century = (julian_date - 2451545.0) / 36525.0; + + sun_elements.ec += math::radians(-1.151e-9) * j2k_day; + sun_elements.w += math::radians(4.70935e-5) * j2k_day; + sun_elements.ma += math::radians(0.9856002585) * j2k_day; + + + ast::orbital_state sun_ecliptic = ast::orbital_elements_to_state(sun_elements, ke_tolerance, ke_iterations); + double3 sun_horizontal = ecliptic_to_horizontal * sun_ecliptic.r; + sun_horizontal.z -= 4.25875e-5; + double3 sun_spherical = ast::rectangular_to_spherical(sun_horizontal); + double2 sun_az_el = {sun_spherical.z - math::pi, sun_spherical.y}; + + std::cout << "new azel: " << math::degrees(sun_az_el.x) << ", " << math::degrees(sun_az_el.y) << std::endl; + // Update horizontal (topocentric) positions of orbiting bodies registry.view().each( [&](auto entity, auto& orbit, auto& transform) { + /* double a = orbit.a; double ec = orbit.ec; double w = orbit.w; @@ -65,6 +95,7 @@ void solar_system::update(double t, double dt) transform.local.translation = math::type_cast(translation); transform.local.rotation = math::type_cast(math::quaternion_cast(rotation)); + */ }); }