From e1d5b6b3cf2ee854008d613eaea7fc95353b0f61 Mon Sep 17 00:00:00 2001 From: "C. J. Howard" Date: Sun, 18 Sep 2022 07:22:59 +0800 Subject: [PATCH] Add ephemeris loader and change orbit system to be ephemeris-based --- src/astro/astro.hpp | 31 --- src/astro/illuminance.hpp | 58 ----- src/entity/commands.cpp | 4 +- src/entity/components/atmosphere.hpp | 3 + src/entity/components/blackbody.hpp | 2 +- src/entity/components/celestial-body.hpp | 32 +-- src/entity/components/orbit.hpp | 27 +- src/entity/systems/astronomy.cpp | 162 +++++++----- src/entity/systems/astronomy.hpp | 14 +- src/entity/systems/orbit.cpp | 82 ++---- src/entity/systems/orbit.hpp | 25 +- src/entity/systems/terrain.cpp | 2 +- src/game/ant/morphogenesis.cpp | 2 +- src/game/state/boot.cpp | 2 +- src/game/state/brood.cpp | 6 +- src/game/state/forage.cpp | 4 +- src/game/state/nuptial-flight.cpp | 87 ++++--- src/game/world.cpp | 31 ++- src/game/world.hpp | 3 + src/geom/view-frustum.hpp | 2 +- src/math/constants.hpp | 66 +---- src/math/interpolation.hpp | 9 +- src/math/matrix-type.hpp | 26 +- src/math/polynomial.hpp | 11 +- src/math/quaternion-type.hpp | 9 +- src/math/transform-type.hpp | 12 +- src/physics/light/exposure.hpp | 114 ++++++++ src/physics/light/light.hpp | 1 + src/physics/light/phase.hpp | 2 +- .../light/vmag.hpp} | 50 +++- .../orbit/ephemeris.hpp} | 25 +- src/physics/orbit/frame.hpp | 10 +- src/physics/orbit/trajectory.hpp | 84 ++++++ src/render/passes/ground-pass.cpp | 2 +- src/render/passes/shadow-map-pass.cpp | 8 +- src/render/passes/sky-pass.cpp | 43 +-- src/render/passes/sky-pass.hpp | 8 +- src/resources/entity-archetype-loader.cpp | 69 ++--- src/resources/ephemeris-loader.cpp | 244 ++++++++++++++++++ src/scene/camera.cpp | 11 +- src/scene/camera.hpp | 9 - src/scene/object.cpp | 2 +- src/utility/bit-math.hpp | 55 ++++ 43 files changed, 950 insertions(+), 499 deletions(-) delete mode 100644 src/astro/astro.hpp delete mode 100644 src/astro/illuminance.hpp create mode 100644 src/physics/light/exposure.hpp rename src/{astro/illuminance.cpp => physics/light/vmag.hpp} (50%) rename src/{astro/constants.hpp => physics/orbit/ephemeris.hpp} (66%) create mode 100644 src/physics/orbit/trajectory.hpp create mode 100644 src/resources/ephemeris-loader.cpp diff --git a/src/astro/astro.hpp b/src/astro/astro.hpp deleted file mode 100644 index 8e6115d..0000000 --- a/src/astro/astro.hpp +++ /dev/null @@ -1,31 +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 . - */ - -#ifndef ANTKEEPER_ASTRO_HPP -#define ANTKEEPER_ASTRO_HPP - -/// Astronomy/astrophysics -namespace astro {} - -#include "apparent-size.hpp" -#include "coordinates.hpp" -#include "constants.hpp" -#include "illuminance.hpp" - -#endif // ANTKEEPER_ASTRO_HPP diff --git a/src/astro/illuminance.hpp b/src/astro/illuminance.hpp deleted file mode 100644 index 257eef1..0000000 --- a/src/astro/illuminance.hpp +++ /dev/null @@ -1,58 +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 . - */ - -#ifndef ANTKEEPER_ASTRO_ILLUMINANCE_HPP -#define ANTKEEPER_ASTRO_ILLUMINANCE_HPP - -namespace astro -{ - -/** - * Converts apparent (visual) magnitude to a brightness factor relative to a 0th magnitude star. - * - * @param mv Illuminance in apparent magnitude. - * @return Relative brightness factor. - * - * @see https://en.wikipedia.org/wiki/Illuminance - */ -double vmag_to_brightness(double mv); - -/** - * Converts apparent (visual) magnitude to lux. - * - * @param mv Illuminance in apparent magnitude. - * @return Illuminance in lux. - * - * @see https://en.wikipedia.org/wiki/Illuminance - */ -double vmag_to_lux(double mv); - -/** - * Converts lux to apparent (visual) magnitude. - * - * @param ev Illuminance in lux. - * @return Illuminance in apparent magnitude. - * - * @see https://en.wikipedia.org/wiki/Illuminance - */ -double lux_to_vmag(double ev); - -} // namespace astro - -#endif // ANTKEEPER_ASTRO_ILLUMINANCE_HPP diff --git a/src/entity/commands.cpp b/src/entity/commands.cpp index 30a0670..dde26bd 100644 --- a/src/entity/commands.cpp +++ b/src/entity/commands.cpp @@ -143,7 +143,7 @@ math::transform get_local_transform(entity::registry& registry, entity::i return component.local; } - return math::identity_transform; + return math::transform::identity; } math::transform get_world_transform(entity::registry& registry, entity::id eid) @@ -154,7 +154,7 @@ math::transform get_world_transform(entity::registry& registry, entity::i return component.world; } - return math::identity_transform; + return math::transform::identity; } void parent(entity::registry& registry, entity::id child, entity::id parent) diff --git a/src/entity/components/atmosphere.hpp b/src/entity/components/atmosphere.hpp index 47994a8..7e9330d 100644 --- a/src/entity/components/atmosphere.hpp +++ b/src/entity/components/atmosphere.hpp @@ -54,6 +54,9 @@ struct atmosphere /// (Dependent) Mie scattering coefficients at sea level. double3 mie_scattering; + + /// Airglow, in lux. + double airglow; }; } // namespace component diff --git a/src/entity/components/blackbody.hpp b/src/entity/components/blackbody.hpp index 24d608c..45932fe 100644 --- a/src/entity/components/blackbody.hpp +++ b/src/entity/components/blackbody.hpp @@ -29,7 +29,7 @@ struct blackbody /// Effective temperature, in Kelvin. double temperature; - /// (Dependent) RGB luminance, in lumens. + /// (Dependent) RGB spectral luminance, in cd/m^2. double3 luminance; }; diff --git a/src/entity/components/celestial-body.hpp b/src/entity/components/celestial-body.hpp index fe5b056..16460cd 100644 --- a/src/entity/components/celestial-body.hpp +++ b/src/entity/components/celestial-body.hpp @@ -32,35 +32,29 @@ struct celestial_body /// Mass of the body, in kilograms. double mass; - /// Right ascension of the body's north pole at epoch, in radians. - double pole_ra; + /// Polynomial coefficients, in descending order of degree, of the right ascension of the body's north pole, in radians, w.r.t. Julian centuries (36525 days) from epoch. + std::vector pole_ra; - /// Declination of the body's north pole at epoch, in radians. - double pole_dec; + /// Polynomial coefficients, in descending order of degree, of the declination of the body's north pole, in radians, w.r.t. Julian centuries (36525 days) from epoch. + std::vector pole_dec; - /// Location of the prime meridian at epoch, as a rotation about the north pole, in radians. - double prime_meridian; + /// Polynomial coefficients, in descending order of degree, of the rotation state of the body's prime meridian, in radians, w.r.t. days from epoch. + std::vector prime_meridian; /* - /// Quadratic e coefficients for the right ascension of the body's north pole, in radians. Right ascension is calculated as `x + y * T + z * T^2`, where `T` is the Julian centuries (36525 days) from epoch. - double3 pole_ra; + /// Polynomial coefficients, in descending order of degree, of the nutation and precession angles, in radians. Angles are calculated as `x + y * d`, where `d` is the days from epoch. + std::vector> nutation_precession_angles; - /// Quadratic coefficients for the declination of the body's north pole, in radians. Declination is calculated as `x + y * T + z * T^2`, where `T` is the Julian centuries (36525 days) from epoch. - double3 pole_ra; - - /// Quadratic coefficients for the rotation state of the body's prime meridian, in radians. Prime meridian rotation is calculated as `x + y * d + z * d^2`, where `d` is the days from epoch. - double3 prime_meridian; - - /// Linear coefficients of the nutation and precession angles, in radians. Angles are calculated as `x + y * d`, where `d` is the days from epoch. - std::vector nutation_precession_angles; + /// Nutation and precession amplitudes of the right ascension of the body's north pole, in radians. std::vector nutation_precession_ra; + + /// Nutation and precession amplitudes of the declination of the body's north pole, in radians. std::vector nutation_precession_dec; + + /// Nutation and precession amplitudes of the rotation state of the body's prime meridian, in radians. std::vector nutation_precession_pm; */ - /// Sidereal rotation period, in rotations per day. - double rotation_period; - /// Geometric albedo double albedo; }; diff --git a/src/entity/components/orbit.hpp b/src/entity/components/orbit.hpp index e19af70..f31f379 100644 --- a/src/entity/components/orbit.hpp +++ b/src/entity/components/orbit.hpp @@ -21,35 +21,24 @@ #define ANTKEEPER_ENTITY_COMPONENT_ORBIT_HPP #include "entity/id.hpp" -#include "physics/orbit/elements.hpp" -#include "physics/orbit/state.hpp" -#include "math/se3.hpp" +#include "math/vector-type.hpp" namespace entity { namespace component { struct orbit { - /// Parent body. + /// Entity ID of the parent orbit. entity::id parent; - /// Keplerian orbital elements at epoch. - physics::orbit::elements elements; + /// Index of the orbit in the ephemeris. + int ephemeris_index; - /// Orbital semi-minor axis (b) at epoch. - double semiminor_axis; + /// Orbit scale, for two-body orbits with one ephemeris item. + double scale; - /// 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 pqw_to_bci; - - /// Orbital Cartesian position of the body in the BCI frame of the parent body. - math::vector3 bci_position; - - /// Orbital Cartesian position of the body in the ICRF frame. - math::vector3 icrf_position; + /// Cartesian position of the orbit, w.r.t. the ICRF frame. + math::vector3 position; }; } // namespace component diff --git a/src/entity/systems/astronomy.cpp b/src/entity/systems/astronomy.cpp index 402f6c5..3bd6c15 100644 --- a/src/entity/systems/astronomy.cpp +++ b/src/entity/systems/astronomy.cpp @@ -33,8 +33,8 @@ #include "physics/atmosphere.hpp" #include "geom/cartesian.hpp" #include "astro/apparent-size.hpp" -#include "astro/illuminance.hpp" #include "geom/solid-angle.hpp" +#include "math/polynomial.hpp" #include namespace entity { @@ -57,7 +57,7 @@ math::vector3 transmittance(T depth_r, T depth_m, T depth_o, const math::vect astronomy::astronomy(entity::registry& registry): updatable(registry), - universal_time(0.0), + time(0.0), time_scale(1.0), reference_entity(entt::null), observer_location{0, 0, 0}, @@ -65,7 +65,8 @@ astronomy::astronomy(entity::registry& registry): sky_light(nullptr), moon_light(nullptr), camera(nullptr), - sky_pass(nullptr) + sky_pass(nullptr), + exposure_offset(0.0) { // Construct transformation which transforms coordinates from ENU to EUS enu_to_eus = math::transformation::se3 @@ -81,9 +82,10 @@ astronomy::astronomy(entity::registry& registry): void astronomy::update(double t, double dt) { double total_illuminance = 0.0; + double sky_light_illuminance = 0.0; // Add scaled timestep to current time - set_universal_time(universal_time + dt * time_scale); + set_time(time + dt * time_scale); // Abort if no reference body if (reference_entity == entt::null) @@ -95,16 +97,24 @@ void astronomy::update(double t, double dt) const entity::component::orbit& reference_orbit = registry.get(reference_entity); const entity::component::celestial_body& reference_body = registry.get(reference_entity); - math::transformation::se3 icrf_to_bci{{-reference_orbit.icrf_position}, math::identity_quaternion}; + math::transformation::se3 icrf_to_bci{{-reference_orbit.position}, math::quaternion::identity}; + + const double days_from_epoch = time; + const double centuries_from_epoch = days_from_epoch / 36525.0; + + // Evaluate reference body orientation polynomials + const double reference_body_pole_ra = math::polynomial::horner(reference_body.pole_ra.begin(), reference_body.pole_ra.end(), centuries_from_epoch); + const double reference_body_pole_dec = math::polynomial::horner(reference_body.pole_dec.begin(), reference_body.pole_dec.end(), centuries_from_epoch); + const double reference_body_prime_meridian = math::polynomial::horner(reference_body.prime_meridian.begin(), reference_body.prime_meridian.end(), days_from_epoch); // Construct transformation from the ICRF frame to the reference body BCBF frame icrf_to_bcbf = physics::orbit::frame::bci::to_bcbf ( - reference_body.pole_ra, - reference_body.pole_dec, - reference_body.prime_meridian + (math::two_pi / reference_body.rotation_period) * universal_time + reference_body_pole_ra, + reference_body_pole_dec, + reference_body_prime_meridian ); - icrf_to_bcbf.t = icrf_to_bcbf.r * -reference_orbit.icrf_position; + icrf_to_bcbf.t = icrf_to_bcbf.r * -reference_orbit.position; icrf_to_enu = icrf_to_bcbf * bcbf_to_enu; icrf_to_eus = icrf_to_enu * enu_to_eus; @@ -118,14 +128,19 @@ void astronomy::update(double t, double dt) return; // Transform orbital Cartesian position (r) from the ICRF frame to the EUS frame - const double3 r_eus = icrf_to_eus * orbit.icrf_position; + const double3 r_eus = icrf_to_eus * orbit.position; + + // Evaluate body orientation polynomials + const double body_pole_ra = math::polynomial::horner(body.pole_ra.begin(), body.pole_ra.end(), centuries_from_epoch); + const double body_pole_dec = math::polynomial::horner(body.pole_dec.begin(), body.pole_dec.end(), centuries_from_epoch); + const double body_prime_meridian = math::polynomial::horner(body.prime_meridian.begin(), body.prime_meridian.end(), days_from_epoch); // Determine body orientation in the ICRF frame math::quaternion rotation_icrf = physics::orbit::frame::bcbf::to_bci ( - body.pole_ra, - body.pole_dec, - body.prime_meridian + (math::two_pi / body.rotation_period) * this->universal_time + body_pole_ra, + body_pole_dec, + body_prime_meridian ).r; // Transform body orientation from the ICRF frame to the EUS frame. @@ -139,11 +154,11 @@ void astronomy::update(double t, double dt) transform.local.scale = {1.0f, 1.0f, 1.0f}; } - + /* if (orbit.parent != entt::null) { // RA-DEC - const double3 r_bci = icrf_to_bci * orbit.icrf_position; + const double3 r_bci = icrf_to_bci * orbit.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; @@ -151,17 +166,17 @@ void astronomy::update(double t, double dt) const double ra = math::degrees(r_bci_spherical.z); // AZ-EL - const double3 r_enu = icrf_to_enu * orbit.icrf_position; + const double3 r_enu = icrf_to_enu * orbit.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; 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; + std::cout << "t: " << this->time << "; ra: " << ra << "; dec: " << dec << std::endl; + std::cout << "t: " << this->time << "; az: " << az << "; el: " << el << std::endl; } - + */ }); // Update blackbody lighting @@ -173,14 +188,21 @@ void astronomy::update(double t, double dt) double3 blackbody_up_icrf = {0, 0, 1}; // 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_position_eus = icrf_to_eus * orbit.position; + double3 blackbody_position_enu = icrf_to_enu * orbit.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_eus) - body.radius; + // Calculate blackbody solid angle + const double blackbody_angular_radius = astro::angular_radius(body.radius, blackbody_distance); + const double blackbody_solid_angle = geom::solid_angle::cone(blackbody_angular_radius); + + // Calculate blackbody illuminance + const double3 blackbody_illuminance = blackbody.luminance * blackbody_solid_angle; + // Init atmospheric transmittance double3 atmospheric_transmittance = {1.0, 1.0, 1.0}; @@ -211,11 +233,20 @@ void astronomy::update(double t, double dt) atmospheric_transmittance = transmittance(optical_depth_r, optical_depth_k, optical_depth_o, reference_atmosphere.rayleigh_scattering, reference_atmosphere.mie_scattering); } + + // Add airglow to sky light illuminance + sky_light_illuminance += reference_atmosphere.airglow; } + // Blackbody illuminance transmitted through the atmosphere + const double3 transmitted_blackbody_illuminance = blackbody_illuminance * atmospheric_transmittance; + + // Add atmosphere-transmitted blackbody illuminance to total illuminance + total_illuminance += (transmitted_blackbody_illuminance.x + transmitted_blackbody_illuminance.y + transmitted_blackbody_illuminance.z) / 3.0; + + // Update sun light if (sun_light != nullptr) { - // Update blackbody light transform sun_light->set_translation(math::normalize(math::type_cast(blackbody_position_eus))); sun_light->set_rotation ( @@ -226,56 +257,49 @@ void astronomy::update(double t, double dt) ) ); - // Calculate blackbody solid angle - const double angular_radius = astro::angular_radius(body.radius, blackbody_distance); - const double solid_angle = geom::solid_angle::cone(angular_radius); - - const double3 blackbody_illuminance = blackbody.luminance * solid_angle; + sun_light->set_color(math::type_cast(transmitted_blackbody_illuminance)); + } + + // Update sky light + if (sky_light != nullptr) + { + // Calculate sky illuminance + double3 blackbody_position_enu_spherical = physics::orbit::frame::enu::spherical(blackbody_position_enu); + const double sky_illuminance = 25000.0 * std::max(0.0, std::sin(blackbody_position_enu_spherical.y)); - // Sun illuminance at the outer atmosphere - float3 sun_illuminance_outer = math::type_cast(blackbody_illuminance); + // Add sky illuminance to sky light illuminance + sky_light_illuminance += sky_illuminance; - // Sun illuminance at sea level - float3 sun_illuminance_inner = math::type_cast(blackbody_illuminance * atmospheric_transmittance); + // Add starlight illuminance to sky light illuminance + const double starlight_illuminance = 0.0002; + sky_light_illuminance += starlight_illuminance; - // Update blackbody light color and intensity - sun_light->set_color(sun_illuminance_inner); - sun_light->set_intensity(1.0f); + // Add sky light illuminance to total illuminance + total_illuminance += sky_light_illuminance; - total_illuminance += (sun_illuminance_inner.x + sun_illuminance_inner.y + sun_illuminance_inner.z) / 3.0; + //std::cout << "sky light illum: " << sky_light_illuminance << std::endl; - // Upload blackbody params to sky pass - if (this->sky_pass) - { - this->sky_pass->set_sun_position(math::type_cast(blackbody_position_eus)); - this->sky_pass->set_sun_illuminance(sun_illuminance_outer, sun_illuminance_inner); - - double blackbody_angular_radius = astro::angular_radius(body.radius, blackbody_distance); - this->sky_pass->set_sun_angular_radius(static_cast(blackbody_angular_radius)); - } + // Update sky light + sky_light->set_color(float3{1.0f, 1.0f, 1.0f} * static_cast(sky_light_illuminance)); } - if (sky_light != nullptr) + // Upload blackbody params to sky pass + if (this->sky_pass) { - double3 blackbody_position_enu_spherical = physics::orbit::frame::enu::spherical(icrf_to_enu * orbit.icrf_position); - - const double starlight_illuminance = 0.0011; - const double sky_illuminance = 25000.0 * std::max(0.0, std::sin(blackbody_position_enu_spherical.y)); - const double sky_light_illuminance = sky_illuminance + starlight_illuminance; - total_illuminance += sky_light_illuminance; - - sky_light->set_color({1.0f, 1.0f, 1.0f}); - sky_light->set_intensity(static_cast(sky_illuminance)); + this->sky_pass->set_sun_position(math::type_cast(blackbody_position_eus)); + this->sky_pass->set_sun_luminance(math::type_cast(blackbody.luminance)); + this->sky_pass->set_sun_illuminance(math::type_cast(blackbody_illuminance)); + this->sky_pass->set_sun_angular_radius(static_cast(blackbody_angular_radius)); } - const double3& blackbody_icrf_position = orbit.icrf_position; + const double3& blackbody_icrf_position = orbit.position; // Update diffuse reflectors this->registry.view().each( [&](entity::id entity_id, const auto& reflector_body, const auto& reflector_orbit, const auto& reflector, const auto& transform) { // Calculate distance to blackbody and direction of incoming light - double3 blackbody_light_direction_icrf = reflector_orbit.icrf_position - blackbody_icrf_position; + double3 blackbody_light_direction_icrf = reflector_orbit.position - blackbody_icrf_position; double blackbody_distance = math::length(blackbody_light_direction_icrf); blackbody_light_direction_icrf = blackbody_light_direction_icrf / blackbody_distance; @@ -288,7 +312,7 @@ void astronomy::update(double t, double dt) // Calculate blackbody illuminance - double3 view_direction_icrf = reflector_orbit.icrf_position - reference_orbit.icrf_position; + double3 view_direction_icrf = reflector_orbit.position - reference_orbit.position; const double reflector_distance = math::length(view_direction_icrf); view_direction_icrf = view_direction_icrf / reflector_distance; @@ -300,7 +324,7 @@ void astronomy::update(double t, double dt) const double3 planet_luminance = (sunlight_illuminance * reference_body.albedo) / math::pi; const double planet_angular_radius = astro::angular_radius(reference_body.radius, reflector_distance); const double planet_solid_angle = geom::solid_angle::cone(planet_angular_radius); - const double planet_phase_factor = math::dot(-view_direction_icrf, math::normalize(reference_orbit.icrf_position - blackbody_icrf_position)) * 0.5 + 0.5; + const double planet_phase_factor = math::dot(-view_direction_icrf, math::normalize(reference_orbit.position - blackbody_icrf_position)) * 0.5 + 0.5; const double3 planetlight_illuminance = planet_luminance * planet_solid_angle * planet_phase_factor; double3 planetlight_direction_eus = math::normalize(icrf_to_eus.r * view_direction_icrf); @@ -309,12 +333,16 @@ void astronomy::update(double t, double dt) const double3 reflected_planetlight_luminance = (planetlight_illuminance * reflector.albedo) / math::pi; const double3 reflected_planetlight_illuminance = reflected_planetlight_luminance * reflector_solid_angle; - std::cout << "sunlight illuminance: " << sunlight_illuminance << std::endl; + /* std::cout << "reflected sunlight illuminance: " << reflected_sunlight_illuminance << std::endl; std::cout << "planetlight illuminance: " << planetlight_illuminance << std::endl; + std::cout << "planet luminance: " << planet_luminance << std::endl; + std::cout << "reflected planetlight luminance: " << reflected_planetlight_luminance << std::endl; std::cout << "reflected planetlight illuminance: " << reflected_planetlight_illuminance << std::endl; std::cout << "reflector phase: " << reflector_phase_factor << std::endl; std::cout << "planet phase: " << planet_phase_factor << std::endl; + */ + if (this->sky_pass) { @@ -332,7 +360,7 @@ void astronomy::update(double t, double dt) float3 reflector_up_eus = math::type_cast(icrf_to_eus.r * double3{0, 0, 1}); double3 reflected_illuminance = reflected_sunlight_illuminance + reflected_planetlight_illuminance; - reflected_illuminance *= std::max(0.0, std::sin(transform.local.translation.y)); + //reflected_illuminance *= std::max(0.0, std::sin(transform.local.translation.y)); total_illuminance += (reflected_illuminance.x + reflected_illuminance.y + reflected_illuminance.z) / 3.0; @@ -397,14 +425,15 @@ void astronomy::update(double t, double dt) { const double calibration = 250.0; const double ev100 = std::log2((total_illuminance * 100.0) / calibration); - // std::cout << "EV100: " << ev100 << std::endl; - camera->set_exposure(static_cast(ev100)); + //std::cout << "LUX: " << total_illuminance << std::endl; + //std::cout << "EV100: " << ev100 << std::endl; + //camera->set_exposure(exposure_offset); } } -void astronomy::set_universal_time(double time) +void astronomy::set_time(double time) { - universal_time = time; + this->time = time; } void astronomy::set_time_scale(double scale) @@ -444,6 +473,11 @@ void astronomy::set_camera(scene::camera* camera) this->camera = camera; } +void astronomy::set_exposure_offset(float offset) +{ + exposure_offset = offset; +} + void astronomy::set_sky_pass(::render::sky_pass* pass) { this->sky_pass = pass; diff --git a/src/entity/systems/astronomy.hpp b/src/entity/systems/astronomy.hpp index 015c3a3..d8a0769 100644 --- a/src/entity/systems/astronomy.hpp +++ b/src/entity/systems/astronomy.hpp @@ -53,14 +53,14 @@ public: virtual void update(double t, double dt); /** - * Sets the current universal time. + * Sets the current time. * - * @param time Universal time, in days. + * @param time Time, in days. */ - void set_universal_time(double time); + void set_time(double time); /** - * Sets the factor by which the timestep `dt` will be scaled before being added to the current universal time. + * Sets the factor by which the timestep `dt` will be scaled before being added to the current time. * * @param scale Factor by which to scale the timestep. */ @@ -84,6 +84,9 @@ public: void set_sky_light(scene::ambient_light* light); void set_moon_light(scene::directional_light* light); void set_camera(scene::camera* camera); + void set_exposure_offset(float offset); + + inline float get_exposure_offset() const { return exposure_offset; }; void set_sky_pass(::render::sky_pass* pass); @@ -93,7 +96,7 @@ private: void update_bcbf_to_enu(); - double universal_time; + double time; double time_scale; entity::id reference_entity; @@ -111,6 +114,7 @@ private: scene::directional_light* moon_light; scene::camera* camera; ::render::sky_pass* sky_pass; + float exposure_offset; }; } // namespace system diff --git a/src/entity/systems/orbit.cpp b/src/entity/systems/orbit.cpp index 003733d..fb8a8ba 100644 --- a/src/entity/systems/orbit.cpp +++ b/src/entity/systems/orbit.cpp @@ -26,83 +26,53 @@ namespace system { orbit::orbit(entity::registry& registry): updatable(registry), - universal_time(0.0), - time_scale(1.0), - ke_iterations(10), - ke_tolerance(1e-6) -{ - registry.on_construct().connect<&orbit::on_orbit_construct>(this); - registry.on_replace().connect<&orbit::on_orbit_replace>(this); -} + ephemeris(nullptr), + time(0.0), + time_scale(1.0) +{} void orbit::update(double t, double dt) { // Add scaled timestep to current time - set_universal_time(universal_time + dt * time_scale); + set_time(time + dt * time_scale); - // Propagate orbits - registry.view().each( - [&](entity::id entity_id, auto& orbit) - { - // 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::anomaly::mean_to_eccentric(orbit.elements.ec, ma, ke_iterations, ke_tolerance); - - // Calculate Cartesian orbital position in the PQW frame - math::vector3 pqw_position = physics::orbit::frame::pqw::cartesian(orbit.elements.ec, orbit.elements.a, ea, orbit.semiminor_axis); - - // Transform orbital position from PQW frame to BCI frame - orbit.bci_position = orbit.pqw_to_bci.transform(pqw_position); - }); + if (!ephemeris) + return; - // Update orbital positions in the ICRF frame + // Calculate positions of ephemeris items, in meters + for (std::size_t i = 0; i < ephemeris->size(); ++i) + positions[i] = (*ephemeris)[i].position(time) * 1000.0; + + // Propagate orbits registry.view().each( - [&](entity::id entity_id, auto& orbit) + [&](entity::id entity_eid, auto& orbit) { - orbit.icrf_position = orbit.bci_position; + orbit.position = positions[orbit.ephemeris_index] * orbit.scale; - entity::id parent = orbit.parent; - while (parent != entt::null) + entity::id parent_id = orbit.parent; + while (parent_id != entt::null) { - const component::orbit& parent_orbit = registry.get(parent); - orbit.icrf_position += parent_orbit.bci_position; - parent = parent_orbit.parent; + const component::orbit& parent_orbit = registry.get(parent_id); + orbit.position += positions[parent_orbit.ephemeris_index] * parent_orbit.scale; + parent_id = parent_orbit.parent; } }); } -void orbit::set_universal_time(double time) +void orbit::set_ephemeris(const physics::orbit::ephemeris* ephemeris) { - universal_time = time; + this->ephemeris = ephemeris; + positions.resize((ephemeris) ? ephemeris->size() : 0); } -void orbit::set_time_scale(double scale) +void orbit::set_time(double time) { - time_scale = scale; + this->time = time; } -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) +void orbit::set_time_scale(double scale) { - 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 - ); + time_scale = scale; } } // namespace system diff --git a/src/entity/systems/orbit.hpp b/src/entity/systems/orbit.hpp index 7b793b1..a6069a6 100644 --- a/src/entity/systems/orbit.hpp +++ b/src/entity/systems/orbit.hpp @@ -24,6 +24,7 @@ #include "utility/fundamental-types.hpp" #include "entity/id.hpp" #include "entity/components/orbit.hpp" +#include "physics/orbit/ephemeris.hpp" namespace entity { namespace system { @@ -46,27 +47,31 @@ public: virtual void update(double t, double dt); /** - * Sets the current universal time. + * Sets the current time. * - * @param time Universal time, in days. + * @param time Time, in days. */ - void set_universal_time(double time); + void set_time(double time); /** - * Sets the factor by which the timestep `dt` will be scaled before being added to the current universal time. + * Sets the factor by which the timestep `dt` will be scaled before being added to the current time. * * @param scale Factor by which to scale the timestep. */ 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); + /** + * Sets the ephemeris used to calculate orbital positions. + * + * @param ephemeris Ephemeris. + */ + void set_ephemeris(const physics::orbit::ephemeris* ephemeris); - double universal_time; +private: + const physics::orbit::ephemeris* ephemeris; + double time; double time_scale; - std::size_t ke_iterations; - double ke_tolerance; + std::vector positions; }; } // namespace system diff --git a/src/entity/systems/terrain.cpp b/src/entity/systems/terrain.cpp index 27220c2..0eb9ee7 100644 --- a/src/entity/systems/terrain.cpp +++ b/src/entity/systems/terrain.cpp @@ -48,7 +48,7 @@ terrain::terrain(entity::registry& registry): max_error(0.0) { // Build set of quaternions to rotate quadtree cube coordinates into BCBF space according to face index - face_rotations[0] = math::quaternion::identity(); // +x + face_rotations[0] = math::quaternion::identity; // +x face_rotations[1] = math::quaternion::rotate_z(math::pi); // -x face_rotations[2] = math::quaternion::rotate_z( math::half_pi); // +y face_rotations[3] = math::quaternion::rotate_z(-math::half_pi); // -y diff --git a/src/game/ant/morphogenesis.cpp b/src/game/ant/morphogenesis.cpp index 576d3d7..25d9a18 100644 --- a/src/game/ant/morphogenesis.cpp +++ b/src/game/ant/morphogenesis.cpp @@ -563,7 +563,7 @@ render::model* build_model + sting_index_count; // Calculate transform from legs space to body space - const math::transform& legs_to_body = math::identity_transform; + const math::transform& legs_to_body = math::transform::identity; // Reskin leg bones std::unordered_set old_foreleg_coxa_l_indices; diff --git a/src/game/state/boot.cpp b/src/game/state/boot.cpp index cc61f4e..fdf47db 100644 --- a/src/game/state/boot.cpp +++ b/src/game/state/boot.cpp @@ -492,7 +492,7 @@ void boot::setup_rendering() ctx.sky_pass = new render::sky_pass(ctx.rasterizer, ctx.hdr_framebuffer, ctx.resource_manager); ctx.sky_pass->set_enabled(false); - ctx.sky_pass->set_magnification(4.0f); + ctx.sky_pass->set_magnification(10.0f); ctx.app->get_event_dispatcher()->subscribe(ctx.sky_pass); ctx.ground_pass = new render::ground_pass(ctx.rasterizer, ctx.hdr_framebuffer, ctx.resource_manager); diff --git a/src/game/state/brood.cpp b/src/game/state/brood.cpp index d5a8d3e..d774a1f 100644 --- a/src/game/state/brood.cpp +++ b/src/game/state/brood.cpp @@ -88,7 +88,7 @@ void setup_nest(game::context* ctx) ctx->entities["shaft"] = shaft_eid; entity::component::transform transform; - transform.local = math::identity_transform; + transform.local = math::transform::identity; transform.world = transform.local; transform.warp = true; @@ -134,7 +134,7 @@ void setup_camera(game::context* ctx) { // Transform entity::component::transform target_transform; - target_transform.local = math::identity_transform; + target_transform.local = math::transform::identity; target_transform.world = target_transform.local; target_transform.warp = true; ctx->entity_registry->assign(target_eid, target_transform); @@ -146,7 +146,7 @@ void setup_camera(game::context* ctx) // Create camera transform component entity::component::transform transform; - transform.local = math::identity_transform; + transform.local = math::transform::identity; transform.world = transform.local; transform.warp = true; ctx->entity_registry->assign(camera_eid, transform); diff --git a/src/game/state/forage.cpp b/src/game/state/forage.cpp index 402599d..e95510b 100644 --- a/src/game/state/forage.cpp +++ b/src/game/state/forage.cpp @@ -130,7 +130,7 @@ void setup_camera(game::context* ctx) { // Transform entity::component::transform target_transform; - target_transform.local = math::identity_transform; + target_transform.local = math::transform::identity; target_transform.world = target_transform.local; target_transform.warp = true; ctx->entity_registry->assign(target_eid, target_transform); @@ -142,7 +142,7 @@ void setup_camera(game::context* ctx) // Create camera transform component entity::component::transform transform; - transform.local = math::identity_transform; + transform.local = math::transform::identity; transform.world = transform.local; transform.warp = true; ctx->entity_registry->assign(camera_eid, transform); diff --git a/src/game/state/nuptial-flight.cpp b/src/game/state/nuptial-flight.cpp index 51c8a8d..9cf6ce9 100644 --- a/src/game/state/nuptial-flight.cpp +++ b/src/game/state/nuptial-flight.cpp @@ -44,6 +44,7 @@ #include "game/ant/breed.hpp" #include "game/ant/morphogenesis.hpp" #include "physics/time/time.hpp" +#include "math/interpolation.hpp" #include using namespace game::ant; @@ -92,7 +93,9 @@ nuptial_flight::nuptial_flight(game::context& ctx): // Create world if not yet created if (ctx.entities.find("planet") == ctx.entities.end()) - { + { + // Create cosmos + game::world::load_ephemeris(ctx); game::world::create_stars(ctx); game::world::create_sun(ctx); game::world::create_em_bary(ctx); @@ -103,10 +106,10 @@ nuptial_flight::nuptial_flight(game::context& ctx): const double utc = 0.0; const double equinox_2000 = physics::time::gregorian::to_ut1(2000, 1, 1, 12, 0, 0.0, utc); const double spring_2000 = physics::time::gregorian::to_ut1(2000, 6, 19, 12, 0, 0.0, utc); - game::world::set_time(ctx, 55.0); + game::world::set_time(ctx, 14622.5); // Freeze time - game::world::set_time_scale(ctx, 0.0); + game::world::set_time_scale(ctx, 60.0); } // Load biome @@ -163,8 +166,8 @@ nuptial_flight::nuptial_flight(game::context& ctx): ctx.entity_registry->assign(boid_eid, model); entity::component::transform transform; - transform.local = math::identity_transform; - transform.world = math::identity_transform; + transform.local = math::transform::identity; + transform.world = math::transform::identity; transform.warp = true; ctx.entity_registry->assign(boid_eid, transform); @@ -216,7 +219,7 @@ void nuptial_flight::setup_camera() { // Transform entity::component::transform target_transform; - target_transform.local = math::identity_transform; + target_transform.local = math::transform::identity; target_transform.world = target_transform.local; target_transform.warp = true; ctx.entity_registry->assign(target_eid, target_transform); @@ -228,7 +231,7 @@ void nuptial_flight::setup_camera() // Create camera transform component entity::component::transform transform; - transform.local = math::identity_transform; + transform.local = math::transform::identity; transform.world = transform.local; transform.warp = true; ctx.entity_registry->assign(camera_eid, transform); @@ -313,27 +316,28 @@ void nuptial_flight::enable_keeper_controls() bool gamepad_invert_pan = false; bool mouse_look_toggle = false; ctx.mouse_look = false; - const double time_scale = 10000.0; + const double time_scale = 60.0; + const double ff_time_scale = 50000.0; if (ctx.config->contains("mouse_tilt_sensitivity")) mouse_tilt_sensitivity = math::radians((*ctx.config)["mouse_tilt_sensitivity"].get()); if (ctx.config->contains("mouse_pan_sensitivity")) mouse_pan_sensitivity = math::radians((*ctx.config)["mouse_pan_sensitivity"].get()); if (ctx.config->contains("mouse_invert_tilt")) - mouse_invert_tilt = math::radians((*ctx.config)["mouse_invert_tilt"].get()); + mouse_invert_tilt = (*ctx.config)["mouse_invert_tilt"].get(); if (ctx.config->contains("mouse_invert_pan")) - mouse_invert_pan = math::radians((*ctx.config)["mouse_invert_pan"].get()); + mouse_invert_pan = (*ctx.config)["mouse_invert_pan"].get(); if (ctx.config->contains("mouse_look_toggle")) - mouse_look_toggle = math::radians((*ctx.config)["mouse_look_toggle"].get()); + mouse_look_toggle = (*ctx.config)["mouse_look_toggle"].get(); if (ctx.config->contains("gamepad_tilt_sensitivity")) gamepad_tilt_sensitivity = math::radians((*ctx.config)["gamepad_tilt_sensitivity"].get()); if (ctx.config->contains("gamepad_pan_sensitivity")) gamepad_pan_sensitivity = math::radians((*ctx.config)["gamepad_pan_sensitivity"].get()); if (ctx.config->contains("gamepad_invert_tilt")) - gamepad_invert_tilt = math::radians((*ctx.config)["gamepad_invert_tilt"].get()); + gamepad_invert_tilt = (*ctx.config)["gamepad_invert_tilt"].get(); if (ctx.config->contains("gamepad_invert_pan")) - gamepad_invert_pan = math::radians((*ctx.config)["gamepad_invert_pan"].get()); + gamepad_invert_pan = (*ctx.config)["gamepad_invert_pan"].get(); const input::control* move_slow = ctx.controls["move_slow"]; const input::control* move_fast = ctx.controls["move_fast"]; @@ -581,30 +585,30 @@ void nuptial_flight::enable_keeper_controls() // Fast-forward ctx.controls["fast_forward"]->set_activated_callback ( - [&ctx = this->ctx, time_scale]() + [&ctx = this->ctx, ff_time_scale]() { - game::world::set_time_scale(ctx, time_scale); + game::world::set_time_scale(ctx, ff_time_scale); } ); ctx.controls["fast_forward"]->set_deactivated_callback ( [&ctx = this->ctx, time_scale]() { - game::world::set_time_scale(ctx, 0.0); + game::world::set_time_scale(ctx, time_scale); } ); ctx.controls["rewind"]->set_activated_callback ( - [&ctx = this->ctx, time_scale]() + [&ctx = this->ctx, ff_time_scale]() { - game::world::set_time_scale(ctx, -time_scale); + game::world::set_time_scale(ctx, -ff_time_scale); } ); ctx.controls["rewind"]->set_deactivated_callback ( [&ctx = this->ctx, time_scale]() { - game::world::set_time_scale(ctx, 0.0); + game::world::set_time_scale(ctx, time_scale); } ); @@ -627,6 +631,26 @@ void nuptial_flight::enable_keeper_controls() ctx.state_machine.emplace(new game::state::pause_menu(ctx)); } ); + + ctx.controls["increase_exposure"]->set_activated_callback + ( + [&ctx = this->ctx]() + { + //ctx.astronomy_system->set_exposure_offset(ctx.astronomy_system->get_exposure_offset() - 1.0f); + ctx.surface_camera->set_exposure(ctx.surface_camera->get_exposure() + 1.0f); + ctx.logger->log("EV100: " + std::to_string(ctx.surface_camera->get_exposure())); + } + ); + + ctx.controls["decrease_exposure"]->set_activated_callback + ( + [&ctx = this->ctx]() + { + //ctx.astronomy_system->set_exposure_offset(ctx.astronomy_system->get_exposure_offset() + 1.0f); + ctx.surface_camera->set_exposure(ctx.surface_camera->get_exposure() - 1.0f); + ctx.logger->log("EV100: " + std::to_string(ctx.surface_camera->get_exposure())); + } + ); } void nuptial_flight::disable_keeper_controls() @@ -651,6 +675,8 @@ void nuptial_flight::disable_keeper_controls() ctx.controls["fast_forward"]->set_activated_callback(nullptr); ctx.controls["rewind"]->set_activated_callback(nullptr); ctx.controls["pause"]->set_activated_callback(nullptr); + ctx.controls["increase_exposure"]->set_activated_callback(nullptr); + ctx.controls["decrease_exposure"]->set_activated_callback(nullptr); } void nuptial_flight::enable_ant_controls() @@ -672,25 +698,26 @@ void nuptial_flight::enable_ant_controls() float gamepad_pan_sensitivity = 1.0f; bool gamepad_invert_tilt = false; bool gamepad_invert_pan = false; - const double time_scale = 10000.0; + const double time_scale = 60.0; + const double ff_time_scale = 50000.0; if (ctx.config->contains("mouse_tilt_sensitivity")) mouse_tilt_sensitivity = math::radians((*ctx.config)["mouse_tilt_sensitivity"].get()); if (ctx.config->contains("mouse_pan_sensitivity")) mouse_pan_sensitivity = math::radians((*ctx.config)["mouse_pan_sensitivity"].get()); if (ctx.config->contains("mouse_invert_tilt")) - mouse_invert_tilt = math::radians((*ctx.config)["mouse_invert_tilt"].get()); + mouse_invert_tilt = (*ctx.config)["mouse_invert_tilt"].get(); if (ctx.config->contains("mouse_invert_pan")) - mouse_invert_pan = math::radians((*ctx.config)["mouse_invert_pan"].get()); + mouse_invert_pan = (*ctx.config)["mouse_invert_pan"].get(); if (ctx.config->contains("gamepad_tilt_sensitivity")) gamepad_tilt_sensitivity = math::radians((*ctx.config)["gamepad_tilt_sensitivity"].get()); if (ctx.config->contains("gamepad_pan_sensitivity")) gamepad_pan_sensitivity = math::radians((*ctx.config)["gamepad_pan_sensitivity"].get()); if (ctx.config->contains("gamepad_invert_tilt")) - gamepad_invert_tilt = math::radians((*ctx.config)["gamepad_invert_tilt"].get()); + gamepad_invert_tilt = (*ctx.config)["gamepad_invert_tilt"].get(); if (ctx.config->contains("gamepad_invert_pan")) - gamepad_invert_pan = math::radians((*ctx.config)["gamepad_invert_pan"].get()); + gamepad_invert_pan = (*ctx.config)["gamepad_invert_pan"].get(); const input::control* move_slow = ctx.controls["move_slow"]; const input::control* move_fast = ctx.controls["move_fast"]; @@ -880,30 +907,30 @@ void nuptial_flight::enable_ant_controls() // Fast-forward ctx.controls["fast_forward"]->set_activated_callback ( - [&ctx = this->ctx, time_scale]() + [&ctx = this->ctx, ff_time_scale]() { - game::world::set_time_scale(ctx, time_scale); + game::world::set_time_scale(ctx, ff_time_scale); } ); ctx.controls["fast_forward"]->set_deactivated_callback ( [&ctx = this->ctx, time_scale]() { - game::world::set_time_scale(ctx, 0.0); + game::world::set_time_scale(ctx, time_scale); } ); ctx.controls["rewind"]->set_activated_callback ( - [&ctx = this->ctx, time_scale]() + [&ctx = this->ctx, ff_time_scale]() { - game::world::set_time_scale(ctx, -time_scale); + game::world::set_time_scale(ctx, -ff_time_scale); } ); ctx.controls["rewind"]->set_deactivated_callback ( [&ctx = this->ctx, time_scale]() { - game::world::set_time_scale(ctx, 0.0); + game::world::set_time_scale(ctx, time_scale); } ); diff --git a/src/game/world.cpp b/src/game/world.cpp index 53299b6..3fdccde 100644 --- a/src/game/world.cpp +++ b/src/game/world.cpp @@ -19,7 +19,7 @@ #include "game/world.hpp" #include "scene/text.hpp" -#include "astro/illuminance.hpp" +#include "physics/light/vmag.hpp" #include "color/color.hpp" #include "entity/components/atmosphere.hpp" #include "entity/components/blackbody.hpp" @@ -38,7 +38,7 @@ #include "gl/vertex-buffer.hpp" #include "physics/light/photometry.hpp" #include "physics/orbit/orbit.hpp" -#include "astro/illuminance.hpp" +#include "physics/orbit/ephemeris.hpp" #include "render/material.hpp" #include "render/model.hpp" #include "render/passes/shadow-map-pass.hpp" @@ -55,6 +55,12 @@ namespace game { namespace world { +void load_ephemeris(game::context& ctx) +{ + // Load ephemeris + ctx.orbit_system->set_ephemeris(ctx.resource_manager->load>("de421.eph")); +} + void create_stars(game::context& ctx) { // Load star catalog @@ -89,12 +95,12 @@ void create_stars(game::context& ctx) vmag = std::stod(catalog_row[3]); bv_color = std::stod(catalog_row[4]); } - catch (const std::exception& e) + catch (const std::exception&) { continue; } - starlight_illuminance += astro::vmag_to_lux(vmag); + starlight_illuminance += physics::light::vmag::to_illuminance(vmag); // Convert right ascension and declination from degrees to radians ra = math::wrap_radians(math::radians(ra)); @@ -113,7 +119,7 @@ void create_stars(game::context& ctx) double3 color_acescg = color::xyz::to_acescg(color_xyz); // Convert apparent magnitude to brightness factor relative to a 0th magnitude star - double brightness = astro::vmag_to_brightness(vmag); + double brightness = physics::light::vmag::to_brightness(vmag); // Scale color by relative brightness double3 scaled_color = color_acescg * brightness; @@ -193,11 +199,12 @@ void create_sun(game::context& ctx) // Create sun directional light scene object scene::directional_light* sun_light = new scene::directional_light(); + sun_light->set_color({0, 0, 0}); + sun_light->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->set_color({0, 0, 0}); sky_light->update_tweens(); // Add sun light scene objects to surface scene @@ -216,9 +223,6 @@ void create_em_bary(game::context& ctx) entity::archetype* em_bary_archetype = ctx.resource_manager->load("em-bary.ent"); entity::id em_bary_eid = em_bary_archetype->create(*ctx.entity_registry); ctx.entities["em_bary"] = em_bary_eid; - - // Assign orbital parent - ctx.entity_registry->get(em_bary_eid).parent = ctx.entities["sun"]; } void create_earth(game::context& ctx) @@ -261,8 +265,7 @@ void create_moon(game::context& ctx) // Create moon directional light scene object scene::directional_light* moon_light = new scene::directional_light(); - moon_light->set_color({1, 1, 1}); - moon_light->set_intensity(1.0f); + moon_light->set_color({0, 0, 0}); moon_light->update_tweens(); // Add moon light scene objects to surface scene @@ -274,8 +277,8 @@ void create_moon(game::context& ctx) void set_time(game::context& ctx, double t) { - ctx.astronomy_system->set_universal_time(t); - ctx.orbit_system->set_universal_time(t); + ctx.astronomy_system->set_time(t); + ctx.orbit_system->set_time(t); } void set_time_scale(game::context& ctx, double scale) diff --git a/src/game/world.hpp b/src/game/world.hpp index 68d4c4f..8e3ccff 100644 --- a/src/game/world.hpp +++ b/src/game/world.hpp @@ -27,6 +27,9 @@ namespace game { /// World creation and manipulation functions. namespace world { +/// Creates an ephemeris. +void load_ephemeris(game::context& ctx); + /// Creates the fixed stars. void create_stars(game::context& ctx); diff --git a/src/geom/view-frustum.hpp b/src/geom/view-frustum.hpp index a6704e8..15fb5db 100644 --- a/src/geom/view-frustum.hpp +++ b/src/geom/view-frustum.hpp @@ -99,7 +99,7 @@ view_frustum::view_frustum(const matrix_type& view_projection): template view_frustum::view_frustum(): - view_frustum(math::identity4x4) + view_frustum(math::matrix4::identity) {} template diff --git a/src/math/constants.hpp b/src/math/constants.hpp index 06dd6af..21899f0 100644 --- a/src/math/constants.hpp +++ b/src/math/constants.hpp @@ -26,73 +26,25 @@ namespace math { -/** - * Pi. - */ -template -constexpr T pi = T(3.14159265358979323846); - -/** - * Pi / 2. - */ +/// Pi. template -constexpr T half_pi = pi * T(0.5); +constexpr T pi = T(3.1415926535897932384626433832795); -/** - * Pi * 2. - */ +/// Pi * 2. template constexpr T two_pi = pi * T(2); -/** - * 2x2 identity matrix. - */ -template -constexpr matrix identity2x2 = -{{ - {1, 0}, - {0, 1} -}}; - -/** - * 3x3 identity matrix. - */ -template -constexpr matrix identity3x3 = -{{ - {1, 0, 0}, - {0, 1, 0}, - {0, 0, 1} -}}; - -/** - * 4x4 identity matrix. - */ +/// Pi * 4. template -constexpr matrix identity4x4 = -{{ - {1, 0, 0, 0}, - {0, 1, 0, 0}, - {0, 0, 1, 0}, - {0, 0, 0, 1} -}}; +constexpr T four_pi = pi * T(4); -/** - * Unit quaternion. - */ +/// Pi / 2. template -constexpr quaternion identity_quaternion = {T(1), T(0), T(0), T(0)}; +constexpr T half_pi = pi * T(0.5); -/** - * Identity transform. - */ +/// 1 / Pi. template -constexpr transform identity_transform = -{ - {T(0), T(0), T(0)}, - {T(1), T(0), T(0), T(0)}, - {T(1), T(1), T(1)} -}; +constexpr T inverse_pi = T(1) / pi; } // namespace math diff --git a/src/math/interpolation.hpp b/src/math/interpolation.hpp index ce70720..5d29f82 100644 --- a/src/math/interpolation.hpp +++ b/src/math/interpolation.hpp @@ -20,12 +20,8 @@ #ifndef ANTKEEPER_MATH_INTERPOLATION_HPP #define ANTKEEPER_MATH_INTERPOLATION_HPP -#include "math/constants.hpp" -#include "math/matrix-type.hpp" -#include "math/quaternion-type.hpp" -#include "math/transform-type.hpp" +#include "math/angles.hpp" #include -#include namespace math { @@ -76,8 +72,7 @@ inline T lerp(const T& x, const T& y, S a) template T lerp_angle(T x, T y, T a) { - T shortest_angle = std::remainder((y - x), two_pi); - return std::remainder(x + shortest_angle * a, two_pi); + return wrap_radians(x + wrap_radians(y - x) * a); } template diff --git a/src/math/matrix-type.hpp b/src/math/matrix-type.hpp index d60d1f1..e00939d 100644 --- a/src/math/matrix-type.hpp +++ b/src/math/matrix-type.hpp @@ -22,6 +22,7 @@ #include "math/vector-type.hpp" #include +#include namespace math { @@ -38,11 +39,35 @@ struct matrix typedef T element_type; typedef vector row_type; row_type columns[N]; + + /// Identity matrix. + static const matrix identity; inline constexpr row_type& operator[](std::size_t i) noexcept { return columns[i]; } inline constexpr const row_type& operator[](std::size_t i) const noexcept { return columns[i]; } }; +template +constexpr vector identity_matrix_column(const std::index_sequence&) +{ + return {(Is == I ? T{1} : T{0})...}; +} + +template +constexpr matrix identity_matrix_rows(const std::index_sequence& is) +{ + return {{identity_matrix_column(is)...}}; +} + +template +constexpr matrix identity_matrix() +{ + return identity_matrix_rows(std::make_index_sequence{}); +} + +template +constexpr matrix matrix::identity = identity_matrix(); + /// 2x2 matrix. template using matrix2 = matrix; @@ -58,4 +83,3 @@ using matrix4 = matrix; } // namespace math #endif // ANTKEEPER_MATH_MATRIX_TYPE_HPP - diff --git a/src/math/polynomial.hpp b/src/math/polynomial.hpp index 50a3753..a489807 100644 --- a/src/math/polynomial.hpp +++ b/src/math/polynomial.hpp @@ -96,12 +96,15 @@ namespace chebyshev { template T evaluate(InputIt first, InputIt last, T x) { - T y = *(first++) * T(0.5) + *(first++) * x; + T y = *(first++); + y += *(first++) * x; - const T x2 = x * T(2); - for (T n2 = T(1), n1 = x, n0; first != last; n2 = n1, n1 = n0) + T n2 = T(1), n1 = x, n0; + x *= T(2); + + for (; first != last; n2 = n1, n1 = n0) { - n0 = x2 * n1 - n2; + n0 = x * n1 - n2; y += *(first++) * n0; } diff --git a/src/math/quaternion-type.hpp b/src/math/quaternion-type.hpp index f895155..14eb95f 100644 --- a/src/math/quaternion-type.hpp +++ b/src/math/quaternion-type.hpp @@ -38,8 +38,8 @@ struct quaternion scalar_type y; scalar_type z; - /// Returns the rotation identity quaternion. - static constexpr quaternion identity(); + /// Rotation identity quaternion. + static const quaternion identity; /** * Returns a quaternion representing a rotation about the x-axis. @@ -67,10 +67,7 @@ struct quaternion }; template -constexpr quaternion quaternion::identity() -{ - return {T(1), T(0), T(0), T(0)}; -}; +constexpr quaternion quaternion::identity = {T(1), T(0), T(0), T(0)}; template quaternion quaternion::rotate_x(scalar_type angle) diff --git a/src/math/transform-type.hpp b/src/math/transform-type.hpp index 249ee6c..dd213dc 100644 --- a/src/math/transform-type.hpp +++ b/src/math/transform-type.hpp @@ -41,9 +41,19 @@ struct transform /// Scale vector vector scale; + + /// Identity transform. + static const transform identity; +}; + +template +constexpr transform transform::identity = +{ + {T(0), T(0), T(0)}, + {T(1), T(0), T(0), T(0)}, + {T(1), T(1), T(1)} }; } // namespace math #endif // ANTKEEPER_MATH_TRANSFORM_TYPE_HPP - diff --git a/src/physics/light/exposure.hpp b/src/physics/light/exposure.hpp new file mode 100644 index 0000000..f74db89 --- /dev/null +++ b/src/physics/light/exposure.hpp @@ -0,0 +1,114 @@ +/* + * 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 . + */ + +#ifndef ANTKEEPER_PHYSICS_LIGHT_EV_HPP +#define ANTKEEPER_PHYSICS_LIGHT_EV_HPP + +#include + +namespace physics { +namespace light { + +/** + * Exposure value. + * + * @see https://en.wikipedia.org/wiki/Exposure_value + */ +namespace ev { + +/** + * Exposure value from luminance. + * + * @param l Luminance, in cd/m^2. + * @param s ISO arithmetic speed. A value of `100` corresponds to ISO 100. + * @param k Reflected-light meter calibration constant. A common value is `12.5`. + * + * @return Exposure value. + */ +template +T from_luminance(T l, T s, T k) +{ + return std::log2((l * s) / k); +} + +/** + * Exposure value from illuminance. + * + * @param e Illuminance, in lux. + * @param s ISO arithmetic speed. A value of `100` corresponds to ISO 100. + * @param c Incident-light meter calibration constant. A common value is `250`. + * + * @return Exposure value. + */ +template +T from_illuminance(T e, T s, T c) +{ + return std::log2((e * s) / c); +} + +/** + * Exposure value from exposure settings. + * + * @param n Relative aperture (f-number). + * @param t Exposure time (shutter speed), in seconds. + * @param s ISO arithmetic speed. A value of `100` corresponds to ISO 100. + * + * @return Exposure value. + */ +template +T from_settings(T n, T t, T s) +{ + return std::log2((n * n) / t * T(100) / s); +} + +/** + * Exposure value to luminance. + * + * @param ev Exposure value. + * @param s ISO arithmetic speed. A value of `100` corresponds to ISO 100. + * @param k Reflected-light meter calibration constant. A common value is `12.5`. + * + * @return Luminance, in cd/m^2. + */ +template +T to_luminance(T ev, T s, T k) +{ + return (k * std::exp2(ev)) / s; +} + +/** + * Exposure value to luminance. + * + * @param ev Exposure value. + * @param s ISO arithmetic speed. A value of `100` corresponds to ISO 100. + * @param k Incident-light meter calibration constant. A common value is `250`. + * + * @return Illuminance, in lux. + */ +template +T to_luminance(T ev, T s, T c) +{ + return (c * std::exp2(ev)) / s; +} + +} // namespace ev +} // namespace light +} // namespace physics + +#endif // ANTKEEPER_PHYSICS_LIGHT_EV_HPP diff --git a/src/physics/light/light.hpp b/src/physics/light/light.hpp index 6e2b044..a9a7e35 100644 --- a/src/physics/light/light.hpp +++ b/src/physics/light/light.hpp @@ -32,5 +32,6 @@ namespace light {} #include "phase.hpp" #include "photometry.hpp" #include "refraction.hpp" +#include "vmag.hpp" #endif // ANTKEEPER_PHYSICS_LIGHT_HPP diff --git a/src/physics/light/phase.hpp b/src/physics/light/phase.hpp index f3a3223..42e28de 100644 --- a/src/physics/light/phase.hpp +++ b/src/physics/light/phase.hpp @@ -50,7 +50,7 @@ template T henyey_greenstein(T mu, T g); /** - * Isotropic phase function. Causes light to be scattered uniformly in all directions. + * Isotropic phase function. */ template constexpr T isotropic() diff --git a/src/astro/illuminance.cpp b/src/physics/light/vmag.hpp similarity index 50% rename from src/astro/illuminance.cpp rename to src/physics/light/vmag.hpp index 3d88192..4558a7c 100644 --- a/src/astro/illuminance.cpp +++ b/src/physics/light/vmag.hpp @@ -17,27 +17,63 @@ * along with Antkeeper source code. If not, see . */ -#include "illuminance.hpp" +#ifndef ANTKEEPER_PHYSICS_LIGHT_VMAG_HPP +#define ANTKEEPER_PHYSICS_LIGHT_VMAG_HPP + #include -namespace astro -{ +namespace physics { +namespace light { -double vmag_to_brightness(double mv) +/// Apparent (visual) magnitude functions. +namespace vmag { + +/** + * Converts apparent magnitude to a brightness factor relative to a 0th magnitude star. + * + * @param mv Apparent magnitude. + * @return Brightness factor relative to a 0th magnitude star. + * + * @see https://en.wikipedia.org/wiki/Illuminance + */ +template +T to_brightness(T mv) { // 100^(1/5) static constexpr double fifth_root_100 = 2.5118864315095801110850320677993; return std::pow(fifth_root_100, -mv); } -double vmag_to_lux(double mv) +/** + * Converts apparent magnitude to illuminance. + * + * @param mv Apparent magnitude. + * @return Illuminance, in lux. + * + * @see https://en.wikipedia.org/wiki/Illuminance + */ +template +T to_illuminance(T mv) { return std::pow(10.0, (-14.18 - mv) * 0.4); } -double lux_to_vmag(double ev) +/** + * Converts illuminance to apparent magnitude. + * + * @param ev Illuminance, in lux. + * @return Apparent magnitude. + * + * @see https://en.wikipedia.org/wiki/Illuminance + */ +template +T from_illuminance(T ev) { return -14.18 - 2.5 * std::log10(ev); } -} // namespace astro +} // namespace vmag +} // namespace light +} // namespace physics + +#endif // ANTKEEPER_PHYSICS_LIGHT_VMAG_HPP diff --git a/src/astro/constants.hpp b/src/physics/orbit/ephemeris.hpp similarity index 66% rename from src/astro/constants.hpp rename to src/physics/orbit/ephemeris.hpp index 7aceb8f..cf2860a 100644 --- a/src/astro/constants.hpp +++ b/src/physics/orbit/ephemeris.hpp @@ -17,16 +17,23 @@ * along with Antkeeper source code. If not, see . */ -#ifndef ANTKEEPER_ASTRO_CONSTANTS_HPP -#define ANTKEEPER_ASTRO_CONSTANTS_HPP +#ifndef ANTKEEPER_PHYSICS_ORBIT_EPHEMERIS_HPP +#define ANTKEEPER_PHYSICS_ORBIT_EPHEMERIS_HPP -namespace astro -{ +#include "physics/orbit/trajectory.hpp" -/// Kilometers per astronomical unit -template -constexpr T km_per_au = 6.68459e-9; +namespace physics { +namespace orbit { -} // namespace astro +/** + * Table of orbital trajectories. + * + * @tparam t Real type. + */ +template +using ephemeris = std::vector>; + +} // namespace orbit +} // namespace physics -#endif // ANTKEEPER_ASTRO_CONSTANTS_HPP +#endif // ANTKEEPER_PHYSICS_ORBIT_EPHEMERIS_HPP diff --git a/src/physics/orbit/frame.hpp b/src/physics/orbit/frame.hpp index 82adad4..c6c0210 100644 --- a/src/physics/orbit/frame.hpp +++ b/src/physics/orbit/frame.hpp @@ -204,7 +204,10 @@ namespace bci { * @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. + * + * @see Archinal, B.A., A’Hearn, M.F., Bowell, E. et al. Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2009. Celest Mech Dyn Astr 109, 101–135 (2011). https://doi.org/10.1007/s10569-010-9320-4 */ template math::transformation::se3 to_bcbf(T ra, T dec, T w) @@ -212,7 +215,7 @@ namespace bci { const math::quaternion r = math::normalize ( math::quaternion::rotate_z(-math::half_pi - ra) * - math::quaternion::rotate_x(-math::half_pi + dec) * + math::quaternion::rotate_x(dec - math::half_pi) * math::quaternion::rotate_z(-w) ); @@ -289,7 +292,10 @@ namespace bcbf { * @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. + * + * @see Archinal, B.A., A’Hearn, M.F., Bowell, E. et al. Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2009. Celest Mech Dyn Astr 109, 101–135 (2011). https://doi.org/10.1007/s10569-010-9320-4 */ template math::transformation::se3 to_bci(T ra, T dec, T w) @@ -298,7 +304,7 @@ namespace bcbf { ( math::quaternion::rotate_z(w) * math::quaternion::rotate_x(math::half_pi - dec) * - math::quaternion::rotate_z(math::half_pi + ra) + math::quaternion::rotate_z(ra + math::half_pi) ); diff --git a/src/physics/orbit/trajectory.hpp b/src/physics/orbit/trajectory.hpp new file mode 100644 index 0000000..d4fcb6a --- /dev/null +++ b/src/physics/orbit/trajectory.hpp @@ -0,0 +1,84 @@ +/* + * 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 . + */ + +#ifndef ANTKEEPER_PHYSICS_ORBIT_TRAJECTORY_HPP +#define ANTKEEPER_PHYSICS_ORBIT_TRAJECTORY_HPP + +#include "math/polynomial.hpp" +#include + +namespace physics { +namespace orbit { + +/** + * Describes the trajectory of an orbit with Chebyshev polynomials. + * + * @tparam t Real type. + */ +template +struct trajectory +{ + /// Start time of the trajectory. + T t0; + + /// End time of the trajectory. + T t1; + + /// Time step duration. + T dt; + + /// Chebyshev polynomial degree. + std::size_t n; + + /// Chebyshev polynomial coefficients. + std::vector a; + + /** + * Calculates the Cartesian position of a trajectory at a given time. + * + * @param t Time, on `[t0, t1)`. + * @return Trajectory position at time @p t. + */ + math::vector position(T t) const; +}; + +template +math::vector trajectory::position(T t) const +{ + t -= t0; + std::size_t i = static_cast(t / dt); + + const T* ax = &a[i * n * 3]; + const T* ay = ax + n; + const T* az = ay + n; + + t = (t / dt - i) * T(2) - T(1); + + math::vector3 r; + r.x = math::polynomial::chebyshev::evaluate(ax, ay, t); + r.y = math::polynomial::chebyshev::evaluate(ay, az, t); + r.z = math::polynomial::chebyshev::evaluate(az, az + n, t); + + return r; +} + +} // namespace orbit +} // namespace physics + +#endif // ANTKEEPER_PHYSICS_ORBIT_TRAJECTORY_HPP diff --git a/src/render/passes/ground-pass.cpp b/src/render/passes/ground-pass.cpp index f49b706..1a8b5c8 100644 --- a/src/render/passes/ground-pass.cpp +++ b/src/render/passes/ground-pass.cpp @@ -85,7 +85,7 @@ void ground_pass::render(const render::context& ctx, render::queue& queue) const float clip_near = camera.get_clip_near_tween().interpolate(ctx.alpha); float clip_far = camera.get_clip_far_tween().interpolate(ctx.alpha); float3 model_scale = float3{1.0f, 1.0f, 1.0f} * (clip_near + clip_far) * 0.5f; - float4x4 model = math::scale(math::identity4x4, model_scale); + float4x4 model = math::scale(math::matrix4::identity, model_scale); float4x4 view = math::resize<4, 4>(math::resize<3, 3>(camera.get_view_tween().interpolate(ctx.alpha))); float4x4 model_view = view * model; float4x4 projection = camera.get_projection_tween().interpolate(ctx.alpha); diff --git a/src/render/passes/shadow-map-pass.cpp b/src/render/passes/shadow-map-pass.cpp index 4a9ed22..81c6ac0 100644 --- a/src/render/passes/shadow-map-pass.cpp +++ b/src/render/passes/shadow-map-pass.cpp @@ -72,13 +72,13 @@ shadow_map_pass::shadow_map_pass(gl::rasterizer* rasterizer, const gl::framebuff skinned_model_view_projection_input = skinned_shader_program->get_input("model_view_projection"); // Calculate bias-tile matrices - float4x4 bias_matrix = math::translate(math::identity4x4, float3{0.5f, 0.5f, 0.5f}) * math::scale(math::identity4x4, float3{0.5f, 0.5f, 0.5f}); - float4x4 tile_scale = math::scale(math::identity4x4, float3{0.5f, 0.5f, 1.0f}); + float4x4 bias_matrix = math::translate(math::matrix4::identity, float3{0.5f, 0.5f, 0.5f}) * math::scale(math::matrix4::identity, float3{0.5f, 0.5f, 0.5f}); + float4x4 tile_scale = math::scale(math::matrix4::identity, float3{0.5f, 0.5f, 1.0f}); for (int i = 0; i < 4; ++i) { float x = static_cast(i % 2) * 0.5f; float y = static_cast(i / 2) * 0.5f; - float4x4 tile_matrix = math::translate(math::identity4x4, float3{x, y, 0.0f}) * tile_scale; + float4x4 tile_matrix = math::translate(math::matrix4::identity, float3{x, y, 0.0f}) * tile_scale; bias_tile_matrices[i] = tile_matrix * bias_matrix; } } @@ -210,7 +210,7 @@ void shadow_map_pass::render(const render::context& ctx, render::queue& queue) c offset.y = std::ceil(offset.y * half_shadow_map_resolution) / half_shadow_map_resolution; // Crop the light view-projection matrix - crop_matrix = math::translate(math::identity4x4, offset) * math::scale(math::identity4x4, scale); + crop_matrix = math::translate(math::matrix4::identity, offset) * math::scale(math::matrix4::identity, scale); cropped_view_projection = crop_matrix * light_view_projection; // Calculate shadow matrix diff --git a/src/render/passes/sky-pass.cpp b/src/render/passes/sky-pass.cpp index 3663656..80873fb 100644 --- a/src/render/passes/sky-pass.cpp +++ b/src/render/passes/sky-pass.cpp @@ -38,7 +38,6 @@ #include "scene/camera.hpp" #include "utility/fundamental-types.hpp" #include "color/color.hpp" -#include "astro/illuminance.hpp" #include "math/interpolation.hpp" #include "geom/cartesian.hpp" #include "geom/spherical.hpp" @@ -71,12 +70,12 @@ sky_pass::sky_pass(gl::rasterizer* rasterizer, const gl::framebuffer* framebuffe cloud_shader_program(nullptr), observer_altitude_tween(0.0f, math::lerp), sun_position_tween(float3{1.0f, 0.0f, 0.0f}, math::lerp), - sun_illuminance_outer_tween(float3{1.0f, 1.0f, 1.0f}, math::lerp), - sun_illuminance_inner_tween(float3{1.0f, 1.0f, 1.0f}, math::lerp), + sun_luminance_tween(float3{0.0f, 0.0f, 0.0f}, math::lerp), + sun_illuminance_tween(float3{0.0f, 0.0f, 0.0f}, math::lerp), icrf_to_eus_translation({0, 0, 0}, math::lerp), - icrf_to_eus_rotation(math::quaternion::identity(), math::nlerp), + icrf_to_eus_rotation(math::quaternion::identity, math::nlerp), moon_position_tween(float3{0, 0, 0}, math::lerp), - moon_rotation_tween(math::quaternion::identity(), math::nlerp), + moon_rotation_tween(math::quaternion::identity, math::nlerp), moon_angular_radius_tween(0.0f, math::lerp), moon_sunlight_direction_tween(float3{0, 0, 0}, math::lerp), moon_sunlight_illuminance_tween(float3{0, 0, 0}, math::lerp), @@ -107,7 +106,7 @@ void sky_pass::render(const render::context& ctx, render::queue& queue) const float clip_near = camera.get_clip_near_tween().interpolate(ctx.alpha); float clip_far = camera.get_clip_far_tween().interpolate(ctx.alpha); float3 model_scale = float3{1.0f, 1.0f, 1.0f} * (clip_near + clip_far) * 0.5f; - float4x4 model = math::scale(math::identity4x4, model_scale); + float4x4 model = math::scale(math::matrix4::identity, model_scale); float4x4 view = math::resize<4, 4>(math::resize<3, 3>(camera.get_view_tween().interpolate(ctx.alpha))); float4x4 model_view = view * model; float4x4 projection = camera.get_projection_tween().interpolate(ctx.alpha); @@ -128,9 +127,9 @@ void sky_pass::render(const render::context& ctx, render::queue& queue) const float3 sun_position = sun_position_tween.interpolate(ctx.alpha); float3 sun_direction = math::normalize(sun_position); - // Interpolate sun color - float3 sun_illuminance_outer = sun_illuminance_outer_tween.interpolate(ctx.alpha); - float3 sun_illuminance_inner = sun_illuminance_inner_tween.interpolate(ctx.alpha); + // Interpolate and expose sun luminance and illuminance + float3 sun_luminance = sun_luminance_tween.interpolate(ctx.alpha) * ctx.exposure; + float3 sun_illuminance = sun_illuminance_tween.interpolate(ctx.alpha) * ctx.exposure; // Draw atmosphere if (sky_model) @@ -156,9 +155,10 @@ void sky_pass::render(const render::context& ctx, render::queue& queue) const if (sun_angular_radius_input) sun_angular_radius_input->upload(sun_angular_radius * magnification); - // Pre-exposure sun color + if (sun_luminance_input) + sun_luminance_input->upload(sun_luminance); if (sun_illuminance_input) - sun_illuminance_input->upload(sun_illuminance_outer * ctx.exposure); + sun_illuminance_input->upload(sun_illuminance); if (scale_height_rm_input) scale_height_rm_input->upload(scale_height_rm); @@ -186,7 +186,7 @@ void sky_pass::render(const render::context& ctx, render::queue& queue) const if (cloud_sun_direction_input) cloud_sun_direction_input->upload(sun_direction); if (cloud_sun_illuminance_input) - cloud_sun_illuminance_input->upload(sun_illuminance_inner); + cloud_sun_illuminance_input->upload(sun_illuminance); if (cloud_camera_position_input) cloud_camera_position_input->upload(ctx.camera_transform.translation); if (cloud_camera_exposure_input) @@ -253,10 +253,10 @@ void sky_pass::render(const render::context& ctx, render::queue& queue) const moon_camera_position_input->upload(ctx.camera_transform.translation); if (moon_sunlight_direction_input) moon_sunlight_direction_input->upload(math::normalize(moon_sunlight_direction_tween.interpolate(ctx.alpha))); - if (moon_planetlight_direction_input) - moon_planetlight_direction_input->upload(math::normalize(moon_planetlight_direction_tween.interpolate(ctx.alpha))); if (moon_sunlight_illuminance_input) moon_sunlight_illuminance_input->upload(moon_sunlight_illuminance_tween.interpolate(ctx.alpha) * ctx.exposure); + if (moon_planetlight_direction_input) + moon_planetlight_direction_input->upload(math::normalize(moon_planetlight_direction_tween.interpolate(ctx.alpha))); if (moon_planetlight_illuminance_input) moon_planetlight_illuminance_input->upload(moon_planetlight_illuminance_tween.interpolate(ctx.alpha) * ctx.exposure); @@ -296,6 +296,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_luminance_input = sky_shader_program->get_input("sun_luminance"); 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"); @@ -429,8 +430,8 @@ void sky_pass::update_tweens() { observer_altitude_tween.update(); sun_position_tween.update(); - sun_illuminance_outer_tween.update(); - sun_illuminance_inner_tween.update(); + sun_luminance_tween.update(); + sun_illuminance_tween.update(); icrf_to_eus_translation.update(); icrf_to_eus_rotation.update(); @@ -459,10 +460,14 @@ void sky_pass::set_sun_position(const float3& position) sun_position_tween[1] = position; } -void sky_pass::set_sun_illuminance(const float3& illuminance_outer, const float3& illuminance_inner) +void sky_pass::set_sun_illuminance(const float3& illuminance) +{ + sun_illuminance_tween[1] = illuminance; +} + +void sky_pass::set_sun_luminance(const float3& luminance) { - sun_illuminance_outer_tween[1] = illuminance_outer; - sun_illuminance_inner_tween[1] = illuminance_inner; + sun_luminance_tween[1] = luminance; } void sky_pass::set_sun_angular_radius(float radius) diff --git a/src/render/passes/sky-pass.hpp b/src/render/passes/sky-pass.hpp index 8072087..7c68277 100644 --- a/src/render/passes/sky-pass.hpp +++ b/src/render/passes/sky-pass.hpp @@ -65,7 +65,8 @@ public: void set_icrf_to_eus(const math::transformation::se3& transformation); void set_sun_position(const float3& position); - void set_sun_illuminance(const float3& illuminance_outer, const float3& illuminance_inner); + void set_sun_luminance(const float3& luminance); + void set_sun_illuminance(const float3& illuminance); void set_sun_angular_radius(float radius); void set_observer_altitude(float altitude); void set_scale_heights(float rayleigh, float mie); @@ -92,6 +93,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_luminance_input; const gl::shader_input* sun_illuminance_input; const gl::shader_input* sun_angular_radius_input; const gl::shader_input* scale_height_rm_input; @@ -157,8 +159,8 @@ private: tween observer_altitude_tween; tween sun_position_tween; - tween sun_illuminance_outer_tween; - tween sun_illuminance_inner_tween; + tween sun_luminance_tween; + tween sun_illuminance_tween; tween icrf_to_eus_translation; tween> icrf_to_eus_rotation; diff --git a/src/resources/entity-archetype-loader.cpp b/src/resources/entity-archetype-loader.cpp index f5bb4f0..49054c3 100644 --- a/src/resources/entity-archetype-loader.cpp +++ b/src/resources/entity-archetype-loader.cpp @@ -39,13 +39,6 @@ static bool load_component_atmosphere(entity::archetype& archetype, const json& element) { entity::component::atmosphere component; - component.exosphere_altitude = 0.0; - component.index_of_refraction = 0.0; - component.rayleigh_density = 0.0; - component.mie_density = 0.0; - component.rayleigh_scale_height = 0.0; - component.mie_scale_height = 0.0; - component.mie_anisotropy = 0.0; if (element.contains("exosphere_altitude")) component.exosphere_altitude = element["exosphere_altitude"].get(); @@ -61,6 +54,8 @@ static bool load_component_atmosphere(entity::archetype& archetype, const json& component.mie_scale_height = element["mie_scale_height"].get(); if (element.contains("mie_anisotropy")) component.mie_anisotropy = element["mie_anisotropy"].get(); + if (element.contains("airglow")) + component.airglow = element["airglow"].get(); archetype.set(component); @@ -98,26 +93,32 @@ static bool load_component_blackbody(entity::archetype& archetype, const json& e static bool load_component_celestial_body(entity::archetype& archetype, const json& element) { entity::component::celestial_body component; - component.radius = 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; - component.albedo = 0.0; if (element.contains("radius")) component.radius = element["radius"].get(); if (element.contains("mass")) component.mass = element["mass"].get(); if (element.contains("pole_ra")) - component.pole_ra = math::radians(element["pole_ra"].get()); + { + component.pole_ra.clear(); + auto& pole_ra_element = element["pole_ra"]; + for (auto it = pole_ra_element.rbegin(); it != pole_ra_element.rend(); ++it) + component.pole_ra.push_back(math::radians(it->get())); + } if (element.contains("pole_dec")) - component.pole_dec = math::radians(element["pole_dec"].get()); + { + component.pole_dec.clear(); + auto& pole_dec_element = element["pole_dec"]; + for (auto it = pole_dec_element.rbegin(); it != pole_dec_element.rend(); ++it) + component.pole_dec.push_back(math::radians(it->get())); + } if (element.contains("prime_meridian")) - component.prime_meridian = math::radians(element["prime_meridian"].get()); - if (element.contains("rotation_period")) - component.rotation_period = element["rotation_period"].get(); + { + component.prime_meridian.clear(); + auto& prime_meridian_element = element["prime_meridian"]; + for (auto it = prime_meridian_element.rbegin(); it != prime_meridian_element.rend(); ++it) + component.prime_meridian.push_back(math::radians(it->get())); + } if (element.contains("albedo")) component.albedo = element["albedo"].get(); @@ -176,28 +177,14 @@ static bool load_component_orbit(entity::archetype& archetype, const json& eleme entity::component::orbit component; component.parent = entt::null; - component.elements.ec = 0.0; - component.elements.a = 0.0; - component.elements.in = 0.0; - component.elements.om = 0.0; - component.elements.w = 0.0; - component.elements.ma = 0.0; - component.mean_motion = 0.0; + component.ephemeris_index = -1; + component.scale = 1.0; + component.position = {0, 0, 0}; - if (element.contains("ec")) - component.elements.ec = element["ec"].get(); - if (element.contains("a")) - component.elements.a = element["a"].get(); - if (element.contains("in")) - component.elements.in = math::radians(element["in"].get()); - if (element.contains("om")) - component.elements.om = math::radians(element["om"].get()); - if (element.contains("w")) - component.elements.w = math::radians(element["w"].get()); - if (element.contains("ma")) - component.elements.ma = math::radians(element["ma"].get()); - if (element.contains("n")) - component.mean_motion = math::radians(element["n"].get()); + if (element.contains("ephemeris_index")) + component.ephemeris_index = element["ephemeris_index"].get(); + if (element.contains("scale")) + component.scale = element["scale"].get(); archetype.set(component); @@ -207,7 +194,7 @@ static bool load_component_orbit(entity::archetype& archetype, const json& eleme static bool load_component_transform(entity::archetype& archetype, const json& element) { entity::component::transform component; - component.local = math::identity_transform; + component.local = math::transform::identity; component.warp = true; if (element.contains("translation")) diff --git a/src/resources/ephemeris-loader.cpp b/src/resources/ephemeris-loader.cpp new file mode 100644 index 0000000..b19710e --- /dev/null +++ b/src/resources/ephemeris-loader.cpp @@ -0,0 +1,244 @@ +/* + * 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 . + */ + +#include "resource-loader.hpp" +#include "physics/orbit/ephemeris.hpp" +#include "utility/bit-math.hpp" +#include +#include + +/// Offset to time data in the JPL DE header, in bytes. +static constexpr std::size_t jpl_de_offset_time = 0xA5C; + +/// Offset to the first coefficient table in the JPL DE header, in bytes. +static constexpr std::size_t jpl_de_offset_table1 = 0xA88; + +/// Offset to the DE version number in the JPL DE header, in bytes. +static constexpr std::size_t jpl_de_offset_denum = 0xB18; + +/// Offset to the second coefficient table in the JPL DE header, in bytes. +static constexpr std::size_t jpl_de_offset_table2 = 0xB1C; + +/// Offset to the third coefficient table in the JPL DE header, in bytes, if the constant limit has not been exceeded. +static constexpr std::size_t jpl_de_offset_table3 = 0xB28; + +/// Mask to detect bytes in the most significant word of the JPL DE version number. +static constexpr std::int32_t jpl_de_denum_endian_mask = 0xFFFF0000; + +/// Number of items in the first coefficient table. +static constexpr std::size_t jpl_de_table1_count = 12; + +/// Number of items in the second coefficient table. +static constexpr std::size_t jpl_de_table2_count = 1; + +/// Number of items in the third coefficient table. +static constexpr std::size_t jpl_de_table3_count = 2; + +/// Maximum number of items in a JPL DE file. +static constexpr std::size_t jpl_de_max_item_count = jpl_de_table1_count + jpl_de_table2_count + jpl_de_table3_count; + +/// Maximum number of constants in the first set of constant names. +static constexpr std::size_t jpl_de_constant_limit = 400; + +/// Length of a constant name, in bytes. +static constexpr std::size_t jpl_de_constant_length = 6; + +/// Enumerated IDs of the JPL DE items. +enum +{ + /// Mercury + jpl_de_id_mercury, + + /// Venus + jpl_de_id_venus, + + /// Earth-Moon barycenter + jpl_de_id_embary, + + /// Mars + jpl_de_id_mars, + + /// Jupiter + jpl_de_id_jupiter, + + /// Saturn + jpl_de_id_saturn, + + /// Uranus + jpl_de_id_uranus, + + /// Neptune + jpl_de_id_neptune, + + /// Pluto + jpl_de_id_pluto, + + /// Moon + jpl_de_id_moon, + + /// Sun + jpl_de_id_sun, + + /// Earth nutation + jpl_de_id_earth_nutation, + + /// Lunar mantle libration + jpl_de_id_luma_libration, + + /// Lunar mantle angular velocity + jpl_de_id_luma_angular_velocity, + + /// TT-TDB + jpl_de_id_tt_tdb +}; + +/// Number of components for each JPL DE item. +static constexpr std::size_t jpl_de_component_count[jpl_de_max_item_count] = +{ + 3, // Mercury: x,y,z (km) + 3, // Venus: x,y,z (km) + 3, // Earth-Moon barycenter: x,y,z (km) + 3, // Mars: x,y,z (km) + 3, // Jupiter: x,y,z (km) + 3, // Saturn: x,y,z (km) + 3, // Uranus: x,y,z (km) + 3, // Neptune: x,y,z (km) + 3, // Pluto: x,y,z (km) + 3, // Moon: x,y,z (km) + 3, // Sun: x,y,z (km) + 2, // Earth nutation: d_psi,d_epsilon (radians) + 3, // Lunar mantle libration: phi,theta,ps (radians) + 3, // Lunar mantle angular velocity: omega_x,omega_y,omega_z (radians/day) + 1 // TT-TDB: t (seconds) +}; + +/// Reads and swaps the byte order of 32-bit numbers. +static PHYSFS_sint64 read_swap32(PHYSFS_File* handle, void* buffer, PHYSFS_uint64 len) +{ + PHYSFS_sint64 status = PHYSFS_readBytes(handle, buffer, len); + for (std::uint32_t* ptr32 = (uint32_t*)buffer; len >= 8; len -= 8, ++ptr32) + *ptr32 = bit::swap32(*ptr32); + return status; +} + +/// Reads and swaps the byte order of 64-bit numbers. +static PHYSFS_sint64 read_swap64(PHYSFS_File* handle, void* buffer, PHYSFS_uint64 len) +{ + PHYSFS_sint64 status = PHYSFS_readBytes(handle, buffer, len); + for (std::uint64_t* ptr64 = (uint64_t*)buffer; len >= 8; len -= 8, ++ptr64) + *ptr64 = bit::swap64(*ptr64); + return status; +} + +template <> +physics::orbit::ephemeris* resource_loader>::load(resource_manager* resource_manager, PHYSFS_File* file, const std::filesystem::path& path) +{ + // Init file reading function pointers + PHYSFS_sint64 (*read32)(PHYSFS_File*, void*, PHYSFS_uint64) = &PHYSFS_readBytes; + PHYSFS_sint64 (*read64)(PHYSFS_File*, void*, PHYSFS_uint64) = &PHYSFS_readBytes; + + // Read DE version number + std::int32_t denum; + PHYSFS_seek(file, jpl_de_offset_denum); + PHYSFS_readBytes(file, &denum, sizeof(std::int32_t)); + + // If file endianness does not match host + if (denum & jpl_de_denum_endian_mask) + { + // Use endian-swapping read functions + read32 = &read_swap32; + read64 = &read_swap64; + } + + // Read ephemeris time + double ephemeris_time[3]; + PHYSFS_seek(file, jpl_de_offset_time); + read64(file, ephemeris_time, sizeof(double) * 3); + + // Make time relative to J2000 epoch + const double epoch = 2451545.0; + ephemeris_time[0] -= epoch; + ephemeris_time[1] -= epoch; + + // Read number of constants + std::int32_t constant_count; + read32(file, &constant_count, sizeof(std::int32_t)); + + // Read first coefficient table + std::int32_t coeff_table[jpl_de_max_item_count][3]; + PHYSFS_seek(file, jpl_de_offset_table1); + read32(file, coeff_table, sizeof(std::int32_t) * jpl_de_table1_count * 3); + + // Read second coefficient table + PHYSFS_seek(file, jpl_de_offset_table2); + read32(file, &coeff_table[jpl_de_table1_count][0], sizeof(std::int32_t) * jpl_de_table2_count * 3); + + // Seek past extra constant names + if (constant_count > jpl_de_constant_limit) + PHYSFS_seek(file, jpl_de_offset_table3 + (constant_count - jpl_de_constant_limit) * jpl_de_constant_length); + + // Read third coefficient table + read32(file, &coeff_table[jpl_de_table1_count + jpl_de_table2_count][0], sizeof(std::int32_t) * jpl_de_table3_count * 3); + + // Calculate number of coefficients per record + std::int32_t record_coeff_count = 0; + for (int i = 0; i < jpl_de_max_item_count; ++i) + { + std::int32_t coeff_count = coeff_table[i][0] + coeff_table[i][1] * coeff_table[i][2] * jpl_de_component_count[i] - 1; + record_coeff_count = std::max(record_coeff_count, coeff_count); + } + + // Calculate record size and record count + std::size_t record_size = record_coeff_count * sizeof(double); + std::size_t record_count = static_cast((ephemeris_time[1] - ephemeris_time[0]) / ephemeris_time[2]); + + // Calculate coefficient strides + std::size_t strides[11]; + for (int i = 0; i < 11; ++i) + strides[i] = coeff_table[i][2] * coeff_table[i][1] * 3; + + // Allocate and resize ephemeris to accommodate items 0-10 + physics::orbit::ephemeris* ephemeris = new physics::orbit::ephemeris(); + ephemeris->resize(11); + + // Init trajectories + for (int i = 0; i < 11; ++i) + { + auto& trajectory = (*ephemeris)[i]; + trajectory.t0 = ephemeris_time[0]; + trajectory.t1 = ephemeris_time[1]; + trajectory.dt = ephemeris_time[2] / static_cast(coeff_table[i][2]); + trajectory.n = coeff_table[i][1]; + trajectory.a.resize(record_count * strides[i]); + } + + // Read coefficients + for (std::size_t i = 0; i < record_count; ++i) + { + // Seek to coefficient record + PHYSFS_seek(file, (i + 2) * record_size + 2 * sizeof(double)); + + for (int j = 0; j < 11; ++j) + { + read64(file, &(*ephemeris)[j].a[i * strides[j]], sizeof(double) * strides[j]); + } + } + + return ephemeris; +} diff --git a/src/scene/camera.cpp b/src/scene/camera.cpp index aae6020..d503e14 100644 --- a/src/scene/camera.cpp +++ b/src/scene/camera.cpp @@ -72,9 +72,9 @@ camera::camera(): clip_far(1.0f, math::lerp), fov(math::half_pi, math::lerp), aspect_ratio(1.0f, math::lerp), - view(math::identity4x4, std::bind(&interpolate_view, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3)), - projection(math::identity4x4, std::bind(&interpolate_projection, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3)), - view_projection(math::identity4x4, std::bind(&interpolate_view_projection, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3)), + view(math::matrix4::identity, std::bind(&interpolate_view, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3)), + projection(math::matrix4::identity, std::bind(&interpolate_projection, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3)), + view_projection(math::matrix4::identity, std::bind(&interpolate_view_projection, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3)), exposure(0.0f, math::lerp) {} @@ -149,11 +149,6 @@ void camera::set_exposure(float ev100) exposure[1] = ev100; } -void camera::set_exposure(float f_number, float speed, float iso) -{ - exposure[1] = std::log2((f_number * f_number) / speed * 100.0f / iso); -} - void camera::set_compositor(render::compositor* compositor) { this->compositor = compositor; diff --git a/src/scene/camera.hpp b/src/scene/camera.hpp index 713a01c..8174780 100644 --- a/src/scene/camera.hpp +++ b/src/scene/camera.hpp @@ -83,15 +83,6 @@ public: * @param ev100 ISO 100 exposure value. */ void set_exposure(float ev100); - - /** - * Sets the camera's ISO 100 exposure value given exposure settings. - * - * @param f_number F-number. - * @param speed Shutter speed, in seconds. - * @param iso ISO sensitivity value. - */ - void set_exposure(float f_number, float speed, float iso); void set_compositor(render::compositor* compositor); void set_composite_index(int index); diff --git a/src/scene/object.cpp b/src/scene/object.cpp index 9777985..0aa7c3f 100644 --- a/src/scene/object.cpp +++ b/src/scene/object.cpp @@ -34,7 +34,7 @@ typename object_base::transform_type object_base::interpolate_transforms(const t object_base::object_base(): active(true), - transform(math::identity_transform, interpolate_transforms), + transform(math::transform::identity, interpolate_transforms), culling_mask(nullptr) {} diff --git a/src/utility/bit-math.hpp b/src/utility/bit-math.hpp index 3166345..008beba 100644 --- a/src/utility/bit-math.hpp +++ b/src/utility/bit-math.hpp @@ -20,6 +20,8 @@ #ifndef ANTKEEPER_BIT_MATH_HPP #define ANTKEEPER_BIT_MATH_HPP +#include + /// Bitwise math namespace bit { @@ -169,6 +171,23 @@ constexpr T segregate(T x) noexcept; template constexpr T swap_adjacent(T x) noexcept; +/// Swaps the byte order of an unsigned 16-bit number. +std::uint16_t swap16(std::uint16_t x) noexcept; + +/// Swaps the byte order of a signed 16-bit number. +std::int16_t swap16(std::int16_t x) noexcept; + +/// Swaps the byte order of an unsigned 32-bit number. +std::uint32_t swap32(std::uint32_t x) noexcept; + +/// Swaps the byte order of a signed 32-bit number. +std::int32_t swap32(std::int32_t x) noexcept; + +/// Swaps the byte order of an unsigned 64-bit number. +std::uint64_t swap64(std::uint64_t x) noexcept; + +/// Swaps the byte order of a signed 64-bit number. +std::int64_t swap64(std::int64_t x) noexcept; template constexpr T compress(T x) noexcept @@ -326,6 +345,42 @@ constexpr T swap_adjacent(T x) noexcept return ((x & T(0xaaaaaaaaaaaaaaaa)) >> 1) | ((x & T(0x5555555555555555)) << 1); } +std::uint16_t swap16(std::uint16_t x) noexcept +{ + return (x << 8) | (x >> 8); +} + +std::int16_t swap16(std::int16_t x) noexcept +{ + return (x << 8) | ((x >> 8) & 0xFF); +} + +std::uint32_t swap32(std::uint32_t x) noexcept +{ + x = ((x << 8) & 0xFF00FF00) | ((x >> 8) & 0xFF00FF); + return (x << 16) | (x >> 16); +} + +std::int32_t swap32(std::int32_t x) noexcept +{ + x = ((x << 8) & 0xFF00FF00) | ((x >> 8) & 0xFF00FF); + return (x << 16) | ((x >> 16) & 0xFFFF); +} + +std::uint64_t swap64(std::uint64_t x) noexcept +{ + x = ((x << 8) & 0xFF00FF00FF00FF00) | ((x >> 8) & 0x00FF00FF00FF00FF); + x = ((x << 16) & 0xFFFF0000FFFF0000) | ((x >> 16) & 0x0000FFFF0000FFFF); + return (x << 32) | (x >> 32); +} + +std::int64_t swap64(std::int64_t x) noexcept +{ + x = ((x << 8) & 0xFF00FF00FF00FF00) | ((x >> 8) & 0x00FF00FF00FF00FF); + x = ((x << 16) & 0xFFFF0000FFFF0000) | ((x >> 16) & 0x0000FFFF0000FFFF); + return (x << 32) | ((x >> 32) & 0xFFFFFFFF); +} + } // namespace bit #endif // ANTKEEPER_BIT_MATH_HPP