From e777f5c63cc86f7e66182b2418faad8d0d002715 Mon Sep 17 00:00:00 2001 From: "C. J. Howard" Date: Fri, 9 Oct 2020 11:19:24 -0700 Subject: [PATCH] Add functions for solving Kepler's equation, add solar system class --- CMakeLists.txt | 1 + src/game/astronomy/celestial-mechanics.cpp | 70 +++++++++++++- src/game/astronomy/celestial-mechanics.hpp | 22 ++--- src/game/components/orbit-component.hpp | 51 +++++++++++ src/game/states/play-state.cpp | 2 +- src/game/systems/solar-system.cpp | 101 +++++++++++++++++++++ src/game/systems/solar-system.hpp | 54 +++++++++++ src/game/systems/weather-system.cpp | 11 ++- 8 files changed, 294 insertions(+), 18 deletions(-) create mode 100644 src/game/components/orbit-component.hpp create mode 100644 src/game/systems/solar-system.cpp create mode 100644 src/game/systems/solar-system.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index b80be9e..18ea3f6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,6 +4,7 @@ option(VERSION_STRING "Project version string" "0.0.0") project(antkeeper VERSION ${VERSION_STRING} LANGUAGES CXX) + # Find dependency packages find_package(dr_wav REQUIRED CONFIG) find_package(stb REQUIRED CONFIG) diff --git a/src/game/astronomy/celestial-mechanics.cpp b/src/game/astronomy/celestial-mechanics.cpp index ac29672..7a63bf7 100644 --- a/src/game/astronomy/celestial-mechanics.cpp +++ b/src/game/astronomy/celestial-mechanics.cpp @@ -70,7 +70,7 @@ double3 approx_moon_ecliptic(double jd) + 0.0007 * std::sin(m1 - m) - 0.0006 * std::sin(d) - 0.0005 * std::sin(m + m1); - + const double latitude = 0.0895 * sin(f) + 0.0049 * std::sin(m1 + f) + 0.0048 * std::sin(m1 - f) @@ -129,4 +129,72 @@ double3x3 approx_moon_ecliptic_rotation(double jd) return rz0 * rx * rz1; } +double3 solve_kepler(const kepler_orbit& orbit, double t) +{ + // Orbital period (Kepler's third law) + double pr = std::sqrt(orbit.a * orbit.a * orbit.a); + + // 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 + { + 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 + }; +} + +double3 solve_kepler(double a, double ec, double w, double ma, double i, double om) +{ + // Semi-minor axis + double b = std::sqrt(1.0 - ec * ec); + + // Eccentric anomaly + double ea = ma + ec * std::sin(ma) * (1.0 + ec * std::cos(ma)); + + // Find radial distance (r) and true anomaly (v) + double x = a * (std::cos(ea) - 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 + { + 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) + }; +} + } // namespace ast diff --git a/src/game/astronomy/celestial-mechanics.hpp b/src/game/astronomy/celestial-mechanics.hpp index a2f938f..f61eaef 100644 --- a/src/game/astronomy/celestial-mechanics.hpp +++ b/src/game/astronomy/celestial-mechanics.hpp @@ -57,23 +57,19 @@ double3 approx_moon_ecliptic(double jd); */ double3x3 approx_moon_ecliptic_rotation(double jd); -struct orbital_elements +struct kepler_orbit { - double a; ///< Semi-major axis (km) - double e; ///< Eccentricity - double w; ///< Argument of periapsis (radians) - double m; ///< Mean anomaly (radians) - double i; ///< Inclination to the ecliptic (radians) - double n; ///< Longitude of the ascending node (radians) + 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) }; -//constexpr orbital_elements j2k_orbit_moon = {384400.0, 0.0554, 318.15, 135.275.16, 125.08}; +double3 solve_kepler(const kepler_orbit& orbit, double t); -struct orbital_state -{ - double3 r; ///< Cartesian position - double3 v; ///< Cartesian velocity -}; +double3 solve_kepler(double a, double ec, double w, double ma, double i, double om); } // namespace ast diff --git a/src/game/components/orbit-component.hpp b/src/game/components/orbit-component.hpp new file mode 100644 index 0000000..6468905 --- /dev/null +++ b/src/game/components/orbit-component.hpp @@ -0,0 +1,51 @@ +/* + * Copyright (C) 2020 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 . + */ + +#ifndef ANTKEEPER_ECS_ORBIT_COMPONENT_HPP +#define ANTKEEPER_ECS_ORBIT_COMPONENT_HPP + +#include "game/astronomy/celestial-mechanics.hpp" + +namespace ecs { + +struct orbit_component +{ + double a; ///< Semi-major axis, a + double d_a; + + double ec; ///< Eccentricity, e + double d_ec; + + double w; ///< Argument of periapsis, w (radians) + double d_w; + + double ma; ///< Mean anomaly, M (radians) + double d_ma; + + double i; ///< Inclination, i (radians) + double d_i; + + double om; ///< Longitude of the ascending node, OMEGA (radians) + double d_om; +}; + +} // namespace ecs + +#endif // ANTKEEPER_ECS_ORBIT_COMPONENT_HPP + diff --git a/src/game/states/play-state.cpp b/src/game/states/play-state.cpp index c3c52e0..efc3349 100644 --- a/src/game/states/play-state.cpp +++ b/src/game/states/play-state.cpp @@ -81,7 +81,7 @@ void play_state_enter(game_context* ctx) sky_pass->set_moon_model(ctx->resource_manager->load("moon.mdl")); ctx->weather_system->set_location(ctx->biome->location[0], ctx->biome->location[1], ctx->biome->location[2]); - ctx->weather_system->set_time(2017, 8, 21, 5, 0, 0.0, -7.0); + ctx->weather_system->set_time(2017, 6, 1, 5, 0, 0.0, -7.0); ctx->weather_system->set_sky_palette(ctx->biome->sky_palette); ctx->weather_system->set_sun_palette(ctx->biome->sun_palette); ctx->weather_system->set_ambient_palette(ctx->biome->ambient_palette); diff --git a/src/game/systems/solar-system.cpp b/src/game/systems/solar-system.cpp new file mode 100644 index 0000000..5d303cd --- /dev/null +++ b/src/game/systems/solar-system.cpp @@ -0,0 +1,101 @@ +/* + * Copyright (C) 2020 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 . + */ + +#include "game/systems/solar-system.hpp" +#include "game/astronomy/celestial-coordinates.hpp" +#include "game/astronomy/celestial-mechanics.hpp" +#include "game/astronomy/celestial-time.hpp" +#include "game/components/orbit-component.hpp" +#include "game/components/transform-component.hpp" + +using namespace ecs; + +static constexpr double seconds_per_day = 24.0 * 60.0 * 60.0; + +solar_system::solar_system(entt::registry& registry): + entity_system(registry), + julian_date(0.0), + time_scale(1.0), + latitude(0.0), + longitude(0.0), + altitude(0.0) +{} + +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); + + // 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; + double ma = orbit.ma; + double i = orbit.i; + double om = orbit.om; + + double3 ecliptic = ast::solve_kepler(a, ec, w, ma, i, om); + double3 horizontal = ecliptic_to_horizontal * ecliptic; + + // Subtract Earth's radius (in AU), for positon of observer + horizontal.z -= 4.25875e-5; + + // Transform into local right-handed coordinates + double3 translation = ast::horizontal_to_right_handed * horizontal; + double3x3 rotation = ast::horizontal_to_right_handed * ecliptic_to_horizontal; + + transform.local.translation = math::type_cast(translation); + transform.local.rotation = math::type_cast(math::quaternion_cast(rotation)); + }); +} + +void solar_system::set_julian_date(double jd) +{ + julian_date = jd; + + // Recalculate obliquity of the ecliptic + ecl = ast::approx_ecliptic_obliquity(julian_date); + + // Recalculate LMST + lmst = ast::jd_to_lmst(julian_date, longitude); + + // Recalculate ecliptic to horizontal transformation matrix + ecliptic_to_horizontal = ast::ecliptic_to_horizontal(ecl, latitude, lmst); +} + +void solar_system::set_time_scale(double scale) +{ + time_scale = scale; +} + +void solar_system::set_observer_location(double latitude, double longitude, double altitude) +{ + this->latitude = latitude; + this->longitude = longitude; + this->altitude = altitude; + + // Recalculate LMST + lmst = ast::jd_to_lmst(julian_date, longitude); + + // Recalculate ecliptic to horizontal transformation matrix + ecliptic_to_horizontal = ast::ecliptic_to_horizontal(ecl, latitude, lmst); +} diff --git a/src/game/systems/solar-system.hpp b/src/game/systems/solar-system.hpp new file mode 100644 index 0000000..1b407e4 --- /dev/null +++ b/src/game/systems/solar-system.hpp @@ -0,0 +1,54 @@ +/* + * Copyright (C) 2020 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 . + */ + +#ifndef ANTKEEPER_SOLAR_SYSTEM_HPP +#define ANTKEEPER_SOLAR_SYSTEM_HPP + +#include "entity-system.hpp" +#include "utility/fundamental-types.hpp" + +class solar_system: + public entity_system +{ +public: + solar_system(entt::registry& registry); + + virtual void update(double t, double dt); + + void set_julian_date(double jd); + + void set_time_scale(double scale); + + void set_observer_location(double latitude, double longitude, double altitude); + +private: + double julian_date; + double time_scale; + double latitude; + double longitude; + double altitude; + + double lmst; + double ecl; + double3x3 ecliptic_to_horizontal; + + double3 sun_position; +}; + +#endif // ANTKEEPER_SOLAR_SYSTEM_HPP diff --git a/src/game/systems/weather-system.cpp b/src/game/systems/weather-system.cpp index f49601c..59df97b 100644 --- a/src/game/systems/weather-system.cpp +++ b/src/game/systems/weather-system.cpp @@ -83,11 +83,16 @@ void weather_system::update(double t, double dt) double3 moon_ecliptic = ast::approx_moon_ecliptic(jd); double3 moon_horizontal = ecliptic_to_horizontal * moon_ecliptic; moon_horizontal.z -= 1.0; // Subtract one earth radius, for position of observer + double3 moon_spherical = ast::rectangular_to_spherical(moon_horizontal); double3 moon_positiond = ast::horizontal_to_right_handed * moon_horizontal; - float2 moon_az_el = {static_cast(moon_spherical.z) - math::pi, static_cast(moon_spherical.y)}; - float3 moon_position = math::normalize(float3{static_cast(moon_positiond.x), static_cast(moon_positiond.y), static_cast(moon_positiond.z)}); - + float2 moon_az_el = {static_cast(moon_spherical.z) - math::pi, static_cast(moon_spherical.y)}; + float3 moon_position = math::normalize(math::type_cast(moon_positiond)); + + //double3 moon_sphericald = ast::rectangular_to_spherical(moon_positiond); + //std::cout << "old azel: " << math::degrees(moon_az_el.x) << ", " << math::degrees(moon_az_el.y) << std::endl; + //std::cout << "new azel: " << math::degrees(moon_sphericald.z) << ", " << math::degrees(moon_sphericald.y) << std::endl; + double3x3 moon_rotation_matrix = ast::horizontal_to_right_handed * ecliptic_to_horizontal; math::quaternion moon_rotationd = math::normalize(math::quaternion_cast(moon_rotation_matrix) * math::angle_axis(math::half_pi, double3{0, 1, 0}) * math::angle_axis(-math::half_pi, double3{0, 0, -1})); math::quaternion moon_rotation =