Add iterative solution to Kepler's equation, improve orbital position calculation

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<double> * (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

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
 

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 <iostream>
 
 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<double>, 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<orbit_component, transform_component>().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<float>(translation);
 		transform.local.rotation = math::type_cast<float>(math::quaternion_cast(rotation));
+		*/
 	});
 }