Improve calculations of orbits and celestial body orientations.

src/animation/ease.hpp

@@ -339,7 +339,6 @@ T ease::in_out_elastic(const T& x, const T& y, S a)
 	if (a < S(0.5))
 	{
 		return math::lerp(x, y, std::pow(S(2), S(20) * a - S(11)) * std::sin(S(15.5334) - S(27.5293) * a));
-		
 	}
 	else
 	{

src/astro/illuminance.hpp

@@ -22,8 +22,7 @@
 
 namespace astro
 {
-	
-	
+
 /**
  * Converts apparent (visual) magnitude to a brightness factor relative to a 0th magnitude star.
  *

src/entity/components/celestial-body.hpp

@@ -26,17 +26,23 @@ namespace component {
 /// A simple celestial body.
 struct celestial_body
 {
-	/// Body radius, in meters.
+	/// Mean radius of the body, in meters.
 	double radius;
 	
-	/// Angle between the body's rotational axis and its orbital axis, in radians.
-	double axial_tilt;
+	/// Mass of the body, in kilograms.
+	double mass;
 	
-	/// Angle of rotation about the body's rotational axis at epoch, in radians.
-	double axial_rotation;
+	/// Right ascension of the body's north pole at epoch, in radians.
+	double pole_ra;
 	
-	/// Angular frequency, in radians per day.
-	double angular_frequency;
+	/// Declination of the body's north pole at epoch, in radians.
+	double pole_dec;
+	
+	/// Location of the prime meridian at epoch, as a rotation about the north pole, in radians.
+	double prime_meridian;
+	
+	/// Sidereal rotation period, in rotations per day.
+	double rotation_period;
 };
 
 } // namespace component

src/entity/components/orbit.hpp

@@ -20,16 +20,36 @@
 #ifndef ANTKEEPER_ENTITY_COMPONENT_ORBIT_HPP
 #define ANTKEEPER_ENTITY_COMPONENT_ORBIT_HPP
 
+#include "entity/id.hpp"
 #include "physics/orbit/elements.hpp"
 #include "physics/orbit/state.hpp"
+#include "math/se3.hpp"
 
 namespace entity {
 namespace component {
 
 struct orbit
 {
+	/// Parent body.
+	entity::id parent;
+	
+	/// Keplerian orbital elements at epoch.
 	physics::orbit::elements<double> elements;
-	physics::orbit::state<double> state;
+	
+	/// Orbital semi-minor axis (b) at epoch.
+	double semiminor_axis;
+	
+	/// Orbital mean motion (n) at epoch.
+	double mean_motion;
+	
+	/// Transformation from the PQW frame to the BCI frame of the parent body.
+	math::transformation::se3<double> pqw_to_bci;
+	
+	/// Orbital Cartesian position of the body in the BCI frame of the parent body.
+	math::vector3<double> bci_position;
+	
+	/// Orbital Cartesian position of the body in the ICRF frame.
+	math::vector3<double> icrf_position;
 };
 
 } // namespace component

src/entity/systems/astronomy.cpp

@@ -22,6 +22,7 @@
 #include "entity/components/blackbody.hpp"
 #include "entity/components/transform.hpp"
 #include "geom/intersection.hpp"
+#include "geom/cartesian.hpp"
 #include "color/color.hpp"
 #include "physics/orbit/orbit.hpp"
 #include "physics/time/ut1.hpp"
@@ -57,22 +58,16 @@ astronomy::astronomy(entity::registry& registry):
 	reference_entity(entt::null),
 	observer_location{0, 0, 0},
 	sun_light(nullptr),
+	sky_light(nullptr),
 	sky_pass(nullptr)
 {
-	// Construct reference frame which transforms coordinates from SEZ to EZS
-	sez_to_ezs = physics::frame<double>
+	// Construct transformation which transforms coordinates from ENU to EUS
+	enu_to_eus = math::transformation::se3<double>
 	{
 		{0, 0, 0},
-		math::normalize
-		(
-			math::quaternion<double>::rotate_x(-math::half_pi<double>) *
-				math::quaternion<double>::rotate_z(-math::half_pi<double>)
-		)
+		math::quaternion<double>::rotate_x(-math::half_pi<double>)
 	};
 	
-	// Construct reference frame which transforms coordinates from EZS to SEZ
-	ezs_to_sez = sez_to_ezs.inverse();
-	
 	registry.on_construct<entity::component::celestial_body>().connect<&astronomy::on_celestial_body_construct>(this);
 	registry.on_replace<entity::component::celestial_body>().connect<&astronomy::on_celestial_body_replace>(this);
 }
@@ -92,47 +87,91 @@ void astronomy::update(double t, double dt)
 	const entity::component::orbit& reference_orbit = registry.get<entity::component::orbit>(reference_entity);
 	const entity::component::celestial_body& reference_body = registry.get<entity::component::celestial_body>(reference_entity);
 	
-	// Determine axial rotation at current time
-	const double reference_axial_rotation = reference_body.axial_rotation + reference_body.angular_frequency * universal_time;
+	math::transformation::se3<double> icrf_to_bci{{-reference_orbit.icrf_position}, math::identity_quaternion<double>};
 	
-	// Construct reference frame which transforms coordinates from inertial space to reference body BCBF space
-	inertial_to_bcbf = physics::orbit::inertial::to_bcbf
+	// Construct transformation from the ICRF frame to the reference body BCBF frame	
+	icrf_to_bcbf = physics::orbit::frame::bci::to_bcbf
 	(
-		reference_orbit.state.r,
-		reference_orbit.elements.i,
-		reference_body.axial_tilt,
-		reference_axial_rotation
+		reference_body.pole_ra,
+		reference_body.pole_dec,
+		reference_body.prime_meridian + (math::two_pi<double> / reference_body.rotation_period) * universal_time
 	);
+	icrf_to_bcbf.t = icrf_to_bcbf.r * -reference_orbit.icrf_position;
 	
-	// Construct reference frame which transforms coordinates from inertial space to reference body topocentric space
-	inertial_to_topocentric = inertial_to_bcbf * bcbf_to_topocentric;
+	icrf_to_enu = icrf_to_bcbf * bcbf_to_enu;
+	icrf_to_eus = icrf_to_enu * enu_to_eus;
 	
 	// Set the transform component translations of orbiting bodies to their topocentric positions
 	registry.view<component::celestial_body, component::orbit, component::transform>().each(
-	[&](entity::id entity_id, const auto& celestial_body, const auto& orbit, auto& transform)
+	[&](entity::id entity_id, const auto& body, const auto& orbit, auto& transform)
 	{
-		// Transform Cartesian position vector (r) from inertial space to topocentric space
-		const math::vector3<double> r_topocentric = inertial_to_topocentric * orbit.state.r;
+		// Skip reference entity
+		if (entity_id == this->reference_entity)
+			return;
+		
+		// Transform orbital Cartesian position (r) from the ICRF frame to the EUS frame
+		const double3 r_eus = icrf_to_eus * orbit.icrf_position;
+		
+		// Determine body orientation in the ICRF frame
+		math::quaternion<double> rotation_icrf = physics::orbit::frame::bcbf::to_bci
+		(
+			body.pole_ra,
+			body.pole_dec,
+			body.prime_meridian + (math::two_pi<double> / body.rotation_period) * this->universal_time
+		).r;
+		
+		// Transform body orientation from the ICRF frame to the EUS frame.
+		math::quaternion<double> rotation_eus = math::normalize(icrf_to_eus.r * rotation_icrf);
 		
 		// Update local transform
-		transform.local.translation = math::type_cast<float>(r_topocentric);
+		if (orbit.parent != entt::null)
+		{
+			transform.local.translation = math::normalize(math::type_cast<float>(r_eus)) * 1000.0f;
+			transform.local.rotation = math::type_cast<float>(rotation_eus);
+			transform.local.scale = {50.0f, 50.0f, 50.0f};
+		}
+		
+		/*
+		if (orbit.parent == entt::null)
+		{
+			// RA-DEC
+			const double3 r_bci = icrf_to_bci * orbit.icrf_position;
+			double3 r_bci_spherical = physics::orbit::frame::bci::spherical(r_bci);
+			if (r_bci_spherical.z < 0.0)
+				r_bci_spherical.z += math::two_pi<double>;
+			const double dec = math::degrees(r_bci_spherical.y);
+			const double ra = math::degrees(r_bci_spherical.z);
+			
+			// AZ-EL
+			const double3 r_enu = icrf_to_enu * orbit.icrf_position;
+			double3 r_enu_spherical = physics::orbit::frame::enu::spherical(r_enu);
+			if (r_enu_spherical.z < 0.0)
+				r_enu_spherical.z += math::two_pi<double>;
+			const double el = math::degrees(r_enu_spherical.y);
+			const double az = math::degrees(r_enu_spherical.z);
+			
+			std::cout << "t: " << this->universal_time << "; ra: " << ra << "; dec: " << dec << std::endl;
+			std::cout << "t: " << this->universal_time << "; az: " << az << "; el: " << el << std::endl;
+		}
+		*/
 	});
 	
 	// Update blackbody lighting
 	registry.view<component::celestial_body, component::orbit, component::blackbody>().each(
-	[&](entity::id entity_id, const auto& celestial_body, const auto& orbit, const auto& blackbody)
+	[&](entity::id entity_id, const auto& body, const auto& orbit, const auto& blackbody)
 	{
 		// Calculate blackbody inertial basis
-		double3 blackbody_forward_inertial = math::normalize(reference_orbit.state.r - orbit.state.r);
-		double3 blackbody_up_inertial = {0, 0, 1};
+		//double3 blackbody_forward_icrf = math::normalize(reference_orbit.icrf_position - orbit.icrf_position);
+		double3 blackbody_up_icrf = {0, 0, 1};
 		
-		// Transform blackbody inertial position and basis into topocentric space
-		double3 blackbody_position_topocentric = inertial_to_topocentric * orbit.state.r;
-		double3 blackbody_forward_topocentric = inertial_to_topocentric.rotation * blackbody_forward_inertial;
-		double3 blackbody_up_topocentric = inertial_to_topocentric.rotation * blackbody_up_inertial;
+		// Transform blackbody ICRF position and basis into EUS frame
+		double3 blackbody_position_eus = icrf_to_eus * orbit.icrf_position;
+		double3 blackbody_position_enu = icrf_to_enu * orbit.icrf_position;
+		double3 blackbody_forward_eus = math::normalize(-blackbody_position_eus);
+		double3 blackbody_up_eus = icrf_to_eus.r * blackbody_up_icrf;
 		
 		// Calculate distance from observer to blackbody
-		double blackbody_distance = math::length(blackbody_position_topocentric) - celestial_body.radius;
+		double blackbody_distance = math::length(blackbody_position_eus) - body.radius;
 		
 		// Calculate blackbody distance attenuation
 		double distance_attenuation = 1.0 / (blackbody_distance * blackbody_distance);
@@ -148,7 +187,7 @@ void astronomy::update(double t, double dt)
 			// Altitude of observer in meters	
 			geom::ray<double> sample_ray;
 			sample_ray.origin = {0, reference_body.radius + observer_location[0], 0};
-			sample_ray.direction = math::normalize(blackbody_position_topocentric);
+			sample_ray.direction = math::normalize(blackbody_position_eus);
 			
 			geom::sphere<double> exosphere;
 			exosphere.center = {0, 0, 0};
@@ -172,49 +211,58 @@ void astronomy::update(double t, double dt)
 		if (sun_light != nullptr)
 		{
 			// Update blackbody light transform
-			sun_light->set_translation(math::normalize(math::type_cast<float>(blackbody_position_topocentric)));
+			sun_light->set_translation(math::normalize(math::type_cast<float>(blackbody_position_eus)));
 			sun_light->set_rotation
 			(
 				math::look_rotation
 				(
-					math::type_cast<float>(blackbody_forward_topocentric),
-					math::type_cast<float>(blackbody_up_topocentric)
+					math::type_cast<float>(blackbody_forward_eus),
+					math::type_cast<float>(blackbody_up_eus)
 				)
 			);
 			
+			// Sun illuminance at the outer atmosphere
+			float3 sun_illuminance_outer = math::type_cast<float>(blackbody.luminous_intensity * distance_attenuation);
 			
-			// Sun color at the outer atmosphere
-			float3 sun_color_outer = math::type_cast<float>(blackbody.luminous_intensity * distance_attenuation);
-			
-			// Sun color at sea level
-			float3 sun_color_inner = math::type_cast<float>(blackbody.luminous_intensity * distance_attenuation * atmospheric_transmittance);
+			// Sun illuminance at sea level
+			float3 sun_illuminance_inner = math::type_cast<float>(blackbody.luminous_intensity * distance_attenuation * atmospheric_transmittance);
 			
 			// Update blackbody light color and intensity
-			sun_light->set_color(sun_color_inner);	
+			sun_light->set_color(sun_illuminance_inner);	
 			sun_light->set_intensity(1.0f);
 			
 			// Upload blackbody params to sky pass
 			if (this->sky_pass)
 			{
-				this->sky_pass->set_sun_position(math::type_cast<float>(blackbody_position_topocentric));
-				this->sky_pass->set_sun_color(sun_color_outer, sun_color_inner);
+				this->sky_pass->set_sun_position(math::type_cast<float>(blackbody_position_eus));
+				this->sky_pass->set_sun_illuminance(sun_illuminance_outer, sun_illuminance_inner);
 				
-				double blackbody_angular_radius = std::asin((celestial_body.radius * 2.0) / (blackbody_distance * 2.0));
+				double blackbody_angular_radius = std::asin((body.radius * 2.0) / (blackbody_distance * 2.0));
 				this->sky_pass->set_sun_angular_radius(static_cast<float>(blackbody_angular_radius));
 			}
 		}
+		
+		if (sky_light != nullptr)
+		{
+			double3 blackbody_position_enu_spherical = physics::orbit::frame::enu::spherical(icrf_to_enu * orbit.icrf_position);
+			
+			double illuminance = 25000.0 * std::max<double>(0.0, std::sin(blackbody_position_enu_spherical.y));
+			
+			sky_light->set_color({1.0f, 1.0f, 1.0f});
+			sky_light->set_intensity(static_cast<float>(illuminance));
+		}
 	});
 	
 	// Update sky pass topocentric frame
 	if (sky_pass != nullptr)
 	{
 		// Upload topocentric frame to sky pass
-		sky_pass->set_topocentric_frame
+		sky_pass->set_icrf_to_eus
 		(
-			physics::frame<float>
+			math::transformation::se3<float>
 			{
-				math::type_cast<float>(inertial_to_topocentric.translation),
-				math::type_cast<float>(inertial_to_topocentric.rotation)
+				math::type_cast<float>(icrf_to_eus.t),
+				math::type_cast<float>(icrf_to_eus.r)
 			}
 		);
 		
@@ -247,13 +295,13 @@ void astronomy::set_time_scale(double scale)
 void astronomy::set_reference_body(entity::id entity_id)
 {
 	reference_entity = entity_id;
-	update_bcbf_to_topocentric();
+	update_bcbf_to_enu();
 }
 
 void astronomy::set_observer_location(const double3& location)
 {
 	observer_location = location;
-	update_bcbf_to_topocentric();
+	update_bcbf_to_enu();
 }
 
 void astronomy::set_sun_light(scene::directional_light* light)
@@ -261,6 +309,11 @@ void astronomy::set_sun_light(scene::directional_light* light)
 	sun_light = light;
 }
 
+void astronomy::set_sky_light(scene::ambient_light* light)
+{
+	sky_light = light;
+}
+
 void astronomy::set_sky_pass(::render::sky_pass* pass)
 {
 	this->sky_pass = pass;
@@ -269,16 +322,16 @@ void astronomy::set_sky_pass(::render::sky_pass* pass)
 void astronomy::on_celestial_body_construct(entity::registry& registry, entity::id entity_id, entity::component::celestial_body& celestial_body)
 {
 	if (entity_id == reference_entity)
-		update_bcbf_to_topocentric();
+		update_bcbf_to_enu();
 }
 
 void astronomy::on_celestial_body_replace(entity::registry& registry, entity::id entity_id, entity::component::celestial_body& celestial_body)
 {
 	if (entity_id == reference_entity)
-		update_bcbf_to_topocentric();
+		update_bcbf_to_enu();
 }
 
-void astronomy::update_bcbf_to_topocentric()
+void astronomy::update_bcbf_to_enu()
 {
 	double radial_distance = observer_location[0];
 	
@@ -288,13 +341,13 @@ void astronomy::update_bcbf_to_topocentric()
 			radial_distance += registry.get<entity::component::celestial_body>(reference_entity).radius;
 	}
 	
-	// Construct reference frame which transforms coordinates from BCBF space to topocentric space
-	bcbf_to_topocentric = physics::orbit::bcbf::to_topocentric
+	// Construct reference frame which transforms coordinates from a BCBF frame to a horizontal frame
+	bcbf_to_enu = physics::orbit::frame::bcbf::to_enu
 	(
 		radial_distance,
 		observer_location[1],
 		observer_location[2]
-	) * sez_to_ezs;
+	);
 }
 
 } // namespace system

src/entity/systems/astronomy.hpp

@@ -23,8 +23,9 @@
 #include "entity/systems/updatable.hpp"
 #include "entity/id.hpp"
 #include "scene/directional-light.hpp"
+#include "scene/ambient-light.hpp"
 #include "utility/fundamental-types.hpp"
-#include "physics/frame.hpp"
+#include "math/se3.hpp"
 #include "render/passes/sky-pass.hpp"
 #include "entity/components/atmosphere.hpp"
 #include "entity/components/celestial-body.hpp"
@@ -79,6 +80,7 @@ public:
 	void set_observer_location(const double3& location);
 	
 	void set_sun_light(scene::directional_light* light);
+	void set_sky_light(scene::ambient_light* light);
 	
 	void set_sky_pass(::render::sky_pass* pass);
 	
@@ -86,7 +88,7 @@ private:
 	void on_celestial_body_construct(entity::registry& registry, entity::id entity_id, entity::component::celestial_body& celestial_body);
 	void on_celestial_body_replace(entity::registry& registry, entity::id entity_id, entity::component::celestial_body& celestial_body);
 	
-	void update_bcbf_to_topocentric();
+	void update_bcbf_to_enu();
 
 	double universal_time;
 	double time_scale;
@@ -95,13 +97,14 @@ private:
 	
 	double3 observer_location;
 	
-	physics::frame<double> inertial_to_bcbf;
-	physics::frame<double> bcbf_to_topocentric;
-	physics::frame<double> inertial_to_topocentric;
-	physics::frame<double> sez_to_ezs;
-	physics::frame<double> ezs_to_sez;
+	math::transformation::se3<double> icrf_to_bcbf;
+	math::transformation::se3<double> bcbf_to_enu;
+	math::transformation::se3<double> icrf_to_enu;
+	math::transformation::se3<double> enu_to_eus;
+	math::transformation::se3<double> icrf_to_eus;
 	
 	scene::directional_light* sun_light;
+	scene::ambient_light* sky_light;
 	::render::sky_pass* sky_pass;
 };
 

src/entity/systems/orbit.cpp

@@ -18,9 +18,8 @@
  */
 
 #include "entity/systems/orbit.hpp"
-#include "entity/components/orbit.hpp"
-#include "entity/id.hpp"
 #include "physics/orbit/orbit.hpp"
+#include <iostream>
 
 namespace entity {
 namespace system {
@@ -31,51 +30,46 @@ orbit::orbit(entity::registry& registry):
 	time_scale(1.0),
 	ke_iterations(10),
 	ke_tolerance(1e-6)
-{}
+{
+	registry.on_construct<entity::component::orbit>().connect<&orbit::on_orbit_construct>(this);
+	registry.on_replace<entity::component::orbit>().connect<&orbit::on_orbit_replace>(this);
+}
 
 void orbit::update(double t, double dt)
 {
 	// Add scaled timestep to current time
 	set_universal_time(universal_time + dt * time_scale);
 	
-	// Update the orbital state of orbiting bodies
+	// Propagate orbits
 	registry.view<component::orbit>().each(
 	[&](entity::id entity_id, auto& orbit)
 	{
-		// Calculate semi-minor axis (b)
-		const double b = physics::orbit::derive_semiminor_axis(orbit.elements.a, orbit.elements.e);
+		// Determine mean anomaly at current time
+		const double ma = orbit.elements.ma + orbit.mean_motion * this->universal_time;
 		
 		// Solve Kepler's equation for eccentric anomaly (E)
-		const double ea = physics::orbit::kepler_ea(orbit.elements.e, orbit.elements.ta, ke_iterations, ke_tolerance);
-		
-		// Calculate radial distance and true anomaly (nu)
-		const double xv = orbit.elements.a * (std::cos(ea) - orbit.elements.e);
-		const double yv = b * std::sin(ea);
-		const double distance = std::sqrt(xv * xv + yv * yv);
-		const double ta = std::atan2(yv, xv);
-		
-		// Calculate Cartesian position (r) in perifocal space
-		const math::vector3<double> r_perifocal = math::quaternion<double>::rotate_z(ta) * math::vector3<double>{distance, 0, 0};
+		const double ea = physics::orbit::anomaly::mean_to_eccentric(orbit.elements.ec, ma, ke_iterations, ke_tolerance);
 		
-		/// @TODO Calculate Cartesian velocity (v) in perifocal space
-		//const math::vector3<double> v_perifocal = ...
+		// Calculate Cartesian orbital position in the PQW frame
+		math::vector3<double> pqw_position = physics::orbit::frame::pqw::cartesian(orbit.elements.ec, orbit.elements.a, ea, orbit.semiminor_axis);
 		
-		// Construct perifocal to inertial reference frame
-		const physics::frame<double> perifocal_to_inertial = physics::orbit::inertial::to_perifocal
-		(
-			{0, 0, 0},
-			orbit.elements.raan,
-			orbit.elements.i,
-			orbit.elements.w
-		).inverse();
-		
-		// Transform orbital state vectors from perifocal space to the parent inertial space
-		const math::vector3<double> r_inertial = perifocal_to_inertial.transform(r_perifocal);
-		//const math::vector3<double> v_inertial = perifocal_frame.transform(v_perifocal);
+		// Transform orbital position from PQW frame to BCI frame
+		orbit.bci_position = orbit.pqw_to_bci.transform(pqw_position);
+	});
+	
+	// Update orbital positions in the ICRF frame
+	registry.view<component::orbit>().each(
+	[&](entity::id entity_id, auto& orbit)
+	{
+		orbit.icrf_position = orbit.bci_position;
 		
-		// Update orbital state of component
-		orbit.state.r = r_inertial;
-		//orbit.state.v = v_inertial;
+		entity::id parent = orbit.parent;
+		while (parent != entt::null)
+		{
+			const component::orbit& parent_orbit = registry.get<component::orbit>(parent);
+			orbit.icrf_position += parent_orbit.bci_position;
+			parent = parent_orbit.parent;
+		}
 	});
 }
 
@@ -89,5 +83,27 @@ void orbit::set_time_scale(double scale)
 	time_scale = scale;
 }
 
+void orbit::on_orbit_construct(entity::registry& registry, entity::id entity_id, entity::component::orbit& component)
+{
+	component.semiminor_axis = physics::orbit::semiminor_axis(component.elements.a, component.elements.ec);
+	component.pqw_to_bci = physics::orbit::frame::pqw::to_bci
+	(
+		component.elements.om,
+		component.elements.in,
+		component.elements.w
+	);
+}
+
+void orbit::on_orbit_replace(entity::registry& registry, entity::id entity_id, entity::component::orbit& component)
+{
+	component.semiminor_axis = physics::orbit::semiminor_axis(component.elements.a, component.elements.ec);
+	component.pqw_to_bci = physics::orbit::frame::pqw::to_bci
+	(
+		component.elements.om,
+		component.elements.in,
+		component.elements.w
+	);
+}
+
 } // namespace system
 } // namespace entity

src/entity/systems/orbit.hpp

@@ -22,6 +22,8 @@
 
 #include "entity/systems/updatable.hpp"
 #include "utility/fundamental-types.hpp"
+#include "entity/id.hpp"
+#include "entity/components/orbit.hpp"
 
 namespace entity {
 namespace system {
@@ -58,6 +60,9 @@ public:
 	void set_time_scale(double scale);
 	
 private:
+	void on_orbit_construct(entity::registry& registry, entity::id entity_id, entity::component::orbit& component);
+	void on_orbit_replace(entity::registry& registry, entity::id entity_id, entity::component::orbit& component);
+	
 	double universal_time;
 	double time_scale;
 	std::size_t ke_iterations;

src/game/state/nuptial-flight.cpp

@@ -97,14 +97,14 @@ nuptial_flight::nuptial_flight(game::context& ctx):
 		game::world::create_moon(ctx);
 		
 		// Set time to solar noon
-		game::world::set_time(ctx, 0.0);
+		game::world::set_time(ctx, 50.0);//107.3);
 		
 		// Freeze time
 		game::world::set_time_scale(ctx, 0.0);
 	}
 	
 	// Load biome
-	game::load::biome(ctx, "desert-scrub.bio");
+	game::load::biome(ctx, "lab.bio");
 	
 	// Setup and enable sky and ground passes
 	ctx.sky_pass->set_enabled(true);
@@ -124,6 +124,20 @@ nuptial_flight::nuptial_flight(game::context& ctx):
 		entity::command::warp_to(*ctx.entity_registry, ruler_10cm_eid, {0, 0, 10});
 	}
 	
+	// Create cocoon
+	{
+		entity::archetype* cocoon_archetype = ctx.resource_manager->load<entity::archetype>("cocoon.ent");
+		auto cocoon_eid = cocoon_archetype->create(*ctx.entity_registry);
+		entity::command::warp_to(*ctx.entity_registry, cocoon_eid, {-10, 0, 10});
+	}
+	
+	// Create larva
+	{
+		entity::archetype* larva_archetype = ctx.resource_manager->load<entity::archetype>("long-neck-larva.ent");
+		auto larva_eid = larva_archetype->create(*ctx.entity_registry);
+		entity::command::warp_to(*ctx.entity_registry, larva_eid, {-13, 0, 10});
+	}
+	
 	// Create keeper if not yet created
 	if (ctx.entities.find("keeper") == ctx.entities.end())
 	{
@@ -152,7 +166,7 @@ nuptial_flight::nuptial_flight(game::context& ctx):
 		locomotion.yaw = 0.0f;
 		ctx.entity_registry->assign<entity::component::locomotion>(boid_eid, locomotion);
 		
-		entity::command::warp_to(*ctx.entity_registry, boid_eid, {0, 2, 0});
+		entity::command::warp_to(*ctx.entity_registry, boid_eid, {0, 1, 0});
 
 		// Set target ant
 		ctx.entities["ant"] = boid_eid;
@@ -264,7 +278,7 @@ void nuptial_flight::setup_camera()
 		ctx.entity_registry->assign<entity::component::constraint_stack>(camera_eid, constraint_stack);
 	}
 	
-	float ev100 = 14.5f;
+	float ev100 = 15.0f;
 	ctx.surface_camera->set_exposure(ev100);
 }
 

src/game/tools.cpp

@@ -95,7 +95,7 @@ entity::id build_time_tool(game::context& ctx)
 		
 		entity::id planet_eid = ctx.entities["planet"];
 		entity::component::celestial_body body = ctx.entity_registry->get<entity::component::celestial_body>(planet_eid);
-		body.axial_rotation = math::radians(360.0f) * ((float)mouse_x / (float)window_w);
+		//body.axial_rotation = math::radians(360.0f) * ((float)mouse_x / (float)window_w);
 		ctx.entity_registry->replace<entity::component::celestial_body>(planet_eid, body);
 	};
 	

src/game/world.cpp

@@ -95,12 +95,8 @@ void create_stars(game::context& ctx)
 		ra = math::wrap_radians(math::radians(ra));
 		dec = math::wrap_radians(math::radians(dec));
 		
-		// Transform spherical equatorial coordinates to rectangular equatorial coordinates
-		double3 position_bci = geom::spherical::to_cartesian(double3{1.0, dec, ra});
-		
-		// Transform coordinates from equatorial space to inertial space
-		physics::frame<double> bci_to_inertial = physics::orbit::inertial::to_bci({0, 0, 0}, 0.0, math::radians(23.4393)).inverse();
-		double3 position_inertial = bci_to_inertial * position_bci;
+		// Convert ICRF coordinates from spherical to Cartesian
+		double3 position = physics::orbit::frame::bci::cartesian(double3{1.0, dec, ra});
 		
 		// Convert color index to color temperature
 		double cct = color::index::bv_to_cct(bv_color);
@@ -118,9 +114,9 @@ void create_stars(game::context& ctx)
 		double3 scaled_color = color_acescg * brightness;
 		
 		// Build vertex
-		*(star_vertex++) = static_cast<float>(position_inertial.x);
-		*(star_vertex++) = static_cast<float>(position_inertial.y);
-		*(star_vertex++) = static_cast<float>(position_inertial.z);
+		*(star_vertex++) = static_cast<float>(position.x);
+		*(star_vertex++) = static_cast<float>(position.y);
+		*(star_vertex++) = static_cast<float>(position.z);
 		*(star_vertex++) = static_cast<float>(scaled_color.x);
 		*(star_vertex++) = static_cast<float>(scaled_color.y);
 		*(star_vertex++) = static_cast<float>(scaled_color.z);
@@ -188,22 +184,23 @@ void create_sun(game::context& ctx)
 	entity::id sun_eid = sun_archetype->create(*ctx.entity_registry);
 	ctx.entities["sun"] = sun_eid;
 	
-	// Create direct sun light scene object
-	scene::directional_light* sun_direct = new scene::directional_light();
+	// Create sun directional light scene object
+	scene::directional_light* sun_light = new scene::directional_light();
 	
-	// Create ambient sun light scene object
-	scene::ambient_light* sun_ambient = new scene::ambient_light();
-	sun_ambient->set_color({1, 1, 1});
-	sun_ambient->set_intensity(30000.0f * 0.0f);
-	sun_ambient->update_tweens();
+	// Create sky ambient light scene object
+	scene::ambient_light* sky_light = new scene::ambient_light();
+	sky_light->set_color({1, 1, 1});
+	sky_light->set_intensity(0.0f);
+	sky_light->update_tweens();
 	
 	// Add sun light scene objects to surface scene
-	ctx.surface_scene->add_object(sun_direct);
-	ctx.surface_scene->add_object(sun_ambient);
+	ctx.surface_scene->add_object(sun_light);
+	ctx.surface_scene->add_object(sky_light);
 	
 	// Pass direct sun light scene object to shadow map pass and astronomy system
-	ctx.surface_shadow_map_pass->set_light(sun_direct);
-	ctx.astronomy_system->set_sun_light(sun_direct);
+	ctx.surface_shadow_map_pass->set_light(sun_light);
+	ctx.astronomy_system->set_sun_light(sun_light);
+	ctx.astronomy_system->set_sky_light(sky_light);
 }
 
 void create_planet(game::context& ctx)
@@ -213,6 +210,9 @@ void create_planet(game::context& ctx)
 	entity::id planet_eid = planet_archetype->create(*ctx.entity_registry);
 	ctx.entities["planet"] = planet_eid;
 	
+	// Assign orbital parent
+	ctx.entity_registry->get<entity::component::orbit>(planet_eid).parent = ctx.entities["sun"];
+	
 	// Assign planetary terrain component
 	entity::component::terrain terrain;
 	terrain.elevation = [](double, double) -> double
@@ -231,9 +231,13 @@ void create_planet(game::context& ctx)
 void create_moon(game::context& ctx)
 {
 	// Create lunar entity
-	entity::id moon_eid = ctx.entity_registry->create();
+	entity::archetype* moon_archetype = ctx.resource_manager->load<entity::archetype>("moon.ent");
+	entity::id moon_eid = moon_archetype->create(*ctx.entity_registry);
 	ctx.entities["moon"] = moon_eid;
 	
+	// Assign orbital parent
+	ctx.entity_registry->get<entity::component::orbit>(moon_eid).parent = ctx.entities["planet"];
+	
 	// Pass moon model to sky pass
 	ctx.sky_pass->set_moon_model(ctx.resource_manager->load<render::model>("moon.mdl"));
 }

src/math/math.hpp

@@ -35,6 +35,7 @@ namespace math {}
 #include "math/quaternion-functions.hpp"
 #include "math/quaternion-operators.hpp"
 
+#include "math/se3.hpp"
 #include "math/transform-type.hpp"
 #include "math/transform-functions.hpp"
 #include "math/transform-operators.hpp"

src/math/se3.hpp

@@ -0,0 +1,139 @@
+/*
+ * Copyright (C) 2021  Christopher J. Howard
+ *
+ * This file is part of Antkeeper source code.
+ *
+ * Antkeeper source code 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.
+ *
+ * Antkeeper source code 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 Antkeeper source code.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef ANTKEEPER_MATH_TRANSFORMATION_SE3_HPP
+#define ANTKEEPER_MATH_TRANSFORMATION_SE3_HPP
+
+#include "math/vector-type.hpp"
+#include "math/vector-operators.hpp"
+#include "math/quaternion-type.hpp"
+#include "math/quaternion-operators.hpp"
+#include "math/quaternion-functions.hpp"
+
+namespace math {
+namespace transformation {
+
+/**
+ * 3-dimensional Euclidean proper rigid transformation in SE(3).
+ *
+ * @tparam T Scalar type.
+ */
+template <class T>
+struct se3
+{
+public:
+	/// Scalar type.
+	typedef T scalar_type;
+	
+	/// Vector type.
+	typedef math::vector3<T> vector_type;
+	
+	/// Quaternion type.
+	typedef math::quaternion<T> quaternion_type;
+	
+	/// Transformation matrix type.
+	typedef math::matrix<T, 4, 4> matrix_type;
+	
+	/// Vector representing the translation component of this SE(3) transformation.
+	vector_type t;
+	
+	/// Quaternion representing the rotation component of this SE(3) transformation.
+	quaternion_type r;
+	
+	/// Returns the inverse of this SE(3) transformation.
+	se3 inverse() const;
+	
+	/// Returns a matrix representation of the SE(3) transformation.
+	matrix_type matrix() const;
+	
+	/**
+	 * Transforms a vector by this SE(3) transformation.
+	 *
+	 * @param x Untransformed vector.
+	 * @return Transformed vector.
+	 */
+	vector_type transform(const vector_type& x) const;
+	
+	/**
+	 * Transforms an SE(3) transformation by this SE(3) transformation.
+	 *
+	 * @param x Other SE(3) transformation.
+	 * @return Frame in this se3's space.
+	 */
+	se3 transform(const se3& x) const;
+	
+	/// @copydoc se3::transform(const vector_type&) const
+	vector_type operator*(const vector_type& x) const;
+	
+	/// @copydoc se3::transform(const se3&) const
+	se3 operator*(const se3& x) const;
+};
+
+template <class T>
+se3<T> se3<T>::inverse() const
+{
+	const quaternion_type inverse_r = math::conjugate(r);
+	const vector_type inverse_t = -(inverse_r * t);
+	return se3{inverse_t, inverse_r};
+}
+
+template <class T>
+typename se3<T>::matrix_type se3<T>::matrix() const
+{
+	matrix_type m = math::resize<4, 4>(math::matrix_cast<T>(r));
+	
+	m[3].x = t.x;
+	m[3].y = t.y;
+	m[3].z = t.z;
+	
+	return m;
+}
+
+template <class T>
+typename se3<T>::vector_type se3<T>::transform(const vector_type& x) const
+{
+	return r * x + t;
+}
+
+template <class T>
+se3<T> se3<T>::transform(const se3& x) const
+{
+	return se3
+	{
+		x.transform(t),
+		math::normalize(x.r * r)
+	};
+}
+
+template <class T>
+typename se3<T>::vector_type se3<T>::operator*(const vector_type& x) const
+{
+	return transform(x);
+}
+
+template <class T>
+se3<T> se3<T>::operator*(const se3& x) const
+{
+	return transform(x);
+}
+
+} // namespace transformation
+} // namespace math
+
+#endif // ANTKEEPER_MATH_TRANSFORMATION_SE3_HPP

src/physics/atmosphere.hpp

@@ -124,7 +124,7 @@ T extinction(T ec, T s)
 template <class T>
 T albedo(T s, T e)
 {
-	return s / t;
+	return s / e;
 }
 
 /**

src/physics/frame.hpp

@@ -1,134 +0,0 @@
-/*
- * Copyright (C) 2021  Christopher J. Howard
- *
- * This file is part of Antkeeper source code.
- *
- * Antkeeper source code 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.
- *
- * Antkeeper source code 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 Antkeeper source code.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#ifndef ANTKEEPER_PHYSICS_FRAME_HPP
-#define ANTKEEPER_PHYSICS_FRAME_HPP
-
-#include "math/math.hpp"
-#include <cmath>
-
-namespace physics {
-
-/**
- * 3-dimensional frame of reference.
- *
- * @tparam T Scalar type.
- */
-template <class T>
-struct frame
-{
-public:
-	/// Scalar type.
-	typedef T scalar_type;
-	
-	/// Vector type.
-	typedef math::vector3<T> vector_type;
-	
-	/// Quaternion type.
-	typedef math::quaternion<T> quaternion_type;
-	
-	/// Transformation matrix type.
-	typedef math::matrix<T, 4, 4> matrix_type;
-	
-	/// Vector representing a translation from the origin of this frame to the origin of the parent frame.	
-	vector_type translation;
-	
-	/// Quaternion representing a rotation from the orientation of this frame to the orientation of the parent frame.
-	quaternion_type rotation;
-	
-	/// Returns the inverse of this frame.
-	frame inverse() const;
-	
-	/// Returns a transformation matrix representation of this frame.
-	matrix_type matrix() const;
-	
-	/**
-	 * Transforms a vector from the parent frame space to this frame's space.
-	 *
-	 * @param v Vector in parent frame space.
-	 * @return Vector in this frame's space.
-	 */
-	vector_type transform(const vector_type& v) const;
-	
-	/**
-	 * Transforms a frame from the parent frame space to this frame's space.
-	 *
-	 * @param f Frame in parent frame space.
-	 * @return Frame in this frame's space.
-	 */
-	frame transform(const frame& f) const;
-	
-	/// @copydoc frame::transform(const vector_type&) const
-	vector_type operator*(const vector_type& v) const;
-	
-	/// @copydoc frame::transform(const frame&) const
-	frame operator*(const frame& f) const;
-};
-
-template <class T>
-frame<T> frame<T>::inverse() const
-{
-	const quaternion_type inverse_rotation = math::conjugate(rotation);
-	const vector_type inverse_translation = -(inverse_rotation * translation);
-	return frame{inverse_translation, inverse_rotation};
-}
-
-template <class T>
-typename frame<T>::matrix_type frame<T>::matrix() const
-{
-	matrix_type m = math::resize<4, 4>(math::matrix_cast<T>(rotation));
-	
-	m[3].x = translation.x;
-	m[3].y = translation.y;
-	m[3].z = translation.z;
-	
-	return m;
-}
-
-template <class T>
-typename frame<T>::vector_type frame<T>::transform(const vector_type& v) const
-{
-	return translation + rotation * v;
-}
-
-template <class T>
-frame<T> frame<T>::transform(const frame& f) const
-{
-	return frame
-	{
-		f.transform(translation),
-		math::normalize(f.rotation * rotation)
-	};
-}
-
-template <class T>
-typename frame<T>::vector_type frame<T>::operator*(const vector_type& v) const
-{
-	return transform(v);
-}
-
-template <class T>
-frame<T> frame<T>::operator*(const frame& f) const
-{
-	return transform(f);
-}
-
-} // namespace physics
-
-#endif // ANTKEEPER_PHYSICS_FRAME_HPP

src/physics/orbit/anomaly.hpp

@@ -0,0 +1,199 @@
+/*
+ * Copyright (C) 2021  Christopher J. Howard
+ *
+ * This file is part of Antkeeper source code.
+ *
+ * Antkeeper source code 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.
+ *
+ * Antkeeper source code 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 Antkeeper source code.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef ANTKEEPER_PHYSICS_ORBIT_ANOMALY_HPP
+#define ANTKEEPER_PHYSICS_ORBIT_ANOMALY_HPP
+
+#include "math/constants.hpp"
+#include <cmath>
+
+namespace physics {
+namespace orbit {
+
+/**
+ * Orbital anomaly functions.
+ */
+namespace anomaly {
+
+/**
+ * Derives the eccentric anomaly given eccentricity and true anomaly.
+ *
+ * @param ec Eccentricity (e).
+ * @param ta True anomaly (nu).
+ * @return Eccentric anomaly (E).
+ */
+template <class T>
+T true_to_eccentric(T ec, T ta);
+
+/**
+ * Derives the mean anomaly given eccentricity and true anomaly.
+ *
+ * @param ec Eccentricity (e).
+ * @param ta True anomaly (nu).
+ * @return Mean anomaly (M).
+ */
+template <class T>
+T true_to_mean(T ec, T ta);
+
+/**
+ * Derives the true anomaly given eccentricity and eccentric anomaly.
+ *
+ * @param ec Eccentricity (e).
+ * @param ea Eccentric anomaly (E).
+ * @return True anomaly (nu).
+ */
+template <class T>
+T eccentric_to_true(T ec, T ea);
+
+/**
+ * Derives the mean anomaly given eccentricity and eccentric anomaly.
+ *
+ * @param ec Eccentricity (e).
+ * @param ea Eccentric anomaly (E).
+ * @return Mean anomaly (M).
+ */
+template <class T>
+T eccentric_to_mean(T ec, T ea);
+
+/**
+ * Iteratively derives the eccentric anomaly given eccentricity and mean anomaly.
+ *
+ * @param ec Eccentricity (e).
+ * @param ma Mean anomaly (M).
+ * @param iterations Maximum number of iterations.
+ * @param tolerance Solution error tolerance.
+ * @return Eccentric anomaly (E).
+ *
+ * @see Murison, Marc. (2006). A Practical Method for Solving the Kepler Equation. 10.13140/2.1.5019.6808.
+ */
+template <class T>
+T mean_to_eccentric(T ec, T ma, std::size_t iterations, T tolerance);
+
+/**
+ * Iteratively derives the true anomaly given eccentricity and mean anomaly.
+ *
+ * @param ec Eccentricity (e).
+ * @param ma Mean anomaly (M).
+ * @param iterations Maximum number of iterations.
+ * @param tolerance Solution error tolerance.
+ * @return True anomaly (nu).
+ */
+template <class T>
+T mean_to_true(T ec, T ma, std::size_t iterations, T tolerance);
+
+template <class T>
+T true_to_eccentric(T ec, T ta)
+{
+	// Parabolic orbit
+	if (ec == T(1))
+		return std::tan(ta * T(0.5));
+	
+	// Hyperbolic orbit
+	if (ec > T(1))
+		return std::acosh((ec + std::cos(ta)) / (T(1) + ec * std::cos(ta))) * ((ta < T(0)) ? T(-1) : T(1));
+	
+	// Elliptic orbit
+	return std::atan2(std::sqrt(T(1) - ec * ec) * std::sin(ta), std::cos(ta) + ec);
+}
+
+template <class T>
+T true_to_mean(T ec, T ta)
+{
+    return eccentric_to_mean(ec, true_to_eccentric(ec, ta));
+}
+
+template <class T>
+T eccentric_to_true(T ec, T ea)
+{
+	// Parabolic orbit
+	if (ec == T(1))
+		return std::atan(ea) * T(2);
+	
+	// Hyperbolic orbit
+	if (ec > T(1))
+		return std::atan(std::sqrt((ec + T(1)) / (ec - T(1))) * std::tanh(ea * T(0.5))) * T(2);
+	
+	// Elliptic orbit
+	return std::atan2(sqrt(T(1) - ec * ec) * std::sin(ea), std::cos(ea) - ec);
+}
+
+template <class T>
+T eccentric_to_mean(T ec, T ea)
+{
+	// Parabolic orbit
+	if (ec == T(1))
+		return (ea * ea * ea) /  T(6) + ea * T(0.5);
+	
+	// Hyperbolic orbit
+	if (ec > T(1))
+		return ec * std::sinh(ea) - ea;
+	
+	// Elliptic orbit
+	return ea - ec * std::sin(ea);
+}
+
+template <class T>
+T mean_to_eccentric(T ec, T ma, std::size_t iterations, T tolerance)
+{
+	// Wrap mean anomaly to `[-pi, pi]`
+	ma = std::remainder(ma, math::two_pi<T>);
+	
+	// Third-order approximation of eccentric anomaly starting value, E0
+	const T t33 = std::cos(ma);
+	const T t34 = ec * ec;
+	const T t35 = t34 * ec;
+	T ea0 = ma + (T(-0.5) * t35 + ec + (t34 + T(1.5) * t33 * t35) * t33) * std::sin(ma);
+	
+	// Iteratively converge E0 and E1
+	for (std::size_t i = 0; i < iterations; ++i)
+	{
+		// Third-order approximation of eccentric anomaly, E1
+		const T t1 = std::cos(ea0);
+		const T t2 = T(-1) + ec * t1;
+		const T t3 = std::sin(ea0);
+		const T t4 = ec * t3;
+		const T t5 = -ea0 + t4 + ma;
+		const T t6 = t5 / (T(0.5) * t5 * t4 / t2 + t2);
+		const T ea1 = ea0 - (t5 / ((T(0.5) * t3 - (T(1) / T(6)) * t1 * t6) * ec * t6 + t2));
+		
+		// Determine solution error
+		const T error = std::abs(ea1 - ea0);
+		
+		// Set E0 to E1
+		ea0 = ea1;
+		
+		// Break if solution is within error tolerance
+		if (error < tolerance)
+			break;
+	}
+	
+	return ea0;
+}
+
+template <class T>
+T mean_to_true(T ec, T ma, std::size_t iterations, T tolerance)
+{
+	eccentric_to_true(ec, mean_to_eccentric(ec, ma, iterations, tolerance));
+}
+
+} // namespace anomaly
+} // namespace orbit
+} // namespace physics
+
+#endif // ANTKEEPER_PHYSICS_ORBIT_ANOMALY_HPP

src/physics/orbit/elements.hpp

@@ -21,7 +21,7 @@
 #define ANTKEEPER_PHYSICS_ORBIT_ELEMENTS_HPP
 
 #include "utility/fundamental-types.hpp"
-#include "physics/orbit/state.hpp"
+#include "math/constants.hpp"
 #include <cmath>
 
 namespace physics {
@@ -39,54 +39,133 @@ struct elements
 	typedef T scalar_type;
 	
 	/// Eccentricity (e).
-	scalar_type e;
+	scalar_type ec;
 	
 	/// Semimajor axis (a).
 	scalar_type a;
 	
 	/// Inclination (i), in radians.
-	scalar_type i;
+	scalar_type in;
 	
 	/// Right ascension of the ascending node (OMEGA), in radians.
-	scalar_type raan;
+	scalar_type om;
 	
 	/// Argument of periapsis (omega), in radians.
 	scalar_type w;
 	
-	/// True anomaly (nu) at epoch, in radians.
-	scalar_type ta;
+	/// Mean anomaly (M) at epoch, in radians.
+	scalar_type ma;
 };
 
+/**
+ * Calculates the period of an elliptical orbit according to Kepler's third law.
+ *
+ * @param a Semimajor axis (a).
+ * @param gm Standard gravitational parameter (GM).
+ * @return Orbital period (T).
+ */
+template <class T>
+T period(T a, T gm);
+
+/**
+ * Calculates the mean motion (n) of an orbit.
+ *
+ * @param a Semimajor axis (a).
+ * @param gm Standard gravitational parameter (GM).
+ * @return Mean motion (n), in radians per unit time.
+ */
+template <class T>
+T mean_motion(T a, T gm);
+
+/**
+ * Calculates the mean motion (n) of an orbit.
+ *
+ * @param t Orbital period (T).
+ * @return Mean motion (n), in radians per unit time.
+ */
+template <class T>
+T mean_motion(T t);
+
+/**
+ * Derives the argument of the periapsis (omega) of an orbit, given the longitude of periapsis (pomega) and longitude of the ascending node (OMEGA).
+ *
+ * @param lp Longitude of the periapsis (pomega), in radians.
+ * @param om Right ascension of the ascending node (OMEGA), in radians.
+ * @return Argument of the periapsis (omega), in radians.
+ */
+template <class T>
+T argument_periapsis(T om, T lp);
+
 /**
  * Derives the longitude of the periapsis (pomega) of an orbit, given the argument of periapsis (omega) and longitude of the ascending node (OMEGA).
  *
- * @param w Argument of periapsis (omega), in radians.
- * @param raan Right ascension of the ascending node (OMEGA), in radians.
+ * @param w Argument of the periapsis (omega), in radians.
+ * @param om Right ascension of the ascending node (OMEGA), in radians.
  * @return Longitude of the periapsis (pomega), in radians.
  */
 template <class T>
-T derive_longitude_periapsis(T w, T raan);
+T longitude_periapsis(T om, T w);
 
 /**
  * Derives the semiminor axis (b) of an orbit, given the semimajor axis (a) and eccentricity (e).
  *
  * @param a Semimajor axis (a).
- * @param e Eccentricity (e).
+ * @param ec Eccentricity (e).
  * @return Semiminor axis (b).
  */
 template <class T>
-T derive_semiminor_axis(T a, T e);
+T semiminor_axis(T a, T ec);
+
+/**
+ * Derives the semi-latus rectum (l) of an orbit, given the semimajor axis (a) and eccentricity (e).
+ *
+ * @param a Semimajor axis (a).
+ * @param ec Eccentricity (e).
+ * @return Semi-latus rectum (l).
+ */
+template <class T>
+T semilatus_rectum(T a, T ec);
+
+template <class T>
+T period(T a, T gm)
+{
+	return math::two_pi<T> * std::sqrt((a * a * a) / gm);
+}
+
+template <class T>
+T mean_motion(T a, T gm)
+{
+	return std::sqrt((a * a * a) / gm);
+}
+
+template <class T>
+T mean_motion(T t)
+{
+	return math::two_pi<T> / t;
+}
+
+template <class T>
+T argument_periapsis(T om, T lp)
+{
+	return lp - om;
+}
+
+template <class T>
+T longitude_periapsis(T om, T w)
+{
+	return w + om;
+}
 
 template <class T>
-T derive_longitude_periapsis(T w, T raan)
+T semiminor_axis(T a, T ec)
 {
-	return w + raan;
+	return a * std::sqrt(T(1) - ec * ec);
 }
 
 template <class T>
-T derive_semiminor_axis(T a, T e)
+T semilatus_rectum(T a, T ec)
 {
-	return a * std::sqrt(T(1) - e * e);
+	return a * (T(1) - ec * ec);
 }
 
 } // namespace orbit

src/physics/orbit/frame.hpp

@@ -0,0 +1,399 @@
+/*
+ * Copyright (C) 2021  Christopher J. Howard
+ *
+ * This file is part of Antkeeper source code.
+ *
+ * Antkeeper source code 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.
+ *
+ * Antkeeper source code 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 Antkeeper source code.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef ANTKEEPER_PHYSICS_ORBIT_FRAME_HPP
+#define ANTKEEPER_PHYSICS_ORBIT_FRAME_HPP
+
+#include "math/se3.hpp"
+#include <cmath>
+
+namespace physics {
+namespace orbit {
+
+/// Orbital reference frames.
+namespace frame {
+
+/// Perifocal (PQW) frame.
+namespace pqw {
+	
+	/**
+	 * Converts PQW coordinates from Cartesian to spherical.
+	 *
+	 * @param v PQW Cartesian coordinates.
+	 * @return PQW spherical coordinates, in the ISO order of radial distance, inclination (radians), and true anomaly (radians).
+	 */
+	template <class T>
+	math::vector3<T> spherical(const math::vector3<T>& v)
+	{
+		const T xx_yy = v.x * v.x + v.y * v.y;
+		
+		return math::vector3<T>
+		{
+			std::sqrt(xx_yy + v.z * v.z),
+			std::atan2(v.z, std::sqrt(xx_yy)),
+			std::atan2(v.y, v.x)
+		};
+	}
+	
+	/**
+	 * Constructs spherical PQW coordinates from Keplerian orbital elements.
+	 *
+	 * @param ec Eccentricity (e).
+	 * @param a Semimajor axis (a).
+	 * @param ea Eccentric anomaly (E), in radians.
+	 * @param b Semiminor axis (b).
+	 * @return PQW spherical coordinates, in the ISO order of radial distance, inclination (radians), and true anomaly (radians).
+	 */
+	template <class T>
+	math::vector3<T> spherical(T ec, T a, T ea, T b)
+	{
+		const T x = a * (std::cos(ea) - ec);
+		const T y = b * std::sin(ea);
+		const T d = std::sqrt(x * x + y * y);
+		const T ta = std::atan2(y, x);
+		
+		return math::vector3<T>{d, T(0), ta};
+	}
+	
+	/**
+	 * Constructs spherical PQW coordinates from Keplerian orbital elements.
+	 *
+	 * @param ec Eccentricity (e).
+	 * @param a Semimajor axis (a).
+	 * @param ea Eccentric anomaly (E), in radians.
+	 * @return PQW spherical coordinates, in the ISO order of radial distance, inclination (radians), and true anomaly (radians).
+	 */
+	template <class T>
+	math::vector3<T> spherical(T ec, T a, T ea)
+	{
+		const T b = a * std::sqrt(T(1) - ec * ec);
+		return spherical<T>(ec, a, ea, b);
+	}
+	
+	/**
+	 * Converts PQW coordinates from spherical to Cartesian.
+	 *
+	 * @param PQW spherical coordinates, in the ISO order of radial distance, inclination (radians), and true anomaly (radians).
+	 * @return PQW Cartesian coordinates.
+	 */
+	template <class T>
+	math::vector3<T> cartesian(const math::vector3<T>& v)
+	{
+		const T x = v[0] * std::cos(v[1]);
+		
+		return math::vector3<T>
+		{
+			x * std::cos(v[2]),
+			x * std::sin(v[2]),
+			v[0] * std::sin(v[1]),
+		};
+	}
+	
+	/**
+	 * Constructs Cartesian PQW coordinates from Keplerian orbital elements.
+	 *
+	 * @param ec Eccentricity (e).
+	 * @param a Semimajor axis (a).
+	 * @param ea Eccentric anomaly (E), in radians.
+	 * @param b Semiminor axis (b).
+	 * @return PQW Cartesian coordinates.
+	 */
+	template <class T>
+	math::vector3<T> cartesian(T ec, T a, T ea, T b)
+	{
+		return cartesian<T>(spherical<T>(ec, a, ea, b));
+	}
+	
+	/**
+	 * Constructs Cartesian PQW coordinates from Keplerian orbital elements.
+	 *
+	 * @param ec Eccentricity (e).
+	 * @param a Semimajor axis (a).
+	 * @param ea Eccentric anomaly (E), in radians.
+	 * @return PQW Cartesian coordinates.
+	 */
+	template <class T>
+	math::vector3<T> cartesian(T ec, T a, T ea)
+	{
+		return cartesian<T>(spherical<T>(ec, a, ea));
+	}
+
+	/**
+	 * Constructs an SE(3) transformation from a PQW frame to a BCI frame.
+	 *
+	 * @param om Right ascension of the ascending node (OMEGA), in radians.
+	 * @param in Orbital inclination (i), in radians.
+	 * @param w Argument of periapsis (omega), in radians.
+	 * @return PQW to BCI transformation.
+	 */
+	template <typename T>
+	math::transformation::se3<T> to_bci(T om, T in, T w)
+	{
+		const math::quaternion<T> r = math::normalize
+		(
+			math::quaternion<T>::rotate_z(om) *
+				math::quaternion<T>::rotate_x(in) *
+				math::quaternion<T>::rotate_z(w)
+		);
+		
+		return math::transformation::se3<T>{{T(0), T(0), T(0)}, r};
+	}
+
+} // namespace pqw
+
+/// Body-centered inertial (BCI) frame.
+namespace bci {
+	
+	/**
+	 * Converts BCI coordinates from spherical to Cartesian.
+	 *
+	 * @param v BCI spherical coordinates, in the ISO order of radial distance, declination (radians), and right ascension (radians).
+	 * @return BCI Cartesian coordinates.
+	 */
+	template <class T>
+	math::vector3<T> cartesian(const math::vector3<T>& v)
+	{
+		const T x = v[0] * std::cos(v[1]);
+		
+		return math::vector3<T>
+		{
+			x * std::cos(v[2]),
+			x * std::sin(v[2]),
+			v[0] * std::sin(v[1]),
+		};
+	}
+	
+	/**
+	 * Converts BCI coordinates from Cartesian to spherical.
+	 *
+	 * @param v BCI Cartesian coordinates.
+	 * @return BCI spherical coordinates, in the ISO order of radial distance, declination (radians), and right ascension (radians).
+	 */
+	template <class T>
+	math::vector3<T> spherical(const math::vector3<T>& v)
+	{
+		const T xx_yy = v.x * v.x + v.y * v.y;
+		
+		return math::vector3<T>
+		{
+			std::sqrt(xx_yy + v.z * v.z),
+			std::atan2(v.z, std::sqrt(xx_yy)),
+			std::atan2(v.y, v.x)
+		};
+	}
+	
+	/**
+	 * Constructs an SE(3) transformation from a BCI frame to a BCBF frame.
+	 *
+	 * @param ra Right ascension of the north pole, in radians.
+	 * @param dec Declination of the north pole, in radians.
+	 * @param w Location of the prime meridian, as a rotation about the north pole, in radians.
+	 * @return BCI to BCBF transformation.
+	 */
+	template <typename T>
+	math::transformation::se3<T> to_bcbf(T ra, T dec, T w)
+	{
+		const math::quaternion<T> r = math::normalize
+		(
+			math::quaternion<T>::rotate_z(-w) *
+				math::quaternion<T>::rotate_x(-math::half_pi<T> + dec) *
+				math::quaternion<T>::rotate_z(-math::half_pi<T> - ra)
+		);
+		
+		return math::transformation::se3<T>{{T(0), T(0), T(0)}, r};
+	}
+	
+	/**
+	 * Constructs an SE(3) transformation from a BCI frame to a PQW frame.
+	 *
+	 * @param om Right ascension of the ascending node (OMEGA), in radians.
+	 * @param in Orbital inclination (i), in radians.
+	 * @param w Argument of periapsis (omega), in radians.
+	 * @return BCI to PQW transformation.
+	 */
+	template <typename T>
+	math::transformation::se3<T> to_pqw(T om, T in, T w)
+	{
+		const math::quaternion<T> r = math::normalize
+		(
+			math::quaternion<T>::rotate_z(-w) *
+				math::quaternion<T>::rotate_x(-in) *
+				math::quaternion<T>::rotate_z(-om)
+		);
+		
+		return math::transformation::se3<T>{{T(0), T(0), T(0)}, r};
+	}
+	
+} // namespace bci
+
+/// Body-centered, body-fixed (BCBF) frame.
+namespace bcbf {
+	
+	/**
+	 * Converts BCBF coordinates from spherical to Cartesian.
+	 *
+	 * @param v BCBF spherical coordinates, in the ISO order of radial distance, latitude (radians), and longitude (radians).
+	 * @return BCBF Cartesian coordinates.
+	 */
+	template <class T>
+	math::vector3<T> cartesian(const math::vector3<T>& v)
+	{
+		const T x = v[0] * std::cos(v[1]);
+		
+		return math::vector3<T>
+		{
+			x * std::cos(v[2]),
+			x * std::sin(v[2]),
+			v[0] * std::sin(v[1]),
+		};
+	}
+	
+	/**
+	 * Converts BCBF coordinates from Cartesian to spherical.
+	 *
+	 * @param v BCBF Cartesian coordinates.
+	 * @return BCBF spherical coordinates, in the ISO order of radial distance, latitude (radians), and longitude (radians).
+	 */
+	template <class T>
+	math::vector3<T> spherical(const math::vector3<T>& v)
+	{
+		const T xx_yy = v.x * v.x + v.y * v.y;
+		
+		return math::vector3<T>
+		{
+			std::sqrt(xx_yy + v.z * v.z),
+			std::atan2(v.z, std::sqrt(xx_yy)),
+			std::atan2(v.y, v.x)
+		};
+	}
+	
+	/**
+	 * Constructs an SE(3) transformation from a BCBF frame to a BCI frame.
+	 *
+	 * @param ra Right ascension of the north pole, in radians.
+	 * @param dec Declination of the north pole, in radians.
+	 * @param w Location of the prime meridian, as a rotation about the north pole, in radians.
+	 * @return BCBF to BCI transformation.
+	 */
+	template <typename T>
+	math::transformation::se3<T> to_bci(T ra, T dec, T w)
+	{
+		const math::quaternion<T> r = math::normalize
+		(
+			math::quaternion<T>::rotate_z(math::half_pi<T> + ra) *
+				math::quaternion<T>::rotate_x(math::half_pi<T> - dec) *
+				math::quaternion<T>::rotate_z(w)
+		);
+		
+		return math::transformation::se3<T>{{T(0), T(0), T(0)}, r};
+	}
+	
+	/**
+	 * Constructs an SE(3) transformation from a BCBF frame to an ENU frame.
+	 *
+	 * @param distance Radial distance of the observer from the center of the body.
+	 * @param latitude Latitude of the observer, in radians.
+	 * @param longitude Longitude of the observer, in radians.
+	 * @return BCBF to ENU transformation.
+	 */
+	template <typename T>
+	math::transformation::se3<T> to_enu(T distance, T latitude, T longitude)
+	{
+		const math::vector3<T> t = {T(0), T(0), -distance};
+		const math::quaternion<T> r = math::normalize
+		(
+			math::quaternion<T>::rotate_x(-math::half_pi<T> + latitude) *
+				math::quaternion<T>::rotate_z(-longitude - math::half_pi<T>)
+		);
+		
+		return math::transformation::se3<T>{t, r};
+	}
+
+} // namespace bcbf
+
+/// East, North, Up (ENU) horizontal frame.
+namespace enu {
+	
+	/**
+	 * Converts ENU coordinates from spherical to Cartesian.
+	 *
+	 * @param v ENU spherical coordinates, in the ISO order of radial distance, elevation (radians), and azimuth (radians).
+	 * @return ENU Cartesian coordinates.
+	 */
+	template <class T>
+	math::vector3<T> cartesian(const math::vector3<T>& v)
+	{
+		const T x = v[0] * std::cos(v[1]);
+		const T y = math::half_pi<T> - v[2];
+		
+		return math::vector3<T>
+		{
+			x * std::cos(y),
+			x * std::sin(y),
+			v[0] * std::sin(v[1]),
+		};
+	}
+	
+	/**
+	 * Converts ENU coordinates from Cartesian to spherical.
+	 *
+	 * @param v ENU Cartesian coordinates.
+	 * @return ENU spherical coordinates, in the ISO order of radial distance, elevation (radians), and azimuth (radians).
+	 */
+	template <class T>
+	math::vector3<T> spherical(const math::vector3<T>& v)
+	{
+		const T xx_yy = v.x * v.x + v.y * v.y;
+		
+		return math::vector3<T>
+		{
+			std::sqrt(xx_yy + v.z * v.z),
+			std::atan2(v.z, std::sqrt(xx_yy)),
+			math::half_pi<T> - std::atan2(v.y, v.x)
+		};
+	}
+	
+	/**
+	 * Constructs an SE(3) transformation from an ENU frame to a BCBF frame.
+	 *
+	 * @param distance Radial distance of the observer from the center of the body.
+	 * @param latitude Latitude of the observer, in radians.
+	 * @param longitude Longitude of the observer, in radians.
+	 * @return ENU to BCBF transformation.
+	 */
+	template <typename T>
+	math::transformation::se3<T> to_bcbf(T distance, T latitude, T longitude)
+	{
+		const math::vector3<T> t = {T(0), T(0), distance};
+		const math::quaternion<T> r = math::normalize
+		(
+			math::quaternion<T>::rotate_z(longitude + math::half_pi<T>) *
+				math::quaternion<T>::rotate_x(math::half_pi<T> - latitude)
+		);
+		
+		return math::transformation::se3<T>{r * t, r};
+	}
+	
+} // namespace enu
+
+} // namespace frame
+} // namespace orbit
+} // namespace physics
+
+#endif // ANTKEEPER_PHYSICS_ORBIT_FRAME_HPP

src/physics/orbit/frames.hpp

@@ -1,131 +0,0 @@
-/*
- * Copyright (C) 2021  Christopher J. Howard
- *
- * This file is part of Antkeeper source code.
- *
- * Antkeeper source code 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.
- *
- * Antkeeper source code 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 Antkeeper source code.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#ifndef ANTKEEPER_PHYSICS_ORBIT_FRAMES_HPP
-#define ANTKEEPER_PHYSICS_ORBIT_FRAMES_HPP
-
-#include "physics/frame.hpp"
-
-namespace physics {
-namespace orbit {
-
-/// Inertial right-handed coordinate system
-namespace inertial {
-
-	/**
-	 * Constucts a reference frame which transforms coordinates from inertial space into perifocal space.
-	 *
-	 * @param focus Cartesian coordinates of the focus of the orbit, in the parent space.
-	 * @param raan Right ascension of the ascending node (OMEGA), in radians.
-	 * @param i Orbital inclination (i), in radians.
-	 * @param w Argument of periapsis (omega), in radians.
-	 * @return Perifocal frame.
-	 */
-	template <typename T>
-	physics::frame<T> to_perifocal(const math::vector3<T>& focus, T raan, T i, T w)
-	{
-		const math::quaternion<T> rotation = math::normalize
-		(
-			math::quaternion<T>::rotate_z(raan) *
-				math::quaternion<T>::rotate_x(i) *
-				math::quaternion<T>::rotate_z(w)
-		);
-		
-		return physics::frame<T>{focus, rotation}.inverse();
-	}
-	
-	/**
-	 * Constructs a reference frame which transforms coordinates from inertial space to body-centered inertial space.
-	 *
-	 * @param r Cartesian position vector (r) of the center of the body, in inertial space.
-	 * @param i Orbital inclination (i), in radians.
-	 * @param axial_tilt Angle between the body's rotational axis and its orbital axis, in radians.
-	 * @return Body-centered inertial frame.
-	 */
-	template <typename T>
-	physics::frame<T> to_bci(const math::vector3<T>& r, T i, T axial_tilt)
-	{
-		return physics::frame<T>{r, math::quaternion<T>::rotate_x(-axial_tilt - i)}.inverse();
-	}
-	
-	/**
-	 * Constructs a reference frame which transforms coordinates from inertial space to body-centered, body-fixed space.
-	 *
-	 * @param r Cartesian position vector (r) of the center of the body, in inertial space.
-	 * @param i Orbital inclination (i), in radians.
-	 * @param axial_tilt Angle between the orbiting body's rotational axis and its orbital axis, in radians.
-	 * @param axial_rotation Angle of rotation about the orbiting body's rotational axis, in radians.
-	 */
-	template <typename T>
-	frame<T> to_bcbf(const math::vector3<T>& r, T i, T axial_tilt, T axial_rotation)
-	{
-		const math::quaternion<T> rotation = math::normalize
-		(
-			math::quaternion<T>::rotate_x(-axial_tilt - i) *
-				math::quaternion<T>::rotate_z(axial_rotation)
-		);
-		
-		return physics::frame<T>{r, rotation}.inverse();
-	}
-
-} // namespace inertial
-
-/// Perifocal right-handed coordinate system in which the x-axis points toward the periapsis of the orbit, the y-axis has a true anomaly of 90 degrees past the periapsis, and the z-axis is the angular momentum vector, which is orthogonal to the orbital plane.
-namespace perifocal {
-
-
-
-} // namespace perifocal
-
-/// Non-inertial body-centered, body-fixed right-handed coordinate system. The x-axis is orthogonal to the intersection of the prime meridian and the equator. The z-axis points toward the positive pole. The y-axis is right-hand orthogonal to the xz-plane.
-namespace bcbf {
-
-	/**
-	 * Constructs a reference frame which transforms coordinates from BCBF space to topocentric space.
-	 *
-	 * @param distance Radial distance of the observer from the center of the body.
-	 * @param latitude Latitude of the observer, in radians.
-	 * @param longitude Longitude of the obserer, in radians.
-	 * @return Topocentric frame.
-	 */
-	template <typename T>
-	physics::frame<T> to_topocentric(T distance, T latitude, T longitude)
-	{
-		const math::vector3<T> translation = {T(0), T(0), distance};
-		
-		const math::quaternion<T> rotation = math::normalize
-		(
-			math::quaternion<T>::rotate_z(longitude) *
-				math::quaternion<T>::rotate_y(math::half_pi<T> - latitude)
-		);
-		
-		return physics::frame<T>{rotation * translation, rotation}.inverse();
-	}
-
-} // namespace bcbf
-
-/// Non-inertial topocentric right-handed coordinate system. Topocentric frames are constructed as south-east-zenith (SEZ) frames; the x-axis points south, the y-axis points east, and the z-axis points toward the zenith (orthogonal to reference spheroid).
-namespace topocentric {
-	
-} // namespace topocentric
-
-} // namespace orbit
-} // namespace physics
-
-#endif // ANTKEEPER_PHYSICS_ORBIT_FRAMES_HPP

src/physics/orbit/kepler.hpp

@@ -1,79 +0,0 @@
-/*
- * Copyright (C) 2021  Christopher J. Howard
- *
- * This file is part of Antkeeper source code.
- *
- * Antkeeper source code 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.
- *
- * Antkeeper source code 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 Antkeeper source code.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#ifndef ANTKEEPER_PHYSICS_ORBIT_KEPLER_HPP
-#define ANTKEEPER_PHYSICS_ORBIT_KEPLER_HPP
-
-#include <cmath>
-
-namespace physics {
-namespace orbit {
-
-/**
- * Iteratively solves Kepler's equation for eccentric anomaly (E).
- *
- * @param ec Eccentricity (e).
- * @param ma Mean anomaly (M).
- * @param iterations Maximum number of iterations.
- * @param tolerance Solution error tolerance.
- * @return Eccentric anomaly (E).
- */
-template <class T>
-T kepler_ea(T ec, T ma, std::size_t iterations, T tolerance = T(0));
-
-/**
- * Solves Kepler's equation for mean anomaly (M).
- *
- * @param ec Eccentricity (e).
- * @param ea Eccentric anomaly (E).
- * @return Mean anomaly (M).
- */
-template <class T>
-T kepler_ma(T ec, T ea);
-
-template <class T>
-T kepler_ea(T ec, T ma, std::size_t iterations, T tolerance)
-{
-	// Approximate eccentric anomaly, E
-	T ea0 = ma + ec * std::sin(ma) * (T(1.0) + ec * std::cos(ma));
-	
-	// Iteratively converge E0 and E1
-	for (std::size_t i = 0; i < iterations; ++i)
-	{
-		const T ea1 = ea0 - (ea0 - ec * std::sin(ea0) - ma) / (T(1.0) - ec * std::cos(ea0));
-		const T error = std::abs(ea1 - ea0);
-		ea0 = ea1;
-		
-		if (error < tolerance)
-			break;
-	}
-	
-	return ea0;
-}
-
-template <class T>
-T kepler_ma(T ec, T ea)
-{
-	return ea - ec * std::sin(ea);
-}
-
-} // namespace orbit
-} // namespace physics
-
-#endif // ANTKEEPER_PHYSICS_ORBIT_KEPLER_HPP

src/physics/orbit/orbit.hpp

@@ -31,9 +31,9 @@ namespace orbit {}
 
 } // namespace physics
 
+#include "physics/orbit/anomaly.hpp"
 #include "physics/orbit/elements.hpp"
-#include "physics/orbit/frames.hpp"
-#include "physics/orbit/kepler.hpp"
+#include "physics/orbit/frame.hpp"
 #include "physics/orbit/state.hpp"
 
 #endif // ANTKEEPER_PHYSICS_ORBIT_HPP

src/render/passes/sky-pass.cpp

@@ -71,10 +71,10 @@ sky_pass::sky_pass(gl::rasterizer* rasterizer, const gl::framebuffer* framebuffe
 	cloud_shader_program(nullptr),
 	observer_altitude_tween(0.0f, math::lerp<float, float>),
 	sun_position_tween(float3{1.0f, 0.0f, 0.0f}, math::lerp<float3, float>),
-	sun_color_outer_tween(float3{1.0f, 1.0f, 1.0f}, math::lerp<float3, float>),
-	sun_color_inner_tween(float3{1.0f, 1.0f, 1.0f}, math::lerp<float3, float>),
-	topocentric_frame_translation({0, 0, 0}, math::lerp<float3, float>),
-	topocentric_frame_rotation(math::quaternion<float>::identity(), math::nlerp<float>)
+	sun_illuminance_outer_tween(float3{1.0f, 1.0f, 1.0f}, math::lerp<float3, float>),
+	sun_illuminance_inner_tween(float3{1.0f, 1.0f, 1.0f}, math::lerp<float3, float>),
+	icrf_to_eus_translation({0, 0, 0}, math::lerp<float3, float>),
+	icrf_to_eus_rotation(math::quaternion<float>::identity(), math::nlerp<float>)
 {}
 
 sky_pass::~sky_pass()
@@ -109,20 +109,20 @@ void sky_pass::render(const render::context& ctx, render::queue& queue) const
 	// Interpolate observer altitude
 	float observer_altitude = observer_altitude_tween.interpolate(ctx.alpha);
 	
-	// Construct tweened inertial to topocentric frame
-	physics::frame<float> topocentric_frame =
+	// Construct tweened ICRF to EUS transformation
+	math::transformation::se3<float> icrf_to_eus =
 	{
-		topocentric_frame_translation.interpolate(ctx.alpha),
-		topocentric_frame_rotation.interpolate(ctx.alpha)
+		icrf_to_eus_translation.interpolate(ctx.alpha),
+		icrf_to_eus_rotation.interpolate(ctx.alpha)
 	};
 	
-	// Get topocentric space direction to sun
+	// Get EUS direction to sun
 	float3 sun_position = sun_position_tween.interpolate(ctx.alpha);
 	float3 sun_direction = math::normalize(sun_position);
 	
 	// Interpolate sun color
-	float3 sun_color_outer = sun_color_outer_tween.interpolate(ctx.alpha);
-	float3 sun_color_inner = sun_color_inner_tween.interpolate(ctx.alpha);
+	float3 sun_illuminance_outer = sun_illuminance_outer_tween.interpolate(ctx.alpha);
+	float3 sun_illuminance_inner = sun_illuminance_inner_tween.interpolate(ctx.alpha);
 	
 	// Draw atmosphere
 	if (sky_model)
@@ -149,8 +149,8 @@ void sky_pass::render(const render::context& ctx, render::queue& queue) const
 			sun_angular_radius_input->upload(sun_angular_radius);
 			
 		// Pre-exposure sun color
-		if (sun_color_input)
-			sun_color_input->upload(sun_color_outer * ctx.exposure);
+		if (sun_illuminance_input)
+			sun_illuminance_input->upload(sun_illuminance_outer * ctx.exposure);
 		
 		if (scale_height_rm_input)
 			scale_height_rm_input->upload(scale_height_rm);
@@ -177,8 +177,8 @@ void sky_pass::render(const render::context& ctx, render::queue& queue) const
 			cloud_model_view_projection_input->upload(model_view_projection);
 		if (cloud_sun_direction_input)
 			cloud_sun_direction_input->upload(sun_direction);
-		if (cloud_sun_color_input)
-			cloud_sun_color_input->upload(sun_color_inner);
+		if (cloud_sun_illuminance_input)
+			cloud_sun_illuminance_input->upload(sun_illuminance_inner);
 		if (cloud_camera_position_input)
 			cloud_camera_position_input->upload(ctx.camera_transform.translation);
 		if (cloud_camera_exposure_input)
@@ -198,7 +198,7 @@ void sky_pass::render(const render::context& ctx, render::queue& queue) const
 	{
 		float star_distance = (clip_near + clip_far) * 0.5f;
 		
-		model = math::resize<4, 4>(math::matrix_cast<float>(topocentric_frame.rotation));
+		model = math::resize<4, 4>(math::matrix_cast<float>(icrf_to_eus.r));
 		model = math::scale(model, {star_distance, star_distance, star_distance});
 		
 		model_view = view * model;
@@ -218,13 +218,12 @@ void sky_pass::render(const render::context& ctx, render::queue& queue) const
 		rasterizer->draw_arrays(*stars_model_vao, stars_model_drawing_mode, stars_model_start_index, stars_model_index_count);
 	}
 	
-	// Draw moon model
 	/*
-	float3 moon_position = {0, 0, 0};
+	// Draw moon model
+	float3 moon_position = moon_position_tween.interpolate(ctx.alpha);
+	float moon_angular_radius = math::radians(2.0f);
 	if (moon_position.y >= -moon_angular_radius)
 	{
-
-		
 		float moon_distance = (clip_near + clip_far) * 0.5f;		
 		float moon_radius = moon_angular_radius * moon_distance;
 		
@@ -284,7 +283,7 @@ void sky_pass::set_sky_model(const model* model)
 
 				observer_altitude_input = sky_shader_program->get_input("observer_altitude");
 				sun_direction_input = sky_shader_program->get_input("sun_direction");
-				sun_color_input = sky_shader_program->get_input("sun_color");
+				sun_illuminance_input = sky_shader_program->get_input("sun_illuminance");
 				sun_angular_radius_input = sky_shader_program->get_input("sun_angular_radius");
 				scale_height_rm_input = sky_shader_program->get_input("scale_height_rm");
 				rayleigh_scattering_input = sky_shader_program->get_input("rayleigh_scattering");
@@ -397,7 +396,7 @@ void sky_pass::set_clouds_model(const model* model)
 			{
 				cloud_model_view_projection_input = cloud_shader_program->get_input("model_view_projection");
 				cloud_sun_direction_input = cloud_shader_program->get_input("sun_direction");
-				cloud_sun_color_input = cloud_shader_program->get_input("sun_color");
+				cloud_sun_illuminance_input = cloud_shader_program->get_input("sun_illuminance");
 				cloud_camera_position_input = cloud_shader_program->get_input("camera.position");
 				cloud_camera_exposure_input = cloud_shader_program->get_input("camera.exposure");
 			}
@@ -413,16 +412,17 @@ void sky_pass::update_tweens()
 {
 	observer_altitude_tween.update();
 	sun_position_tween.update();
-	sun_color_outer_tween.update();
-	sun_color_inner_tween.update();
-	topocentric_frame_translation.update();
-	topocentric_frame_rotation.update();
+	sun_illuminance_outer_tween.update();
+	sun_illuminance_inner_tween.update();
+	icrf_to_eus_translation.update();
+	icrf_to_eus_rotation.update();
+	moon_position_tween.update();
 }
 
-void sky_pass::set_topocentric_frame(const physics::frame<float>& frame)
+void sky_pass::set_icrf_to_eus(const math::transformation::se3<float>& transformation)
 {
-	topocentric_frame_translation[1] = frame.translation;
-	topocentric_frame_rotation[1] = frame.rotation;
+	icrf_to_eus_translation[1] = transformation.t;
+	icrf_to_eus_rotation[1] = transformation.r;
 }
 
 void sky_pass::set_sun_position(const float3& position)
@@ -430,10 +430,10 @@ void sky_pass::set_sun_position(const float3& position)
 	sun_position_tween[1] = position;
 }
 
-void sky_pass::set_sun_color(const float3& color_outer, const float3& color_inner)
+void sky_pass::set_sun_illuminance(const float3& illuminance_outer, const float3& illuminance_inner)
 {
-	sun_color_outer_tween[1] = color_outer;
-	sun_color_inner_tween[1] = color_inner;
+	sun_illuminance_outer_tween[1] = illuminance_outer;
+	sun_illuminance_inner_tween[1] = illuminance_inner;
 }
 
 void sky_pass::set_sun_angular_radius(float radius)
@@ -469,6 +469,11 @@ void sky_pass::set_atmosphere_radii(float inner, float outer)
 	atmosphere_radii.z = outer * outer;
 }
 
+void sky_pass::set_moon_position(const float3& position)
+{
+	moon_position_tween[1] = position;
+}
+
 void sky_pass::handle_event(const mouse_moved_event& event)
 {
 	mouse_position = {static_cast<float>(event.x), static_cast<float>(event.y)};

src/render/passes/sky-pass.hpp

@@ -32,7 +32,7 @@
 #include "gl/vertex-array.hpp"
 #include "gl/texture-2d.hpp"
 #include "gl/drawing-mode.hpp"
-#include "physics/frame.hpp"
+#include "math/se3.hpp"
 #include "scene/object.hpp"
 
 class resource_manager;
@@ -60,16 +60,18 @@ public:
 	void set_stars_model(const model* model);
 	void set_clouds_model(const model* model);
 	
-	void set_topocentric_frame(const physics::frame<float>& frame);
+	void set_icrf_to_eus(const math::transformation::se3<float>& transformation);
 	
 	void set_sun_position(const float3& position);
-	void set_sun_color(const float3& color_outer, const float3& color_inner);
+	void set_sun_illuminance(const float3& illuminance_outer, const float3& illuminance_inner);
 	void set_sun_angular_radius(float radius);
 	void set_observer_altitude(float altitude);
 	void set_scale_heights(float rayleigh, float mie);
 	void set_scattering_coefficients(const float3& r, const float3& m);
 	void set_mie_anisotropy(float g);
 	void set_atmosphere_radii(float inner, float outer);
+	
+	void set_moon_position(const float3& position);
 
 private:
 	virtual void handle_event(const mouse_moved_event& event);
@@ -82,7 +84,7 @@ private:
 	const gl::shader_input* exposure_input;
 	const gl::shader_input* observer_altitude_input;
 	const gl::shader_input* sun_direction_input;
-	const gl::shader_input* sun_color_input;
+	const gl::shader_input* sun_illuminance_input;
 	const gl::shader_input* sun_angular_radius_input;
 	const gl::shader_input* scale_height_rm_input;
 	const gl::shader_input* rayleigh_scattering_input;
@@ -132,7 +134,7 @@ private:
 	gl::shader_program* cloud_shader_program;
 	const gl::shader_input* cloud_model_view_projection_input;
 	const gl::shader_input* cloud_sun_direction_input;
-	const gl::shader_input* cloud_sun_color_input;
+	const gl::shader_input* cloud_sun_illuminance_input;
 	const gl::shader_input* cloud_camera_position_input;
 	const gl::shader_input* cloud_camera_exposure_input;
 
@@ -143,10 +145,13 @@ private:
 	
 	tween<float> observer_altitude_tween;
 	tween<float3> sun_position_tween;
-	tween<float3> sun_color_outer_tween;
-	tween<float3> sun_color_inner_tween;
-	tween<float3> topocentric_frame_translation;
-	tween<math::quaternion<float>> topocentric_frame_rotation;
+	tween<float3> sun_illuminance_outer_tween;
+	tween<float3> sun_illuminance_inner_tween;
+	tween<float3> icrf_to_eus_translation;
+	tween<math::quaternion<float>> icrf_to_eus_rotation;
+	
+	tween<float3> moon_position_tween;
+
 	
 	float sun_angular_radius;
 	float2 scale_height_rm;

src/resources/entity-archetype-loader.cpp

@@ -31,6 +31,7 @@
 #include "entity/components/celestial-body.hpp"
 #include "entity/archetype.hpp"
 #include "entity/ebt.hpp"
+#include "physics/orbit/elements.hpp"
 #include "resources/json.hpp"
 #include <stdexcept>
 
@@ -97,18 +98,24 @@ static bool load_component_celestial_body(entity::archetype& archetype, const js
 {
 	entity::component::celestial_body component;
 	component.radius = 0.0;
-	component.axial_tilt = 0.0;
-	component.axial_rotation = 0.0;
-	component.angular_frequency = 0.0;
+	component.mass = 0.0;
+	component.pole_ra = 0.0;
+	component.pole_dec = 0.0;
+	component.prime_meridian = 0.0;
+	component.rotation_period = 0.0;
 	
 	if (element.contains("radius"))
 		component.radius = element["radius"].get<double>();
-	if (element.contains("axial_tilt"))
-		component.axial_tilt = math::radians(element["axial_tilt"].get<double>());
-	if (element.contains("axial_rotation"))
-		component.axial_rotation = math::radians(element["axial_rotation"].get<double>());
-	if (element.contains("angular_frequency"))
-		component.angular_frequency = math::radians(element["angular_frequency"].get<double>());
+	if (element.contains("mass"))
+		component.mass = element["mass"].get<double>();
+	if (element.contains("pole_ra"))
+		component.pole_ra = math::radians(element["pole_ra"].get<double>());
+	if (element.contains("pole_dec"))
+		component.pole_dec = math::radians(element["pole_dec"].get<double>());
+	if (element.contains("prime_meridian"))
+		component.prime_meridian = math::radians(element["prime_meridian"].get<double>());
+	if (element.contains("rotation_period"))
+		component.rotation_period = element["rotation_period"].get<double>();
 	
 	archetype.set<entity::component::celestial_body>(component);
 
@@ -151,26 +158,30 @@ static bool load_component_orbit(entity::archetype& archetype, const json& eleme
 {
 	entity::component::orbit component;
 	
-	component.elements.e = 0.0;
+	component.parent = entt::null;
+	component.elements.ec = 0.0;
 	component.elements.a = 0.0;
-	component.elements.i = 0.0;
-	component.elements.raan = 0.0;
+	component.elements.in = 0.0;
+	component.elements.om = 0.0;
 	component.elements.w = 0.0;
-	component.elements.ta = 0.0;
+	component.elements.ma = 0.0;
+	component.mean_motion = 0.0;
 	
-	if (element.contains("e"))
-		component.elements.e = element["e"].get<double>();
+	if (element.contains("ec"))
+		component.elements.ec = element["ec"].get<double>();
 	if (element.contains("a"))
 		component.elements.a = element["a"].get<double>();
-	if (element.contains("i"))
-		component.elements.i = math::radians(element["i"].get<double>());
-	if (element.contains("raan"))
-		component.elements.raan = math::radians(element["raan"].get<double>());
+	if (element.contains("in"))
+		component.elements.in = math::radians(element["in"].get<double>());
+	if (element.contains("om"))
+		component.elements.om = math::radians(element["om"].get<double>());
 	if (element.contains("w"))
 		component.elements.w = math::radians(element["w"].get<double>());
-	if (element.contains("ta"))
-		component.elements.ta = math::radians(element["ta"].get<double>());
-
+	if (element.contains("ma"))
+		component.elements.ma = math::radians(element["ma"].get<double>());
+	if (element.contains("n"))
+		component.mean_motion = math::radians(element["n"].get<double>());
+	
 	archetype.set<entity::component::orbit>(component);
 
 	return true;