From 7477552eea617d6f4c845ec39b6a9e93395bbe99 Mon Sep 17 00:00:00 2001 From: "C. J. Howard" Date: Wed, 19 May 2021 23:10:13 +0800 Subject: [PATCH] Add physics namespace, move orbit-related functions into orbit namespace within physics namespace, add physics reference frame class, add functions to construct orbital reference frames --- src/astro/astro.hpp | 2 - src/astro/orbit.cpp | 78 ------- src/astro/orbit.hpp | 75 ------ src/coordinates/ecliptic.hpp | 214 ----------------- src/coordinates/equatorial.hpp | 208 ----------------- src/coordinates/horizontal.hpp | 220 ------------------ src/coordinates/rectangular.hpp | 175 -------------- .../components/celestial-body-component.hpp | 8 +- src/ecs/components/orbit-component.hpp | 8 +- src/ecs/systems/astronomy-system.cpp | 128 ++++++++-- src/ecs/systems/solar-system.cpp | 13 +- src/game/states/play-state.cpp | 18 +- src/geom/cartesian.hpp | 58 +++++ src/geom/geom.hpp | 2 + src/{coordinates => geom}/spherical.hpp | 22 +- src/math/matrix-functions.hpp | 69 ++++++ src/math/quaternion-type.hpp | 54 ++++- src/physics/frame.hpp | 134 +++++++++++ src/physics/orbit/elements.hpp | 95 ++++++++ src/physics/orbit/frames.hpp | 131 +++++++++++ src/physics/orbit/kepler.hpp | 79 +++++++ src/physics/orbit/orbit.hpp | 39 ++++ src/physics/orbit/state.hpp | 52 +++++ .../coordinates.hpp => physics/physics.hpp} | 19 +- src/renderer/passes/sky-pass.cpp | 20 +- 25 files changed, 865 insertions(+), 1056 deletions(-) delete mode 100644 src/astro/orbit.cpp delete mode 100644 src/astro/orbit.hpp delete mode 100644 src/coordinates/ecliptic.hpp delete mode 100644 src/coordinates/equatorial.hpp delete mode 100644 src/coordinates/horizontal.hpp delete mode 100644 src/coordinates/rectangular.hpp create mode 100644 src/geom/cartesian.hpp rename src/{coordinates => geom}/spherical.hpp (72%) create mode 100644 src/physics/frame.hpp create mode 100644 src/physics/orbit/elements.hpp create mode 100644 src/physics/orbit/frames.hpp create mode 100644 src/physics/orbit/kepler.hpp create mode 100644 src/physics/orbit/orbit.hpp create mode 100644 src/physics/orbit/state.hpp rename src/{coordinates/coordinates.hpp => physics/physics.hpp} (68%) diff --git a/src/astro/astro.hpp b/src/astro/astro.hpp index d030ec4..8e6115d 100644 --- a/src/astro/astro.hpp +++ b/src/astro/astro.hpp @@ -26,8 +26,6 @@ namespace astro {} #include "apparent-size.hpp" #include "coordinates.hpp" #include "constants.hpp" -#include "coordinates.hpp" #include "illuminance.hpp" -#include "orbit.hpp" #endif // ANTKEEPER_ASTRO_HPP diff --git a/src/astro/orbit.cpp b/src/astro/orbit.cpp deleted file mode 100644 index c157384..0000000 --- a/src/astro/orbit.cpp +++ /dev/null @@ -1,78 +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 . - */ - -#include "orbit.hpp" -#include "math/angles.hpp" -#include - -namespace astro -{ - -double solve_kepler(double ec, double ma, double tolerance, std::size_t iterations) -{ - // Approximate eccentric anomaly, E - double e0 = ma + ec * std::sin(ma) * (1.0 + ec * std::cos(ma)); - - // Iteratively converge E0 and E1 - for (std::size_t i = 0; i < iterations; ++i) - { - double e1 = e0 - (e0 - ec * std::sin(e0) - ma) / (1.0 - ec * std::cos(e0)); - double error = std::abs(e1 - e0); - e0 = e1; - - if (error < tolerance) - break; - } - - return e0; -} - -double3 orbital_elements_to_ecliptic(const orbital_elements& elements, double ke_tolerance, std::size_t ke_iterations) -{ - // Calculate semi-minor axis, b - double b = elements.a * std::sqrt(1.0 - elements.ec * elements.ec); - - // Solve Kepler's equation for eccentric anomaly, E - double ea = solve_kepler(elements.ec, elements.ma, ke_tolerance, ke_iterations); - - // Calculate radial distance, r; and true anomaly, v - double xv = elements.a * (std::cos(ea) - elements.ec); - double yv = b * std::sin(ea); - double r = std::sqrt(xv * xv + yv * yv); - double v = std::atan2(yv, xv); - - // Calculate true longitude, l - double l = elements.w + v; - - // Transform vector (r, 0, 0) from local coordinates to ecliptic coordinates - // = Rz(-omega) * Rx(-i) * Rz(-l) * r - double cos_om = std::cos(elements.om); - double sin_om = std::sin(elements.om); - double cos_i = std::cos(elements.i); - double cos_l = std::cos(l); - double sin_l = std::sin(l); - return double3 - { - r * (cos_om * cos_l - sin_om * sin_l * cos_i), - r * (sin_om * cos_l + cos_om * sin_l * cos_i), - r * sin_l * std::sin(elements.i) - }; -} - -} // namespace astro diff --git a/src/astro/orbit.hpp b/src/astro/orbit.hpp deleted file mode 100644 index 594a8a5..0000000 --- a/src/astro/orbit.hpp +++ /dev/null @@ -1,75 +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_ORBIT_HPP -#define ANTKEEPER_ASTRO_ORBIT_HPP - -#include "utility/fundamental-types.hpp" - -namespace astro -{ - -/** - * Contains six orbital elements which describe a Keplerian orbit. - */ -struct orbital_elements -{ - double ec; ///< Eccentricity, e - double a; ///< Semi-major axis, a - double i; ///< Inclination, i (radians) - double om; ///< Longitude of the ascending node, OMEGA (radians) - double w; ///< Argument of periapsis, w (radians) - double ma; ///< Mean anomaly, M (radians) -}; - -/** - * Orbital state vectors. - */ -struct orbital_state -{ - double3 r; ///< Cartesian position, r. - double3 v; ///< Cartesian velocity, v. -}; - -/** - * Iteratively solves Kepler's equation for eccentric anomaly, E. - * - * @param ec Eccentricity, e. - * @param ma Mean anomaly, M (radians). - * @param tolerance Tolerance of solution. - * @param iterations Maximum number of iterations. - * @return Eccentric anomaly. - */ -double solve_kepler(double ec, double ma, double tolerance, std::size_t iterations); - -/** - * Calculates the ecliptic rectangular orbital position from Keplerian orbital elements. - * - * @note Only works for elliptical orbits. - * - * @param elements Orbital elements. - * @param ke_tolerance Kepler's equation tolerance. - * @param ke_iterations Kepler's equation iterations. - * @return Orbital state. - */ -double3 orbital_elements_to_ecliptic(const orbital_elements& elements, double ke_tolerance, std::size_t ke_iterations); - -} // namespace astro - -#endif // ANTKEEPER_ASTRO_ORBIT_HPP diff --git a/src/coordinates/ecliptic.hpp b/src/coordinates/ecliptic.hpp deleted file mode 100644 index 0c2cba0..0000000 --- a/src/coordinates/ecliptic.hpp +++ /dev/null @@ -1,214 +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_COORDINATES_ECLIPTIC_HPP -#define ANTKEEPER_COORDINATES_ECLIPTIC_HPP - -#include "coordinates/rectangular.hpp" -#include "coordinates/spherical.hpp" -#include "utility/fundamental-types.hpp" -#include - -namespace coordinates { - -namespace rectangular { - - /// Rectangular coordinate system with the plane of the Earth's orbit as the fundamental plane. This is a right-handed coordinate system with the x-axis pointing to the vernal equinox, the y-axis pointing east, and the z-axis is the north orbital pole. - namespace ecliptic { - - /** - * Produces a matrix which rotates rectangular coordinates from ecliptic space into equatorial space. - * - * @param ecl Obliquity of the ecliptic, in radians. - * @return Rotation matrix. - * - * @see coordinates::rectangular::equatorial - */ - template - math::matrix3 to_equatorial(T ecl); - - /** - * Rotates rectangular coordinates from ecliptic space into equatorial space. - * - * @param v Rectangular coordinates in ecliptic space. - * @param ecl Obliquity of the ecliptic, in radians. - * @return Rectangular coordinates in equatorial space. - * - * @see coordinates::rectangular::equatorial - */ - template - math::vector3 to_equatorial(const math::vector3& v, T ecl); - - /** - * Produces a matrix which rotates rectangular coordinates from ecliptic space into local horizontal space. - * - * @param ecl Obliquity of the ecliptic, in radians. - * @param lat Observer's latitude, in radians. - * @param lst Local sidereal time, in radians. - * @return Rotation matrix. - * - * @see coordinates::rectangular::horizontal - */ - template - math::matrix3 to_horizontal(T ecl, T lat, T lst); - - /** - * Rotates rectangular coordinates from ecliptic space into local horizontal space. - * - * @param v Rectangular coordinates in ecliptic space. - * @param ecl Obliquity of the ecliptic, in radians. - * @param lat Observer's latitude, in radians. - * @param lst Local sidereal time, in radians. - * @return Rectangular coordinates in local horizontal space. - * - * @see coordinates::rectangular::horizontal - */ - template - math::vector3 to_horizontal(const math::vector3& v, T ecl, T lat, T lst); - - } // namespace ecliptic - -} // namespace rectangular - -namespace spherical { - - /// Spherical ecliptic coordinate system. - namespace ecliptic { - - /** - * Rotates spherical coordinates from ecliptic space into equatorial space. - * - * @param v Spherical coordinates in ecliptic space, in the ISO order of radial distance, ecliptic latitude (radians), and ecliptic longitude (radians). - * @param ecl Obliquity of the ecliptic, in radians. - * @return Spherical coordinates in equatorial space, in the ISO order of radial distance, declination (radians), and right ascension (radians). - * - * @see coordinates::spherical::equatorial - */ - template - math::vector3 to_equatorial(const math::vector3& v, T ecl); - - /** - * Rotates spherical coordinates from ecliptic space into local horizontal space. - * - * @param v Spherical coordinates in ecliptic space, in the ISO order of radial distance, ecliptic latitude (radians), and ecliptic longitude (radians). - * @param ecl Obliquity of the ecliptic, in radians. - * @param lat Observer's latitude, in radians. - * @param lst Local sidereal time, in radians. - * @return Spherical coordinates in local horizontal space, in the ISO order of radial distance, altitude (radians), and azimuth (radians). - * - * @see coordinates::spherical::horizontal - */ - template - math::vector3 to_horizontal(const math::vector3& v, T ecl, T lat, T lst); - - } // namespace ecliptic - -} // namespace spherical - -namespace rectangular { - - namespace ecliptic { - - template - math::matrix3 to_equatorial(T ecl) - { - const T c_ecl = std::cos(ecl); - const T s_ecl = std::sin(ecl); - - return math::matrix3 - { - T(1.0), T(0.0), T(0.0), - T(0.0), c_ecl, s_ecl, - T(0.0), -s_ecl, c_ecl - }; - } - - template - math::vector3 to_equatorial(const math::vector3& v, T ecl) - { - return to_equatorial(ecl) * v; - } - - template - math::matrix3 to_horizontal(T ecl, T lat, T lst) - { - const T c_ecl = std::cos(ecl); - const T s_ecl = std::sin(ecl); - const T c_lat = std::cos(lat); - const T s_lat = std::sin(lat); - const T c_lst = std::cos(lst); - const T s_lst = std::sin(lst); - - return math::matrix3 - { - s_lat * c_lst, s_lst, c_lat * c_lst, - s_lat * s_lst * c_ecl - c_lat * s_ecl, -c_lst * c_ecl, c_lat * s_lst * c_ecl + s_lat * s_ecl, - s_lat * s_lst * -s_ecl - c_lat * c_ecl, c_lst * s_ecl, c_lat * s_lst * -s_ecl + s_lat * c_ecl - }; - } - - template - math::vector3 to_horizontal(const math::vector3& v, T ecl, T lat, T lst) - { - return to_horizontal(ecl, lat, lst) * v; - } - - } // namespace ecliptic - -} // namespace rectangular - -namespace spherical { - - namespace ecliptic { - - template - math::vector3 to_equatorial(const math::vector3& v, T ecl) - { - return coordinates::rectangular::to_spherical - ( - coordinates::rectangular::ecliptic::to_equatorial - ( - coordinates::spherical::to_rectangular(v), - ecl - ) - ); - } - - template - math::vector3 to_horizontal(const math::vector3& v, T ecl, T lat, T lst) - { - return coordinates::rectangular::to_spherical - ( - coordinates::rectangular::ecliptic::to_horizontal - ( - coordinates::spherical::to_rectangular(v), - ecl, - lat, - lst - ) - ); - } - - } // namespace ecliptic - -} // namespace spherical - -} // namespace coordinates - -#endif // ANTKEEPER_COORDINATES_ECLIPTIC_HPP diff --git a/src/coordinates/equatorial.hpp b/src/coordinates/equatorial.hpp deleted file mode 100644 index cf0c8ac..0000000 --- a/src/coordinates/equatorial.hpp +++ /dev/null @@ -1,208 +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_COORDINATES_EQUATORIAL_HPP -#define ANTKEEPER_COORDINATES_EQUATORIAL_HPP - -#include "coordinates/rectangular.hpp" -#include "coordinates/spherical.hpp" -#include "utility/fundamental-types.hpp" -#include - -namespace coordinates { - -namespace rectangular { - - /// Rectangular coordinate system with the Earth's equator as the fundamental plane. This is a right-handed coordinate system with the x-axis pointing to the vernal equinox, the y-axis pointing east, and the z-axis is the north celestial pole. - namespace equatorial { - - /** - * Produces a matrix which rotates rectangular coordinates from equatorial space into ecliptic space. - * - * @param ecl Obliquity of the ecliptic, in radians. - * @return Rotation matrix. - * - * @see coordinates::rectangular::ecliptic - */ - template - math::matrix3 to_ecliptic(T ecl); - - /** - * Rotates rectangular coordinates from equatorial space into ecliptic space. - * - * @param v Rectangular coordinates in equatorial space. - * @param ecl Obliquity of the ecliptic, in radians. - * @return Rectangular coordinates in ecliptic space. - * - * @see coordinates::rectangular::ecliptic - */ - template - math::vector3 to_ecliptic(const math::vector3& v, T ecl); - - /** - * Produces a matrix which rotates rectangular coordinates from equatorial space into local horizontal space. - * - * @param lat Observer's latitude, in radians. - * @param lst Local sidereal time, in radians. - * @return Rotation matrix. - * - * @see coordinates::rectangular::horizontal - */ - template - math::matrix3 to_horizontal(T lat, T lst); - - /** - * Rotates rectangular coordinates from equatorial space into local horizontal space. - * - * @param v Rectangular coordinates in equatorial space. - * @param lat Observer's latitude, in radians. - * @param lst Local sidereal time, in radians. - * @return Rectangular coordinates in local horizontal space. - * - * @see coordinates::rectangular::horizontal - */ - template - math::vector3 to_horizontal(const math::vector3& v, T lat, T lst); - - } // namespace equatorial - -} // namespace rectangular - -namespace spherical { - - /// Spherical equatorial coordinate system. - namespace equatorial { - - /** - * Rotates spherical coordinates from equatorial space into ecliptic space. - * - * @param v Spherical coordinates in equatorial space, in the ISO order of radial distance, declination (radians), and right ascension (radians). - * @param ecl Obliquity of the ecliptic, in radians. - * @return Spherical coordinates in ecliptic space, in the ISO order of radial distance, ecliptic latitude (radians), and ecliptic longitude (radians). - * - * @see coordinates::spherical::ecliptic - */ - template - math::vector3 to_ecliptic(const math::vector3& v, T ecl); - - /** - * Rotates spherical coordinates from equatorial space into local horizontal space. - * - * @param v Spherical coordinates in equatorial space, in the ISO order of radial distance, declination (radians), and right ascension (radians). - * @param lat Observer's latitude, in radians. - * @param lst Local sidereal time, in radians. - * @return Spherical coordinates in local horizontal space, in the ISO order of radial distance, altitude (radians), and azimuth (radians). - * - * @see coordinates::spherical::horizontal - */ - template - math::vector3 to_horizontal(const math::vector3& v, T lat, T lst); - - } // namespace equatorial - -} // namespace spherical - -namespace rectangular { - - namespace equatorial { - - template - math::matrix3 to_ecliptic(T ecl) - { - const T c_ecl = std::cos(ecl); - const T s_ecl = std::sin(ecl); - - return math::matrix3 - { - T(1.0), T(0.0), T(0.0), - T(0.0), c_ecl, -s_ecl, - T(0.0), s_ecl, c_ecl - }; - } - - template - math::vector3 to_ecliptic(const math::vector3& v, T ecl) - { - return to_ecliptic(ecl) * v; - } - - template - math::matrix3 to_horizontal(T lat, T lst) - { - const T c_lat = std::cos(lat); - const T s_lat = std::sin(lat); - const T c_lst = std::cos(lst); - const T s_lst = std::sin(lst); - - return math::matrix3 - { - s_lat * c_lst, s_lst, c_lat * c_lst, - s_lat * s_lst, -c_lst, c_lat * s_lst, - -c_lat, T(0.0), s_lat - }; - } - - template - math::vector3 to_horizontal(const math::vector3& v, T lat, T lst) - { - return to_horizontal(lat, lst) * v; - } - - } // namespace equatorial - -} // namespace rectangular - -namespace spherical { - - namespace equatorial { - - template - math::vector3 to_ecliptic(const math::vector3& v, T ecl) - { - return coordinates::rectangular::to_spherical - ( - coordinates::rectangular::equatorial::to_ecliptic - ( - coordinates::spherical::to_rectangular(v), - ecl - ) - ); - } - - template - math::vector3 to_horizontal(const math::vector3& v, T lat, T lst) - { - return coordinates::rectangular::to_spherical - ( - coordinates::rectangular::equatorial::to_horizontal - ( - coordinates::spherical::to_rectangular(v), - lat, - lst - ) - ); - } - - } // namepace equatorial - -} // namespace spherical - -} // namespace coordinates - -#endif // ANTKEEPER_COORDINATES_EQUATORIAL_HPP diff --git a/src/coordinates/horizontal.hpp b/src/coordinates/horizontal.hpp deleted file mode 100644 index d64e648..0000000 --- a/src/coordinates/horizontal.hpp +++ /dev/null @@ -1,220 +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_COORDINATES_HORIZONTAL_HPP -#define ANTKEEPER_COORDINATES_HORIZONTAL_HPP - -#include "coordinates/rectangular.hpp" -#include "coordinates/spherical.hpp" -#include "utility/fundamental-types.hpp" -#include - -namespace coordinates { - -namespace rectangular { - - /// Rectangular local horizontal coordinate system in which the x-axis points north, the y-axis points east, and the z-axis points to the vertical. - namespace horizontal { - - /** - * Produces a matrix which rotates rectangular coordinates from local horizontal space into equatorial space. - * - * @param lat Observer's latitude, in radians. - * @param lst Local sidereal time, in radians. - * @return Rotation matrix. - * - * @see coordinates::rectangular::equatorial - */ - template - math::matrix3 to_equatorial(T lat, T lst); - - /** - * Rotates rectangular coordinates from local horizontal space into equatorial space. - * - * @param v Rectangular coordinates in local horizontal space. - * @param lat Observer's latitude, in radians. - * @param lst Local sidereal time, in radians. - * @return Rectangular coordinates in equatorial space. - * - * @see coordinates::rectangular::equatorial - */ - template - math::vector3 to_equatorial(const math::vector3& v, T lat, T lst); - - /** - * Produces a matrix which rotates rectangular coordinates from local horizontal space into ecliptic space. - * - * @param ecl Obliquity of the ecliptic, in radians. - * @param lat Observer's latitude, in radians. - * @param lst Local sidereal time, in radians. - * @return Rotation matrix. - * - * @see coordinates::rectangular::ecliptic - */ - template - math::matrix3 to_ecliptic(T ecl, T lat, T lst); - - /** - * Rotates rectangular coordinates from local horizontal space into ecliptic space. - * - * @param v Rectangular coordinates in local horizontal space. - * @param ecl Obliquity of the ecliptic, in radians. - * @param lat Observer's latitude, in radians. - * @param lst Local sidereal time, in radians. - * @return Rectangular coordinates in ecliptic space. - * - * @see coordinates::rectangular::ecliptic - */ - template - math::vector3 to_ecliptic(const math::vector3& v, T ecl, T lat, T lst); - - } // namespace horizontal - -} // namespace rectangular - -namespace spherical { - - /// Spherical local horizontal coordinate system. - namespace horizontal { - - /** - * Rotates spherical coordinates from local horizontal space into equatorial space. - * - * @param v Spherical coordinates in local horizontal space, in the ISO order of radial distance, altitude (radians), and azimuth (radians). - * @param lat Observer's latitude, in radians. - * @param lst Local sidereal time, in radians. - * @return Spherical coordinates in equatorial space, in the ISO order of radial distance, declination (radians), and right ascension (radians). - * - * @see coordinates::spherical::equatorial - */ - template - math::vector3 to_equatorial(const math::vector3& v, T lat, T lst); - - /** - * Rotates spherical coordinates from local horizontal space into ecliptic space. - * - * @param v Spherical coordinates in local horizontal space, in the ISO order of radial distance, altitude (radians), and azimuth (radians). - * @param ecl Obliquity of the ecliptic, in radians. - * @param lat Observer's latitude, in radians. - * @param lst Local sidereal time, in radians. - * @return Spherical coordinates in ecliptic space, in the ISO order of radial distance, ecliptic latitude (radians), and ecliptic longitude (radians). - * - * @see coordinates::spherical::ecliptic - */ - template - math::vector3 to_ecliptic(const math::vector3& v, T ecl, T lat, T lst); - - } // namespace horizontal - -} // namespace spherical - -namespace rectangular { - - namespace horizontal { - - template - math::matrix3 to_equatorial(T lat, T lst) - { - const T c_lat = std::cos(lat); - const T s_lat = std::sin(lat); - const T c_lst = std::cos(lst); - const T s_lst = std::sin(lst); - - return math::matrix3 - { - c_lst * s_lat, s_lst * s_lat, -c_lat, - s_lst, c_lst, T(0.0), - c_lst * c_lat, s_lst * c_lat, s_lat - }; - } - - template - math::vector3 to_equatorial(const math::vector3& v, T lat, T lst) - { - return to_equatorial(lat, lst) * v; - } - - template - math::matrix3 to_ecliptic(T ecl, T lat, T lst) - { - const T c_ecl = std::cos(ecl); - const T s_ecl = std::sin(ecl); - const T c_lat = std::cos(lat); - const T s_lat = std::sin(lat); - const T c_lst = std::cos(lst); - const T s_lst = std::sin(lst); - - return math::matrix3 - { - s_lat * c_lst, s_lat * s_lst * c_ecl - c_lat * s_ecl, s_lat * s_lst * -s_ecl - c_lat * c_ecl, - s_lst, -c_lst * c_ecl, c_lst * s_ecl, - c_lat * c_lst, c_lat * s_lst * c_ecl + s_lat * s_ecl, c_lat * s_lst * -s_ecl + s_lat * c_ecl - }; - } - - template - math::vector3 to_ecliptic(const math::vector3& v, T ecl, T lat, T lst) - { - return to_ecliptic(ecl, lat, lst) * v; - } - - } // namespace horizontal - -} // namespace rectangular - -namespace spherical { - - namespace horizontal { - - template - math::vector3 to_equatorial(const math::vector3& v, T lat, T lst) - { - return coordinates::rectangular::to_spherical - ( - coordinates::rectangular::horizontal::to_equatorial - ( - coordinates::spherical::to_rectangular(v), - lat, - lst - ) - ); - } - - template - math::vector3 to_ecliptic(const math::vector3& v, T ecl, T lat, T lst) - { - return coordinates::rectangular::to_spherical - ( - coordinates::rectangular::horizontal::to_ecliptic - ( - coordinates::spherical::to_rectangular(v), - ecl, - lat, - lst - ) - ); - } - - } // namespace horizontal - -} // namespace spherical - -} // namespace coordinates - -#endif // ANTKEEPER_COORDINATES_HORIZONTAL_HPP diff --git a/src/coordinates/rectangular.hpp b/src/coordinates/rectangular.hpp deleted file mode 100644 index 382e537..0000000 --- a/src/coordinates/rectangular.hpp +++ /dev/null @@ -1,175 +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_COORDINATES_RECTANGULAR_HPP -#define ANTKEEPER_COORDINATES_RECTANGULAR_HPP - -#include "utility/fundamental-types.hpp" -#include - -namespace coordinates { - -/// Rectangular (Cartesian) coordinate systems. -namespace rectangular { - -/** - * Produces a matrix which rotates rectangular coordinates about the x-axis by a given angle. - * - * @param angle Angle of rotation, in radians. - * @return Rotation matrix. - */ -template -math::matrix3 rotate_x(T angle); - -/** - * Rotates rectangular coordinates about the x-axis. - * - * @param v Rectangular coordinates to rotate. - * @param angle Angle of rotation, in radians. - * @return Rotated rectangular coordinates. - */ -template -math::vector3 rotate_x(const math::vector3& v, T angle); - -/** - * Produces a matrix which rotates rectangular coordinates about the y-axis by a given angle. - * - * @param angle Angle of rotation, in radians. - * @return Rotation matrix. - */ -template -math::matrix3 rotate_y(T angle); - -/** - * Rotates rectangular coordinates about the y-axis. - * - * @param v Rectangular coordinates to rotate. - * @param angle Angle of rotation, in radians. - * @return Rotated rectangular coordinates. - */ -template -math::vector3 rotate_y(const math::vector3& v, T angle); - -/** - * Produces a matrix which rotates rectangular coordinates about the z-axis by a given angle. - * - * @param angle Angle of rotation, in radians. - * @return Rotation matrix. - */ -template -math::matrix3 rotate_z(T angle); - -/** - * Rotates rectangular coordinates about the z-axis. - * - * @param v Rectangular coordinates to rotate. - * @param angle Angle of rotation, in radians. - * @return Rotated rectangular coordinates. - */ -template -math::vector3 rotate_z(const math::vector3& v, T angle); - -/** - * Converts rectangular coordinates to spherical coordinates. - * - * @param v Rectangular coordinates. - * @return Spherical coordinates, in the ISO order of radial distance, polar angle (radians), and azimuthal angle (radians). - * - * @see coordinates::spherical - */ -template -math::vector3 to_spherical(const math::vector3& v); - -template -math::vector3 to_spherical(const math::vector3& v) -{ - const T xx_yy = v.x * v.x + v.y * v.y; - - return math::vector3 - { - std::sqrt(xx_yy + v.z * v.z), - std::atan2(v.z, std::sqrt(xx_yy)), - std::atan2(v.y, v.x) - }; -} - -template -math::matrix3 rotate_x(T angle) -{ - const T c = std::cos(angle); - const T s = std::sin(angle); - - return math::matrix3 - { - T(1), T(0), T(0), - T(0), c, -s, - T(0), s, c - }; -} - -template -math::vector3 rotate_x(const math::vector3& v, T angle) -{ - return rotate_x(angle) * v; -} - -template -math::matrix3 rotate_y(T angle) -{ - const T c = std::cos(angle); - const T s = std::sin(angle); - - return math::matrix3 - { - c, T(0), s, - T(0), T(1), T(0), - -s, T(0), c - }; -} - -template -math::vector3 rotate_y(const math::vector3& v, T angle) -{ - return rotate_y(angle) * v; -} - -template -math::matrix3 rotate_z(T angle) -{ - const T c = std::cos(angle); - const T s = std::sin(angle); - - return math::matrix3 - { - c, -s, T(0), - s, c, T(0), - T(0), T(0), T(1) - }; -} - -template -math::vector3 rotate_z(const math::vector3& v, T angle) -{ - return rotate_z(angle) * v; -} - -} // namespace rectangular -} // namespace coordinates - -#endif // ANTKEEPER_COORDINATES_RECTANGULAR_HPP diff --git a/src/ecs/components/celestial-body-component.hpp b/src/ecs/components/celestial-body-component.hpp index 3119ce0..51cbc59 100644 --- a/src/ecs/components/celestial-body-component.hpp +++ b/src/ecs/components/celestial-body-component.hpp @@ -20,7 +20,7 @@ #ifndef ANTKEEPER_ECS_CELESTIAL_BODY_COMPONENT_HPP #define ANTKEEPER_ECS_CELESTIAL_BODY_COMPONENT_HPP -#include "astro/orbit.hpp" +#include "physics/orbit/elements.hpp" #include "utility/fundamental-types.hpp" #include "math/quaternion-type.hpp" @@ -28,9 +28,9 @@ namespace ecs { struct celestial_body_component { - astro::orbital_elements orbital_elements; - astro::orbital_elements orbital_rate; - astro::orbital_state orbital_state; + physics::orbit::elements orbital_elements; + physics::orbit::elements orbital_rate; + physics::orbit::state orbital_state; double3 position; double3 velocity; diff --git a/src/ecs/components/orbit-component.hpp b/src/ecs/components/orbit-component.hpp index 5f933f6..865b07a 100644 --- a/src/ecs/components/orbit-component.hpp +++ b/src/ecs/components/orbit-component.hpp @@ -20,15 +20,15 @@ #ifndef ANTKEEPER_ECS_ORBIT_COMPONENT_HPP #define ANTKEEPER_ECS_ORBIT_COMPONENT_HPP -#include "astro/orbit.hpp" +#include "physics/orbit/elements.hpp" namespace ecs { struct orbit_component { - astro::orbital_elements elements; - astro::orbital_elements rate; - astro::orbital_state state; + physics::orbit::elements elements; + physics::orbit::elements rate; + physics::orbit::elements state; }; } // namespace ecs diff --git a/src/ecs/systems/astronomy-system.cpp b/src/ecs/systems/astronomy-system.cpp index a6470e5..3dd6d23 100644 --- a/src/ecs/systems/astronomy-system.cpp +++ b/src/ecs/systems/astronomy-system.cpp @@ -18,12 +18,13 @@ */ #include "ecs/systems/astronomy-system.hpp" -#include "coordinates/coordinates.hpp" #include "astro/apparent-size.hpp" #include "ecs/components/celestial-body-component.hpp" #include "ecs/components/transform-component.hpp" #include "renderer/passes/sky-pass.hpp" #include "color/color.hpp" +#include "physics/orbit/orbit.hpp" +#include "geom/cartesian.hpp" #include namespace ecs { @@ -49,48 +50,127 @@ void astronomy_system::update(double t, double dt) // Add scaled timestep to current time set_universal_time(universal_time + dt * days_per_timestep); + set_universal_time(0.0); + // Update horizontal (topocentric) positions of intrasolar celestial bodies registry.view().each( [&](ecs::entity entity, auto& body, auto& transform) { + double time_correction = observer_location[2] / (math::two_pi / 24.0); + double local_jd = universal_time + time_correction / 24.0 - 0.5; + double local_time = (local_jd - std::floor(local_jd)) * 24.0; + double local_lst = local_time / 24.0f * math::two_pi; + // Transform orbital position from ecliptic space to horizontal space - double3 horizontal = ecliptic_to_horizontal * body.orbital_state.r; + //double3 horizontal = ecliptic_to_horizontal * body.orbital_state.r; + double3 horizontal = ecliptic_to_horizontal * double3{1, 0, 0}; // Subtract observer's radial distance (planet radius + observer's altitude) - horizontal.z -= observer_location[0]; + //horizontal.z -= observer_location[0]; - // Convert rectangular horizontal coordinates to spherical - double3 spherical = coordinates::rectangular::to_spherical(horizontal); + // Convert Cartesian horizontal coordinates to spherical + double3 spherical = geom::cartesian::to_spherical(horizontal); // Find angular radius double angular_radius = astro::find_angular_radius(body.radius, spherical[0]); // Transform into local coordinates - const double3x3 horizontal_to_local = - { - 0.0, 0.0, -1.0, - 1.0, 0.0, 0.0, - 0.0, 1.0, 0.0 - }; + const double3x3 horizontal_to_local = math::rotate_x(-math::half_pi) * math::rotate_z(-math::half_pi); double3 translation = horizontal_to_local * horizontal; double3x3 rotation = horizontal_to_local * ecliptic_to_horizontal; + // Set local transform of transform component transform.local.translation = math::type_cast(translation); - transform.local.rotation = math::type_cast(math::quaternion_cast(rotation)); + transform.local.rotation = math::normalize(math::type_cast(math::quaternion_cast(rotation))); transform.local.scale = math::type_cast(double3{body.radius, body.radius, body.radius}); if (sun_light != nullptr) { - math::quaternion sun_azimuth_rotation = math::angle_axis(static_cast(spherical.z), float3{0, 1, 0}); - math::quaternion sun_elevation_rotation = math::angle_axis(static_cast(spherical.y), float3{-1, 0, 0}); - math::quaternion sun_az_el_rotation = math::normalize(sun_azimuth_rotation * sun_elevation_rotation); + const double universal_time_cy = universal_time * 2.7397e-5; + const double3 solar_system_barycenter = {0, 0, 0}; + + physics::orbit::elements earth_elements; + earth_elements.a = 1.00000261 + 0.00000562 * universal_time_cy; + earth_elements.e = 0.01671123 + -0.00004392 * universal_time_cy; + + earth_elements.i = math::radians(-0.00001531) + math::radians(-0.01294668) * universal_time_cy; + earth_elements.raan = 0.0; + + const double earth_elements_mean_longitude = math::radians(100.46457166) + math::radians(35999.37244981) * universal_time_cy; + const double earth_elements_longitude_perihelion = math::radians(102.93768193) + math::radians(0.32327364) * universal_time_cy; + + earth_elements.w = earth_elements_longitude_perihelion - earth_elements.raan; + earth_elements.ta = earth_elements_mean_longitude - earth_elements_longitude_perihelion; + + + // Calculate semi-minor axis, b + double b = physics::orbit::derive_semiminor_axis(earth_elements.a, earth_elements.e); + + // Solve Kepler's equation for eccentric anomaly (E) + double ea = physics::orbit::kepler_ea(earth_elements.e, earth_elements.ta, 10, 1e-6); + + // Calculate radial distance, r; and true anomaly, v + double xv = earth_elements.a * (std::cos(ea) - earth_elements.e); + double yv = b * std::sin(ea); + double r = std::sqrt(xv * xv + yv * yv); + double ta = std::atan2(yv, xv); + + // Position of the body in perifocal space + const math::vector3 earth_position_pqw = math::quaternion::rotate_z(ta) * math::vector3{r, 0, 0}; - //sun_az_el_rotation = math::angle_axis((float)universal_time * math::two_pi, float3{1, 0, 0}); - // - //sun_light->look_at({0, 0, 0}, {0, -1, 0}, {0, 0, 1}); + const double earth_axial_tilt = math::radians(23.45); + const double earth_axial_rotation = math::two_pi * (0.7790572732640 + 1.00273781191135448 * universal_time); + const double earth_radius_au = 4.2635e-5; + + const double observer_altitude = earth_radius_au; + const double observer_latitude = math::radians(0.0); + const double observer_longitude = math::radians(0.0); + + const physics::frame earth_inertial_to_pqw = physics::orbit::inertial::to_perifocal(solar_system_barycenter, earth_elements.raan, earth_elements.i, earth_elements.w); + const math::vector3 earth_position_inertial = earth_inertial_to_pqw.inverse() * earth_position_pqw; + + const math::vector3 sun_position_intertial = math::vector3{0, 0, 0}; + + const physics::frame earth_inertial_to_bci = physics::orbit::inertial::to_bci(earth_position_inertial, earth_elements.i, earth_axial_tilt); + const physics::frame earth_inertial_to_bcbf = physics::orbit::inertial::to_bcbf(earth_position_inertial, earth_elements.i, earth_axial_tilt, earth_axial_rotation); + const physics::frame earth_bcbf_to_topo = physics::orbit::bcbf::to_topocentric(observer_altitude, observer_latitude, observer_longitude); + + const math::vector3 sun_position_earth_bci = earth_inertial_to_bci * sun_position_intertial; + const math::vector3 sun_position_earth_bcbf = earth_inertial_to_bcbf * sun_position_intertial; + const math::vector3 sun_position_earth_topo = earth_bcbf_to_topo * sun_position_earth_bcbf; + + const math::vector3 sun_radec = geom::cartesian::to_spherical(sun_position_earth_bci); + const math::vector3 sun_azel = geom::cartesian::to_spherical(sun_position_earth_topo); + + const double sun_az = sun_azel.z; + const double sun_el = sun_azel.y; + + double sun_ra = sun_radec.z; + const double sun_dec = sun_radec.y; + + if (sun_ra < 0.0) + sun_ra += math::two_pi; + + std::cout << "ra: " << (sun_ra / math::two_pi * 24.0) << "; dec: " << math::degrees(sun_dec) << std::endl; + std::cout << "az: " << math::degrees(math::pi - sun_az) << "; el: " << math::degrees(sun_el) << std::endl; + + + float az = spherical.z; + float el = spherical.y; + + if (az < 0.0f) + az += math::two_pi; + + //std::cout << "local: " << translation << std::endl; + //std::cout << "az: " << math::degrees(az) << "; "; + //std::cout << "el: " << math::degrees(el) << std::endl; + + math::quaternion sun_azimuth_rotation = math::angle_axis(static_cast(spherical.z), float3{0, 1, 0}); + math::quaternion sun_elevation_rotation = math::angle_axis(static_cast(spherical.y), float3{1, 0, 0}); + math::quaternion sun_az_el_rotation = math::normalize(sun_azimuth_rotation * sun_elevation_rotation); // Set sun color float cct = 3000.0f + std::sin(spherical.y) * 5000.0f; @@ -104,14 +184,16 @@ void astronomy_system::update(double t, double dt) sun_light->set_intensity(intensity); - sun_light->set_translation({0, 500, 0}); - //sun_light->set_rotation(math::look_rotation(math::normalize(transform.local.translation), {0, 1, 0})); - sun_light->set_rotation(sun_az_el_rotation); + //sun_light->set_translation({0, 500, 0}); + sun_light->set_translation(transform.local.translation); + //sun_light->set_rotation(transform.local.rotation); + //sun_light->set_rotation(sun_az_el_rotation); //sun_light->set_rotation(sun_elevation_rotation); + sun_light->set_rotation(math::look_rotation(math::normalize(-transform.local.translation), {0, 0, -1})); if (this->sky_pass) { - this->sky_pass->set_sun_coordinates(sun_az_el_rotation * float3{0, 0, -1}, {static_cast(spherical.z), static_cast(spherical.y)}); + this->sky_pass->set_sun_coordinates(transform.local.rotation * float3{0, 0, -1}, {static_cast(spherical.z), static_cast(spherical.y)}); } } }); @@ -186,7 +268,7 @@ void astronomy_system::update_sidereal_time() void astronomy_system::update_ecliptic_to_horizontal() { - ecliptic_to_horizontal = coordinates::rectangular::ecliptic::to_horizontal(obliquity, observer_location[1], lst); + //ecliptic_to_horizontal = coordinates::rectangular::ecliptic::to_horizontal(obliquity, observer_location[1], lst); } } // namespace ecs diff --git a/src/ecs/systems/solar-system.cpp b/src/ecs/systems/solar-system.cpp index 322ade7..59966ab 100644 --- a/src/ecs/systems/solar-system.cpp +++ b/src/ecs/systems/solar-system.cpp @@ -18,9 +18,6 @@ */ #include "ecs/systems/solar-system.hpp" -#include "coordinates/coordinates.hpp" -#include "astro/orbit.hpp" -#include "astro/constants.hpp" #include "ecs/components/celestial-body-component.hpp" #include "ecs/entity.hpp" @@ -45,16 +42,16 @@ void solar_system::update(double t, double dt) registry.view().each( [&](ecs::entity entity, auto& body) { - astro::orbital_elements elements = body.orbital_elements; + auto elements = body.orbital_elements; elements.a += body.orbital_rate.a * universal_time; - elements.ec += body.orbital_rate.ec * universal_time; + elements.e += body.orbital_rate.e * universal_time; elements.w += body.orbital_rate.w * universal_time; - elements.ma += body.orbital_rate.ma * universal_time; + elements.ta += body.orbital_rate.ta * universal_time; elements.i += body.orbital_rate.i * universal_time; - elements.om += body.orbital_rate.om * universal_time; + elements.raan += body.orbital_rate.raan * universal_time; // Calculate ecliptic orbital position - body.orbital_state.r = astro::orbital_elements_to_ecliptic(elements, ke_tolerance, ke_iterations); + //body.orbital_state.r = astro::orbital_elements_to_ecliptic(elements, ke_tolerance, ke_iterations); }); } diff --git a/src/game/states/play-state.cpp b/src/game/states/play-state.cpp index 321cb7b..15e9ce2 100644 --- a/src/game/states/play-state.cpp +++ b/src/game/states/play-state.cpp @@ -119,18 +119,18 @@ void play_state_enter(game_context* ctx) { ecs::celestial_body_component sun_body; sun_body.orbital_elements.a = 1.0; - sun_body.orbital_elements.ec = 0.016709; + sun_body.orbital_elements.e = 0.016709; sun_body.orbital_elements.w = math::radians(282.9404); - sun_body.orbital_elements.ma = math::radians(356.0470); + sun_body.orbital_elements.ta = math::radians(356.0470); sun_body.orbital_elements.i = 0.0; - sun_body.orbital_elements.om = 0.0; + sun_body.orbital_elements.raan = 0.0; sun_body.orbital_rate.a = 0.0; - sun_body.orbital_rate.ec = -1.151e-9; + sun_body.orbital_rate.e = -1.151e-9; sun_body.orbital_rate.w = math::radians(4.70935e-5); - sun_body.orbital_rate.ma = math::radians(0.9856002585); + sun_body.orbital_rate.ta = math::radians(0.9856002585); sun_body.orbital_rate.i = 0.0; - sun_body.orbital_rate.om = 0.0; + sun_body.orbital_rate.raan = 0.0; ecs::transform_component sun_transform; sun_transform.local = math::identity_transform; @@ -161,7 +161,7 @@ void play_state_enter(game_context* ctx) ctx->solar_system->set_universal_time(0.0); - ctx->astronomy_system->set_observer_location(double3{4.26352e-5, ctx->biome->location[0], ctx->biome->location[1]}); + ctx->astronomy_system->set_observer_location(double3{4.26352e-5, math::radians(10.0f), 0.0f}); ctx->astronomy_system->set_universal_time(0.0); ctx->astronomy_system->set_obliquity(math::radians(23.4393)); ctx->astronomy_system->set_axial_rotation_at_epoch(math::radians(280.4606)); @@ -228,7 +228,7 @@ void play_state_enter(game_context* ctx) // Create ant-hill auto ant_hill_entity = ant_hill_archetype->create(ecs_registry); - ecs::command::place(ecs_registry, ant_hill_entity, {0, 0}); + ecs::command::place(ecs_registry, ant_hill_entity, {0, -40}); // Creat nest auto nest_entity = nest_archetype->create(ecs_registry); @@ -346,7 +346,7 @@ void play_state_enter(game_context* ctx) { auto larva = larva_archetype->create(ecs_registry); ecs::command::assign_render_layers(ecs_registry, larva, 1); - ecs::command::warp_to(ecs_registry, larva, {50, 0.1935f, 10}); + ecs::command::warp_to(ecs_registry, larva, {50, 0.1935f, 0}); //auto& transform = ecs_registry.get(larva_entity); //transform.transform = math::identity_transform; //transform.transform.translation = nest->get_shaft_position(*central_shaft, central_shaft->depth[1]); diff --git a/src/geom/cartesian.hpp b/src/geom/cartesian.hpp new file mode 100644 index 0000000..360990e --- /dev/null +++ b/src/geom/cartesian.hpp @@ -0,0 +1,58 @@ +/* + * 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_GEOM_CARTESIAN_HPP +#define ANTKEEPER_GEOM_CARTESIAN_HPP + +#include "utility/fundamental-types.hpp" +#include + +namespace geom { + +/// Functions which operate on Cartesian (rectangular) coordinates. +namespace cartesian { + +/** + * Converts Cartesian (rectangular) coordinates to spherical coordinates. + * + * @param v Cartesian coordinates. + * @return Spherical coordinates, in the ISO order of radial distance, polar angle (radians), and azimuthal angle (radians). + * + * @see geom::coordinates::spherical + */ +template +math::vector3 to_spherical(const math::vector3& v); + +template +math::vector3 to_spherical(const math::vector3& v) +{ + const T xx_yy = v.x * v.x + v.y * v.y; + + return math::vector3 + { + std::sqrt(xx_yy + v.z * v.z), + std::atan2(v.z, std::sqrt(xx_yy)), + std::atan2(v.y, v.x) + }; +} + +} // namespace cartesian +} // namespace geom + +#endif // ANTKEEPER_GEOM_CARTESIAN_HPP diff --git a/src/geom/geom.hpp b/src/geom/geom.hpp index 9e2ab23..8c7c731 100644 --- a/src/geom/geom.hpp +++ b/src/geom/geom.hpp @@ -26,6 +26,7 @@ namespace geom {} #include "aabb.hpp" #include "bounding-volume.hpp" #include "convex-hull.hpp" +#include "cartesian.hpp" #include "csg.hpp" #include "intersection.hpp" #include "marching-cubes.hpp" @@ -39,6 +40,7 @@ namespace geom {} #include "ray.hpp" #include "sdf.hpp" #include "sphere.hpp" +#include "spherical.hpp" #include "view-frustum.hpp" #endif // ANTKEEPER_GEOM_HPP diff --git a/src/coordinates/spherical.hpp b/src/geom/spherical.hpp similarity index 72% rename from src/coordinates/spherical.hpp rename to src/geom/spherical.hpp index 649439c..a2dda81 100644 --- a/src/coordinates/spherical.hpp +++ b/src/geom/spherical.hpp @@ -17,30 +17,30 @@ * along with Antkeeper source code. If not, see . */ -#ifndef ANTKEEPER_COORDINATES_SPHERICAL_HPP -#define ANTKEEPER_COORDINATES_SPHERICAL_HPP +#ifndef ANTKEEPER_GEOM_SPHERICAL_HPP +#define ANTKEEPER_GEOM_SPHERICAL_HPP #include "utility/fundamental-types.hpp" #include -namespace coordinates { +namespace geom { -/// Spherical coordinate systems. +/// Functions which operate on spherical coordinates. namespace spherical { /** - * Converts spherical coordinates to rectangular coordinates. + * Converts spherical coordinates to Cartesian (rectangular) coordinates. * * @param v Spherical coordinates, in the ISO order of radial distance, polar angle (radians), and azimuthal angle (radians). - * @return Rectangular coordinates. + * @return Cartesian coordinates. * - * @see coordinates::rectangular + * @see geom::coordinates::cartesian */ template -math::vector3 to_rectangular(const math::vector3& v); +math::vector3 to_cartesian(const math::vector3& v); template -math::vector3 to_rectangular(const math::vector3& v) +math::vector3 to_cartesian(const math::vector3& v) { const T x = v[0] * std::cos(v[1]); @@ -53,6 +53,6 @@ math::vector3 to_rectangular(const math::vector3& v) } } // namespace spherical -} // namespace coordinates +} // namespace geom -#endif // ANTKEEPER_COORDINATES_SPHERICAL_HPP +#endif // ANTKEEPER_GEOM_SPHERICAL_HPP diff --git a/src/math/matrix-functions.hpp b/src/math/matrix-functions.hpp index e76822e..acdc5f4 100644 --- a/src/math/matrix-functions.hpp +++ b/src/math/matrix-functions.hpp @@ -251,6 +251,33 @@ matrix resize(const matrix& m); template matrix rotate(const matrix& m, T angle, const vector& axis); +/** + * Produces a matrix which rotates Cartesian coordinates about the x-axis by a given angle. + * + * @param angle Angle of rotation, in radians. + * @return Rotation matrix. + */ +template +matrix3 rotate_x(T angle); + +/** + * Produces a matrix which rotates Cartesian coordinates about the y-axis by a given angle. + * + * @param angle Angle of rotation, in radians. + * @return Rotation matrix. + */ +template +matrix3 rotate_y(T angle); + +/** + * Produces a matrix which rotates Cartesian coordinates about the z-axis by a given angle. + * + * @param angle Angle of rotation, in radians. + * @return Rotation matrix. + */ +template +matrix3 rotate_z(T angle); + /** * Scales a matrix. * @@ -734,6 +761,48 @@ matrix rotate(const matrix& m, T angle, const vector& ax return mul(m, rotation); } +template +matrix3 rotate_x(T angle) +{ + const T c = std::cos(angle); + const T s = std::sin(angle); + + return matrix3 + { + T(1), T(0), T(0), + T(0), c, s, + T(0), -s, c + }; +} + +template +matrix3 rotate_y(T angle) +{ + const T c = std::cos(angle); + const T s = std::sin(angle); + + return matrix3 + { + c, T(0), -s, + T(0), T(1), T(0), + s, T(0), c + }; +} + +template +matrix3 rotate_z(T angle) +{ + const T c = std::cos(angle); + const T s = std::sin(angle); + + return matrix3 + { + c, s, T(0), + -s, c, T(0), + T(0), T(0), T(1) + }; +} + template matrix scale(const matrix& m, const vector& v) { diff --git a/src/math/quaternion-type.hpp b/src/math/quaternion-type.hpp index 24a7e2a..f895155 100644 --- a/src/math/quaternion-type.hpp +++ b/src/math/quaternion-type.hpp @@ -20,6 +20,8 @@ #ifndef ANTKEEPER_MATH_QUATERNION_TYPE_HPP #define ANTKEEPER_MATH_QUATERNION_TYPE_HPP +#include + namespace math { /** @@ -35,9 +37,59 @@ struct quaternion scalar_type x; scalar_type y; scalar_type z; + + /// Returns the rotation identity quaternion. + static constexpr quaternion identity(); + + /** + * Returns a quaternion representing a rotation about the x-axis. + * + * @param angle Angle of rotation, in radians. + * @return Quaternion representing an x-axis rotation. + */ + static quaternion rotate_x(scalar_type angle); + + /** + * Returns a quaternion representing a rotation about the y-axis. + * + * @param angle Angle of rotation, in radians. + * @return Quaternion representing an y-axis rotation. + */ + static quaternion rotate_y(scalar_type angle); + + /** + * Returns a quaternion representing a rotation about the z-axis. + * + * @param angle Angle of rotation, in radians. + * @return Quaternion representing an z-axis rotation. + */ + static quaternion rotate_z(scalar_type angle); +}; + +template +constexpr quaternion quaternion::identity() +{ + return {T(1), T(0), T(0), T(0)}; }; +template +quaternion quaternion::rotate_x(scalar_type angle) +{ + return {std::cos(angle * T(0.5)), std::sin(angle * T(0.5)), T(0), T(0)}; +} + +template +quaternion quaternion::rotate_y(scalar_type angle) +{ + return {std::cos(angle * T(0.5)), T(0), std::sin(angle * T(0.5)), T(0)}; +} + +template +quaternion quaternion::rotate_z(scalar_type angle) +{ + return {std::cos(angle * T(0.5)), T(0), T(0), std::sin(angle * T(0.5))}; +} + } // namespace math #endif // ANTKEEPER_MATH_QUATERNION_TYPE_HPP - diff --git a/src/physics/frame.hpp b/src/physics/frame.hpp new file mode 100644 index 0000000..a368152 --- /dev/null +++ b/src/physics/frame.hpp @@ -0,0 +1,134 @@ +/* + * 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_FRAME_HPP +#define ANTKEEPER_PHYSICS_FRAME_HPP + +#include "utility/fundamental-types.hpp" +#include + +namespace physics { + +/** + * 3-dimensional frame of reference. + * + * @tparam T Scalar type. + */ +template +struct frame +{ +public: + /// Scalar type. + typedef T scalar_type; + + /// Vector type. + typedef math::vector3 vector_type; + + /// Quaternion type. + typedef math::quaternion quaternion_type; + + /// Transformation matrix type. + typedef math::matrix matrix_type; + + /// Vector representing a translation from the origin of this frame to the origin of the parent frame. + vector_type translation; + + /// Quaternion representing a rotation from the orientation of this frame to the orientation of the parent frame. + quaternion_type rotation; + + /// Returns the inverse of this frame. + frame inverse() const; + + /// Returns a transformation matrix representation of this frame. + matrix_type matrix() const; + + /** + * Transforms a vector from the parent frame space to this frame's space. + * + * @param v Vector in parent frame space. + * @return Vector in this frame's space. + */ + vector_type transform(const vector_type& v) const; + + /** + * Transforms a frame from the parent frame space to this frame's space. + * + * @param f Frame in parent frame space. + * @return Frame in this frame's space. + */ + frame transform(const frame& f) const; + + /// @copydoc frame::transform(const vector_type&) const + vector_type operator*(const vector_type& v) const; + + /// @copydoc frame::transform(const frame&) const + frame operator*(const frame& f) const; +}; + +template +frame frame::inverse() const +{ + const quaternion_type inverse_rotation = math::conjugate(rotation); + const vector_type inverse_translation = -(inverse_rotation * translation); + return frame{inverse_translation, inverse_rotation}; +} + +template +typename frame::matrix_type frame::matrix() const +{ + matrix_type m = math::resize<4, 4>(math::matrix_cast(rotation)); + + m[3].x = translation.x; + m[3].y = translation.y; + m[3].z = translation.z; + + return m; +} + +template +typename frame::vector_type frame::transform(const vector_type& v) const +{ + return translation + rotation * v; +} + +template +frame frame::transform(const frame& f) const +{ + return frame + { + transform(f.translation), + math::normalize(rotation * f.rotation) + }; +} + +template +typename frame::vector_type frame::operator*(const vector_type& v) const +{ + return transform(v); +} + +template +frame frame::operator*(const frame& f) const +{ + return transform(f); +} + +} // namespace physics + +#endif // ANTKEEPER_PHYSICS_FRAME_HPP diff --git a/src/physics/orbit/elements.hpp b/src/physics/orbit/elements.hpp new file mode 100644 index 0000000..2e26b0e --- /dev/null +++ b/src/physics/orbit/elements.hpp @@ -0,0 +1,95 @@ +/* + * 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_ELEMENTS_HPP +#define ANTKEEPER_PHYSICS_ORBIT_ELEMENTS_HPP + +#include "utility/fundamental-types.hpp" +#include "physics/orbit/state.hpp" +#include + +namespace physics { +namespace orbit { + +/** + * Set of six Keplerian elements required to uniquely identify an orbit. + * + * @tparam T Scalar type. + */ +template +struct elements +{ + /// Scalar type. + typedef T scalar_type; + + /// Eccentricity (e). + scalar_type e; + + /// Semimajor axis (a). + scalar_type a; + + /// Inclination (i), in radians. + scalar_type i; + + /// Right ascension of the ascending node (OMEGA), in radians. + scalar_type raan; + + /// Argument of periapsis (omega), in radians. + scalar_type w; + + /// True anomaly (nu) at epoch, in radians. + scalar_type ta; +}; + +/** + * Derives the longitude of the periapsis (pomega) of an orbit, given the argument of periapsis (omega) and longitude of the ascending node (OMEGA). + * + * @param w Argument of periapsis (omega), in radians. + * @param raan Right ascension of the ascending node (OMEGA), in radians. + * @return Longitude of the periapsis (pomega), in radians. + */ +template +T derive_longitude_periapsis(T w, T raan); + +/** + * Derives the semiminor axis (b) of an orbit, given the semimajor axis (a) and eccentricity (e). + * + * @param a Semimajor axis (a). + * @param e Eccentricity (e). + * @return Semiminor axis (b). + */ +template +T derive_semiminor_axis(T a, T e); + +template +T derive_longitude_periapsis(T w, T raan) +{ + return w + raan; +} + +template +T derive_semiminor_axis(T a, T e) +{ + return a * std::sqrt(T(1) - e * e); +} + +} // namespace orbit +} // namespace physics + +#endif // ANTKEEPER_PHYSICS_ORBIT_ELEMENTS_HPP diff --git a/src/physics/orbit/frames.hpp b/src/physics/orbit/frames.hpp new file mode 100644 index 0000000..74f6f57 --- /dev/null +++ b/src/physics/orbit/frames.hpp @@ -0,0 +1,131 @@ +/* + * 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_FRAMES_HPP +#define ANTKEEPER_PHYSICS_ORBIT_FRAMES_HPP + +#include "physics/frame.hpp" + +namespace physics { +namespace orbit { + +/// Inertial right-handed coordinate system +namespace inertial { + + /** + * Constucts a reference frame which transforms coordinates from inertial space into perifocal space. + * + * @param focus Cartesian coordinates of the focus of the orbit, in the parent space. + * @param raan Right ascension of the ascending node (OMEGA), in radians. + * @param i Orbital inclination (i), in radians. + * @param w Argument of periapsis (omega), in radians. + * @return Perifocal frame. + */ + template + physics::frame to_perifocal(const math::vector3& focus, T raan, T i, T w) + { + const math::quaternion rotation = math::normalize + ( + math::quaternion::rotate_z(raan) * + math::quaternion::rotate_x(i) * + math::quaternion::rotate_z(w) + ); + + return physics::frame{focus, rotation}.inverse(); + } + + /** + * Constructs a reference frame which transforms coordinates from inertial space to body-centered inertial space. + * + * @param r Cartesian position vector (r) of the center of the body, in inertial space. + * @param i Orbital inclination (i), in radians. + * @param axial_tilt Angle between the body's rotational axis and its orbital axis, in radians. + * @return Body-centered inertial frame. + */ + template + physics::frame to_bci(const math::vector3& r, T i, T axial_tilt) + { + return physics::frame{r, math::quaternion::rotate_x(-axial_tilt - i)}.inverse(); + } + + /** + * Constructs a reference frame which transforms coordinates from inertial space to body-centered, body-fixed space. + * + * @param r Cartesian position vector (r) of the center of the body, in inertial space. + * @param i Orbital inclination (i), in radians. + * @param axial_tilt Angle between the orbiting body's rotational axis and its orbital axis, in radians. + * @param axial_rotation Angle of rotation about the orbiting body's rotational axis, in radians. + */ + template + frame to_bcbf(const math::vector3& r, T i, T axial_tilt, T axial_rotation) + { + const math::quaternion rotation = math::normalize + ( + math::quaternion::rotate_x(-axial_tilt - i) * + math::quaternion::rotate_z(axial_rotation) + ); + + return physics::frame{r, rotation}.inverse(); + } + +} // namespace inertial + +/// Perifocal right-handed coordinate system in which the x-axis points toward the periapsis of the orbit, the y-axis has a true anomaly of 90 degrees past the periapsis, and the z-axis is the angular momentum vector, which is orthogonal to the orbital plane. +namespace perifocal { + + + +} // namespace perifocal + +/// Non-inertial body-centered, body-fixed right-handed coordinate system. The x-axis is orthogonal to the intersection of the prime meridian and the equator. The z-axis points toward the positive pole. The y-axis is right-hand orthogonal to the xz-plane. +namespace bcbf { + + /** + * Constructs a reference frame which transforms coordinates from BCBF space to topocentric space. + * + * @param distance Radial distance of the observer from the center of the body. + * @param latitude Latitude of the observer, in radians. + * @param longitude Longitude of the obserer, in radians. + * @return Topocentric frame. + */ + template + physics::frame to_topocentric(T distance, T latitude, T longitude) + { + const math::vector3 translation = {T(0), T(0), distance}; + + const math::quaternion rotation = math::normalize + ( + math::quaternion::rotate_z(longitude) * + math::quaternion::rotate_y(math::half_pi - latitude) + ); + + return physics::frame{rotation * translation, rotation}.inverse(); + } + +} // namespace bcbf + +/// Non-inertial topocentric right-handed coordinate system. Topocentric frames are constructed as south-east-zenith (SEZ) frames; the x-axis points south, the y-axis points east, and the z-axis points toward the zenith (orthogonal to reference spheroid). +namespace topocentric { + +} // namespace topocentric + +} // namespace orbit +} // namespace physics + +#endif // ANTKEEPER_PHYSICS_ORBIT_FRAMES_HPP diff --git a/src/physics/orbit/kepler.hpp b/src/physics/orbit/kepler.hpp new file mode 100644 index 0000000..637648f --- /dev/null +++ b/src/physics/orbit/kepler.hpp @@ -0,0 +1,79 @@ +/* + * 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_KEPLER_HPP +#define ANTKEEPER_PHYSICS_ORBIT_KEPLER_HPP + +#include + +namespace physics { +namespace orbit { + +/** + * Iteratively solves Kepler's equation for eccentric anomaly (E). + * + * @param ec Eccentricity (e). + * @param ma Mean anomaly (M). + * @param iterations Maximum number of iterations. + * @param tolerance Solution error tolerance. + * @return Eccentric anomaly (E). + */ +template +T kepler_ea(T ec, T ma, std::size_t iterations, T tolerance = T(0)); + +/** + * Solves Kepler's equation for mean anomaly (M). + * + * @param ec Eccentricity (e). + * @param ea Eccentric anomaly (E). + * @return Mean anomaly (M). + */ +template +T kepler_ma(T ec, T ea); + +template +T kepler_ea(T ec, T ma, std::size_t iterations, T tolerance) +{ + // Approximate eccentric anomaly, E + T ea0 = ma + ec * std::sin(ma) * (T(1.0) + ec * std::cos(ma)); + + // Iteratively converge E0 and E1 + for (std::size_t i = 0; i < iterations; ++i) + { + const T ea1 = ea0 - (ea0 - ec * std::sin(ea0) - ma) / (T(1.0) - ec * std::cos(ea0)); + const T error = std::abs(ea1 - ea0); + ea0 = ea1; + + if (error < tolerance) + break; + } + + return ea0; +} + +template +T kepler_ma(T ec, T ea) +{ + return ea - ec * std::sin(ea); +} + +} // namespace orbit +} // namespace physics + +#endif // ANTKEEPER_PHYSICS_ORBIT_KEPLER_HPP diff --git a/src/physics/orbit/orbit.hpp b/src/physics/orbit/orbit.hpp new file mode 100644 index 0000000..3871e40 --- /dev/null +++ b/src/physics/orbit/orbit.hpp @@ -0,0 +1,39 @@ +/* + * 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_HPP +#define ANTKEEPER_PHYSICS_ORBIT_HPP + +namespace physics { + +/** + * Orbital mechanics. + * + * @see Curtis, H. D. (2005). *Orbital mechanics for engineering students*. Amsterdam: Elsevier Butterworth Heinemann. + */ +namespace orbit {} + +} // namespace physics + +#include "physics/orbit/elements.hpp" +#include "physics/orbit/frames.hpp" +#include "physics/orbit/kepler.hpp" +#include "physics/orbit/state.hpp" + +#endif // ANTKEEPER_PHYSICS_ORBIT_HPP diff --git a/src/physics/orbit/state.hpp b/src/physics/orbit/state.hpp new file mode 100644 index 0000000..2e514c1 --- /dev/null +++ b/src/physics/orbit/state.hpp @@ -0,0 +1,52 @@ +/* + * 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_STATE_HPP +#define ANTKEEPER_PHYSICS_ORBIT_STATE_HPP + +#include "utility/fundamental-types.hpp" + +namespace physics { +namespace orbit { + +/** + * Pair of orbital state Cartesian position (r) and velocity (v) vectors. + * + * @tparam T Scalar type. + */ +template +struct state +{ + /// Scalar type. + typedef T scalar_type; + + /// Vector type. + typedef math::vector3 vector_type; + + /// Cartesian orbital position vector (r). + vector_type r; + + /// Cartesian orbital velocity vector (v). + vector_type v; +}; + +} // namespace orbit +} // namespace physics + +#endif // ANTKEEPER_PHYSICS_ORBIT_STATE_HPP diff --git a/src/coordinates/coordinates.hpp b/src/physics/physics.hpp similarity index 68% rename from src/coordinates/coordinates.hpp rename to src/physics/physics.hpp index 4751438..4f04a13 100644 --- a/src/coordinates/coordinates.hpp +++ b/src/physics/physics.hpp @@ -17,18 +17,13 @@ * along with Antkeeper source code. If not, see . */ -#ifndef ANTKEEPER_COORDINATES_HPP -#define ANTKEEPER_COORDINATES_HPP +#ifndef ANTKEEPER_PHYSICS_HPP +#define ANTKEEPER_PHYSICS_HPP -#include "utility/fundamental-types.hpp" +/// Physics +namespace physics {} -/// Functions for converting between coordinate systems. -namespace coordinates {} +#include "frame.hpp" +#include "orbit/orbit.hpp" -#include "ecliptic.hpp" -#include "equatorial.hpp" -#include "horizontal.hpp" -#include "rectangular.hpp" -#include "spherical.hpp" - -#endif // ANTKEEPER_COORDINATES_HPP +#endif // ANTKEEPER_PHYSICS_HPP diff --git a/src/renderer/passes/sky-pass.cpp b/src/renderer/passes/sky-pass.cpp index c15c033..652bba2 100644 --- a/src/renderer/passes/sky-pass.cpp +++ b/src/renderer/passes/sky-pass.cpp @@ -40,7 +40,8 @@ #include "color/color.hpp" #include "astro/illuminance.hpp" #include "math/interpolation.hpp" -#include "coordinates/coordinates.hpp" +#include "geom/cartesian.hpp" +#include "geom/spherical.hpp" #include #include #include @@ -103,13 +104,13 @@ sky_pass::sky_pass(gl::rasterizer* rasterizer, const gl::framebuffer* framebuffe } catch (const std::exception& e) {} - + // Convert right ascension and declination from degrees to radians ra = math::wrap_radians(math::radians(ra)); dec = math::wrap_radians(math::radians(dec)); // Transform spherical equatorial coordinates to rectangular equatorial coordinates - double3 position = coordinates::spherical::to_rectangular(double3{1.0, dec, ra}); + double3 position = geom::spherical::to_cartesian(double3{1.0, dec, ra}); // Convert color index to color temperature double cct = color::index::bv_to_cct(bv_color); @@ -289,24 +290,19 @@ void sky_pass::render(render_context* context) const { float star_distance = (clip_near + clip_far) * 0.5f; - double lat = math::radians(30.0); + double lat = math::radians(1.0); double lst = time_of_day / 24.0f * math::two_pi; //std::cout << "lst: " << lst << std::endl; + /* double3x3 equatorial_to_horizontal = coordinates::rectangular::equatorial::to_horizontal(lat, lst); - const double3x3 horizontal_to_local = - { - 0.0, 0.0, -1.0, - 1.0, 0.0, 0.0, - 0.0, 1.0, 0.0 - }; + const double3x3 horizontal_to_local = coordinates::rectangular::rotate_x(-math::half_pi) * coordinates::rectangular::rotate_z(-math::half_pi); double3x3 rotation = horizontal_to_local * equatorial_to_horizontal; - //math::quaternion rotation_y = math::angle_axis(time_of_day / 24.0f * -math::two_pi, {0, 1, 0}); - //math::quaternion rotation_x = math::angle_axis(-math::half_pi + math::radians(30.0f), {1, 0, 0}); model = math::type_cast(math::scale(math::resize<4, 4>(rotation), double3{star_distance, star_distance, star_distance}));; + */ //math::transform star_transform; //star_transform.translation = {0.0, 0.0, 0.0};