From 1ef5ff08b02accbe231fcdb6859b6fd7a3aa3a41 Mon Sep 17 00:00:00 2001 From: "C. J. Howard" Date: Wed, 7 Oct 2020 22:08:51 -0700 Subject: [PATCH] Add celestial time functions --- CMakeLists.txt | 1 - src/game/astronomy/celestial-coordinates.cpp | 24 +- src/game/astronomy/celestial-coordinates.hpp | 24 +- src/game/astronomy/celestial-mechanics.cpp | 104 +++++++++ src/game/astronomy/celestial-mechanics.hpp | 32 +++ src/game/astronomy/celestial-time.cpp | 62 ++++++ src/game/astronomy/celestial-time.hpp | 66 ++++++ src/game/states/play-state.cpp | 18 +- src/game/systems/weather-system.cpp | 222 ++----------------- src/game/systems/weather-system.hpp | 2 +- 10 files changed, 318 insertions(+), 237 deletions(-) create mode 100644 src/game/astronomy/celestial-time.cpp create mode 100644 src/game/astronomy/celestial-time.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 18ea3f6..b80be9e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,7 +4,6 @@ option(VERSION_STRING "Project version string" "0.0.0") project(antkeeper VERSION ${VERSION_STRING} LANGUAGES CXX) - # Find dependency packages find_package(dr_wav REQUIRED CONFIG) find_package(stb REQUIRED CONFIG) diff --git a/src/game/astronomy/celestial-coordinates.cpp b/src/game/astronomy/celestial-coordinates.cpp index 194b798..7a9aa43 100644 --- a/src/game/astronomy/celestial-coordinates.cpp +++ b/src/game/astronomy/celestial-coordinates.cpp @@ -58,14 +58,14 @@ double3x3 ecliptic_to_equatorial(double ecl) }; } -double3x3 ecliptic_to_horizontal(double ecl, double lat, double lst) +double3x3 ecliptic_to_horizontal(double ecl, double lat, double lmst) { const double c_ecl = std::cos(ecl); const double s_ecl = std::sin(ecl); const double c_lat = std::cos(lat); const double s_lat = std::sin(lat); - const double c_lst = std::cos(lst); - const double s_lst = std::sin(lst); + const double c_lst = std::cos(lmst); + const double s_lst = std::sin(lmst); return double3x3 { @@ -88,12 +88,12 @@ double3x3 equatorial_to_ecliptic(double ecl) }; } -double3x3 equatorial_to_horizontal(double lat, double lst) +double3x3 equatorial_to_horizontal(double lat, double lmst) { const double c_lat = std::cos(lat); const double s_lat = std::sin(lat); - const double c_lst = std::cos(lst); - const double s_lst = std::sin(lst); + const double c_lst = std::cos(lmst); + const double s_lst = std::sin(lmst); return double3x3 { @@ -103,12 +103,12 @@ double3x3 equatorial_to_horizontal(double lat, double lst) }; } -double3x3 horizontal_to_equatorial(double lat, double lst) +double3x3 horizontal_to_equatorial(double lat, double lmst) { const double c_lat = std::cos(lat); const double s_lat = std::sin(lat); - const double c_lst = std::cos(lst); - const double s_lst = std::sin(lst); + const double c_lst = std::cos(lmst); + const double s_lst = std::sin(lmst); return double3x3 { @@ -118,14 +118,14 @@ double3x3 horizontal_to_equatorial(double lat, double lst) }; } -double3x3 horizontal_to_ecliptic(double ecl, double lat, double lst) +double3x3 horizontal_to_ecliptic(double ecl, double lat, double lmst) { const double c_ecl = std::cos(ecl); const double s_ecl = std::sin(ecl); const double c_lat = std::cos(lat); const double s_lat = std::sin(lat); - const double c_lst = std::cos(lst); - const double s_lst = std::sin(lst); + const double c_lst = std::cos(lmst); + const double s_lst = std::sin(lmst); return double3x3 { diff --git a/src/game/astronomy/celestial-coordinates.hpp b/src/game/astronomy/celestial-coordinates.hpp index 602cfc9..7b401a3 100644 --- a/src/game/astronomy/celestial-coordinates.hpp +++ b/src/game/astronomy/celestial-coordinates.hpp @@ -54,10 +54,10 @@ double3x3 ecliptic_to_equatorial(double ecl); * * @param ecl Obliquity of the ecliptic, in radians. * @param lat Observer's latitude, in radians. - * @param lst Local sidereal time. + * @param lmst Local mean sidereal time. * @return Transformation matrix. */ -double3x3 ecliptic_to_horizontal(double ecl, double lat, double lst); +double3x3 ecliptic_to_horizontal(double ecl, double lat, double lmst); /** * Produces a matrix which transforms rectangular coordinates from equatorial space to ecliptic space. @@ -71,29 +71,37 @@ double3x3 equatorial_to_ecliptic(double ecl); * Produces a matrix which transforms rectangular coordinates from equatorial space to horizontal space. * * @param lat Observer's latitude, in radians. - * @param lst Local sidereal time. + * @param lmst Local mean sidereal time. * @return Transformation matrix. */ -double3x3 equatorial_to_horizontal(double lat, double lst); +double3x3 equatorial_to_horizontal(double lat, double lmst); /** * Produces a matrix which transforms rectangular coordinates from horizontal space to equatorial space. * * @param lat Observer's latitude, in radians. - * @param lst Local sidereal time. + * @param lmst Local mean sidereal time. * @return Transformation matrix. */ -double3x3 horizontal_to_equatorial(double lat, double lst); +double3x3 horizontal_to_equatorial(double lat, double lmst); /** * Produces a matrix which transforms rectangular coordinates from horizontal space to ecliptic space. * * @param ecl Obliquity of the ecliptic, in radians. * @param lat Observer's latitude, in radians. - * @param lst Local sidereal time. + * @param lmst Local mean sidereal time. * @return Transformation matrix. */ -double3x3 horizontal_to_ecliptic(double ecl, double lat, double lst); +double3x3 horizontal_to_ecliptic(double ecl, double lat, double lmst); + +/// Matrix which transforms coordinates from horizontal space to a right-handed coordinate system. +constexpr double3x3 horizontal_to_right_handed = +{ + 0.0, 0.0, 1.0, + 1.0, 0.0, 0.0, + 0.0, -1.0, 0.0 +}; } // namespace ast diff --git a/src/game/astronomy/celestial-mechanics.cpp b/src/game/astronomy/celestial-mechanics.cpp index 5a12a81..ac29672 100644 --- a/src/game/astronomy/celestial-mechanics.cpp +++ b/src/game/astronomy/celestial-mechanics.cpp @@ -18,11 +18,115 @@ */ #include "celestial-mechanics.hpp" +#include "math/angles.hpp" #include namespace ast { +double approx_ecliptic_obliquity(double jd) +{ + return math::radians(23.4393 - 3.563e-7 * (jd - 2451545.0)); +} + +double3 approx_sun_ecliptic(double jd) +{ + const double t = (jd - 2451545.0) / 36525.0; + const double m = 6.24 + 628.302 * t; + + const double longitude = 4.895048 + 628.331951 * t + (0.033417 - 0.000084 * t) * std::sin(m) + 0.000351 * std::sin(m * 2.0); + const double latitude = 0.0; + const double distance = 1.000140 - (0.016708 - 0.000042 * t) * std::cos(m) - 0.000141 * std::cos(m * 2.0); + + double3 ecliptic; + ecliptic.x = distance * std::cos(longitude) * std::cos(latitude); + ecliptic.y = distance * std::sin(longitude) * std::cos(latitude); + ecliptic.z = distance * std::sin(latitude); + + return ecliptic; +} +double3 approx_moon_ecliptic(double jd) +{ + const double t = (jd - 2451545.0) / 36525.0; + const double l1 = 3.8104 + 8399.7091 * t; + const double m1 = 2.3554 + 8328.6911 * t; + const double m = 6.2300 + 628.3019 * t; + const double d = 5.1985 + 7771.3772 * t; + const double d2 = d * 2.0; + const double f = 1.6280 + 8433.4663 * t; + + const double longitude = l1 + + 0.1098 * std::sin(m1) + + 0.0222 * std::sin(d2 - m1) + + 0.0115 * std::sin(d2) + + 0.0037 * std::sin(m1 * 2.0) + - 0.0032 * std::sin(m) + - 0.0020 * std::sin(d2) + + 0.0010 * std::sin(d2 - m1 * 2.0) + + 0.0010 * std::sin(d2 - m - m1) + + 0.0009 * std::sin(d2 + m1) + + 0.0008 * std::sin(d2 - m) + + 0.0007 * std::sin(m1 - m) + - 0.0006 * std::sin(d) + - 0.0005 * std::sin(m + m1); + + const double latitude = 0.0895 * sin(f) + + 0.0049 * std::sin(m1 + f) + + 0.0048 * std::sin(m1 - f) + + 0.0030 * std::sin(d2 - f) + + 0.0010 * std::sin(d2 + f - m1) + + 0.0008 * std::sin(d2 - f - m1) + + 0.0006 * std::sin(d2 + f); + + const double r = 1.0 / (0.016593 + + 0.000904 * std::cos(m1) + + 0.000166 * std::cos(d2 - m1) + + 0.000137 * std::cos(d2) + + 0.000049 * std::cos(m1 * 2.0) + + 0.000015 * std::cos(d2 + m1) + + 0.000009 * std::cos(d2 - m)); + + double3 ecliptic; + ecliptic.x = r * std::cos(longitude) * std::cos(latitude); + ecliptic.y = r * std::sin(longitude) * std::cos(latitude); + ecliptic.z = r * std::sin(latitude); + + return ecliptic; +} + +double3x3 approx_moon_ecliptic_rotation(double jd) +{ + const double t = (jd - 2451545.0) / 36525.0; + const double l1 = 3.8104 + 8399.7091 * t; + const double f = 1.6280 + 8433.4663 * t; + + const double az0 = f + math::pi; + const double ax = 0.026920; + const double az1 = l1 - f; + + double3x3 rz0 = + { + cos(az0), -sin(az0), 0, + sin(az0), cos(az0), 0, + 0, 0, 1 + }; + + double3x3 rx = + { + 1, 0, 0, + 0, cos(ax), -sin(ax), + 0, sin(ax), cos(ax) + }; + + double3x3 rz1 = + { + cos(az1), -sin(az1), 0, + sin(az1), cos(az1), 0, + 0, 0, 1 + }; + + return rz0 * rx * rz1; +} } // namespace ast diff --git a/src/game/astronomy/celestial-mechanics.hpp b/src/game/astronomy/celestial-mechanics.hpp index b4ae188..a2f938f 100644 --- a/src/game/astronomy/celestial-mechanics.hpp +++ b/src/game/astronomy/celestial-mechanics.hpp @@ -25,6 +25,38 @@ namespace ast { +/** + * Approximates the obliquity of the ecliptic. + * + * @param jd Julian date. + * @return Obliquity of the ecliptic, in radians. + */ +double approx_ecliptic_obliquity(double jd); + +/** + * Approximates the ecliptic coordinates of the Earth's sun. + * + * @param jd Julian date. + * @return Ecliptic rectangular geocentric coordinates of the Earth's sun, with distance in AU. + */ +double3 approx_sun_ecliptic(double jd); + +/** + * Approximates the ecliptic coordinates of the Earth's moon. + * + * @param jd Julian date. + * @return Ecliptic rectangular geocentric coordinates of the Earth's moon, with distance in Earth radii. + */ +double3 approx_moon_ecliptic(double jd); + +/** + * Approximates the ecliptic rotation of the Earth's moon. + * + * @param jd Julian date. + * @return Rotation matrix representing the moon's rotation in ecliptic space. + */ +double3x3 approx_moon_ecliptic_rotation(double jd); + struct orbital_elements { double a; ///< Semi-major axis (km) diff --git a/src/game/astronomy/celestial-time.cpp b/src/game/astronomy/celestial-time.cpp new file mode 100644 index 0000000..c67a23f --- /dev/null +++ b/src/game/astronomy/celestial-time.cpp @@ -0,0 +1,62 @@ +/* + * Copyright (C) 2020 Christopher J. Howard + * + * This file is part of Antkeeper source code. + * + * Antkeeper source code is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * Antkeeper source code is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with Antkeeper source code. If not, see . + */ + +#include "celestial-time.hpp" +#include "math/angles.hpp" +#include + +namespace ast +{ + +double ut_to_jd(int year, int month, int day, int hour, int minute, double second) +{ + if (month < 3) + { + month += 12; + year -= 1; + } + + const signed long Y = year; + const signed long M = month; + const signed long D = day; + const signed long JDN = (1461*(Y+4800+(M-14)/12))/4+(367*(M-2-12*((M-14)/12)))/12-(3*((Y+4900+(M-14)/12)/100))/4+D-32075; + + const double h = static_cast(hour - 12) / 24.0; + const double m = static_cast(minute) / 1440.0; + const double s = second / 86400.0; + + return static_cast(JDN) + h + m + s; +} + +double jd_to_gmst(double jd) +{ + return math::wrap_radians(4.894961212 + 6.300388098 * (jd - 2451545.0)); +} + +double jd_to_lmst(double jd, double longitude) +{ + return gmst_to_lmst(jd_to_gmst(jd), longitude); +} + +double gmst_to_lmst(double gmst, double longitude) +{ + return gmst + longitude; +} + +} // namespace ast diff --git a/src/game/astronomy/celestial-time.hpp b/src/game/astronomy/celestial-time.hpp new file mode 100644 index 0000000..6f3a179 --- /dev/null +++ b/src/game/astronomy/celestial-time.hpp @@ -0,0 +1,66 @@ +/* + * Copyright (C) 2020 Christopher J. Howard + * + * This file is part of Antkeeper source code. + * + * Antkeeper source code is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * Antkeeper source code is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with Antkeeper source code. If not, see . + */ + +#ifndef ANTKEEPER_CELESTIAL_TIME_HPP +#define ANTKEEPER_CELESTIAL_TIME_HPP + +namespace ast +{ + +/** + * Calculates the Julian date from a universal time. + * + * @param year Gregorian year, with 1 BC as 0. + * @param month Gregorian month, on [1, 12]. + * @param day Gregorian day, on [1, 31]. + * @param hour Hour, on [0, 23]. + * @param minute Minute, on [0, 59]. + * @param second Second on [0, 60). + * @return Julian date with fractional day. + */ +double ut_to_jd(int year, int month, int day, int hour, int minute, double second); + +/** + * Calculates the Greenwich mean sidereal time (GMST) from a Julian date. + * + * @param jd Julian date. + * @return GMST, in radians. + */ +double jd_to_gmst(double jd); + +/** + * Calculates local mean sidereal time (LMST) from a Julian date. + * + * @param jd Julian date. + * @param longitude Longitude of the observer, in radians. + */ +double jd_to_lmst(double jd, double longitude); + +/** + * Calculates local mean sidereal time (LMST) from Greenwich mean sidereal time (GMST). + * + * @param gmst Greenwich mean sidereal time, in radians. + * @param longitude Longitude of the observer, in radians. + * @return Local mean sidereal time, in radians. + */ +double gmst_to_lmst(double gmst, double longitude); + +} // namespace ast + +#endif // ANTKEEPER_CELESTIAL_TIME_HPP diff --git a/src/game/states/play-state.cpp b/src/game/states/play-state.cpp index 5a755b4..c3c52e0 100644 --- a/src/game/states/play-state.cpp +++ b/src/game/states/play-state.cpp @@ -81,7 +81,7 @@ void play_state_enter(game_context* ctx) sky_pass->set_moon_model(ctx->resource_manager->load("moon.mdl")); ctx->weather_system->set_location(ctx->biome->location[0], ctx->biome->location[1], ctx->biome->location[2]); - ctx->weather_system->set_time(2017, 8, 21, 5, 0, 0, -7.0); + ctx->weather_system->set_time(2017, 8, 21, 5, 0, 0.0, -7.0); ctx->weather_system->set_sky_palette(ctx->biome->sky_palette); ctx->weather_system->set_sun_palette(ctx->biome->sun_palette); ctx->weather_system->set_ambient_palette(ctx->biome->ambient_palette); @@ -118,7 +118,7 @@ void play_state_enter(game_context* ctx) marker_archetype->assign(ecs_registry, ctx->marker_entity); container_archetype->assign(ecs_registry, ctx->container_entity); twig_archetype->assign(ecs_registry, ctx->twig_entity); - + // Create flashlight and light cone, set light cone parent to flashlight, and move both to underworld scene flashlight_archetype->assign(ecs_registry, ctx->flashlight_entity); auto flashlight_light_cone = flashlight_light_cone_archetype->create(ecs_registry); @@ -137,7 +137,7 @@ void play_state_enter(game_context* ctx) //ec::bind_transform(ecs_registry, lens_light_cone, ctx->lens_entity); ec::parent(ecs_registry, lens_light_cone, ctx->lens_entity); - + // Hide inactive tools ec::assign_render_layers(ecs_registry, ctx->forceps_entity, 0); @@ -241,19 +241,19 @@ void play_state_enter(game_context* ctx) ecs_registry.assign_or_replace(ctx->focal_point_entity, focal_point_transform); ecs_registry.assign_or_replace(ctx->focal_point_entity, focal_point_follow); ecs_registry.assign_or_replace(ctx->focal_point_entity, focal_point_snap); - + // Setup camera ctx->overworld_camera->look_at({0, 0, 1}, {0, 0, 0}, {0, 1, 0}); ctx->camera_system->set_camera(ctx->overworld_camera); - + auto ant_head = ant_head_archetype->create(ecs_registry); ec::place(ecs_registry, ant_head, {50, 0}); - + ctx->overworld_scene->update_tweens(); - + // Allocate a nest nest* nest = new ::nest(); - + // Setup initial nest parameters float tunnel_radius = 1.15f; nest->set_tunnel_radius(tunnel_radius); @@ -275,7 +275,7 @@ void play_state_enter(game_context* ctx) chamber.outer_radius = 10.0f; central_shaft->chambers.push_back(chamber); } - + // Dig nest shafts float shift = 0.1f; for (int i = 0; i < 800; ++i) diff --git a/src/game/systems/weather-system.cpp b/src/game/systems/weather-system.cpp index 1bbeae9..f49601c 100644 --- a/src/game/systems/weather-system.cpp +++ b/src/game/systems/weather-system.cpp @@ -26,6 +26,8 @@ #include "utility/gamma.hpp" #include "resources/image.hpp" #include "game/astronomy/celestial-coordinates.hpp" +#include "game/astronomy/celestial-mechanics.hpp" +#include "game/astronomy/celestial-time.hpp" #include #include @@ -33,190 +35,6 @@ static constexpr double hours_per_day = 24.0; static constexpr double minutes_per_day = hours_per_day * 60.0; static constexpr double seconds_per_day = minutes_per_day * 60.0; -/** - * - * @param year Gregorian year - * @param month Month (1 = January, 12 = December) - * @param day Day (1-31) - * @param time Universal time in decimal hours. - */ -static double julian_day(int year, int month, int day, double time) -{ - if (month < 3) - { - month += 12; - year -= 1; - } - - double y = static_cast(year); - double m = static_cast(month); - double d = static_cast(day); - - return std::floor(365.25 * y) + std::floor(30.6001 * (m + 1.0)) - 15.0 + 1720996.5 + d + time; -} - -/** - * Calculates the ecliptic rectangular geocentric coordinates of the sun, with distance in AU. - */ -double3 calculate_sun_ecliptic(double jd) -{ - const double t = (jd - 2451545.0) / 36525.0; - const double m = 6.24 + 628.302 * t; - - const double longitude = 4.895048 + 628.331951 * t + (0.033417 - 0.000084 * t) * std::sin(m) + 0.000351 * std::sin(m * 2.0); - const double latitude = 0.0; - const double distance = 1.000140 - (0.016708 - 0.000042 * t) * std::cos(m) - 0.000141 * std::cos(m * 2.0); - - double3 ecliptic; - ecliptic.x = distance * std::cos(longitude) * std::cos(latitude); - ecliptic.y = distance * std::sin(longitude) * std::cos(latitude); - ecliptic.z = distance * std::sin(latitude); - - return ecliptic; -} - -/** - * Calculates the ecliptic rectangular geocentric coordinates of the moon, with distance in Earth radii. - */ -double3 calculate_moon_ecliptic(double jd) -{ - const double t = (jd - 2451545.0) / 36525.0; - const double l1 = 3.8104 + 8399.7091 * t; - const double m1 = 2.3554 + 8328.6911 * t; - const double m = 6.2300 + 628.3019 * t; - const double d = 5.1985 + 7771.3772 * t; - const double d2 = d * 2.0; - const double f = 1.6280 + 8433.4663 * t; - - const double longitude = l1 - + 0.1098 * std::sin(m1) - + 0.0222 * std::sin(d2 - m1) - + 0.0115 * std::sin(d2) - + 0.0037 * std::sin(m1 * 2.0) - - 0.0032 * std::sin(m) - - 0.0020 * std::sin(d2) - + 0.0010 * std::sin(d2 - m1 * 2.0) - + 0.0010 * std::sin(d2 - m - m1) - + 0.0009 * std::sin(d2 + m1) - + 0.0008 * std::sin(d2 - m) - + 0.0007 * std::sin(m1 - m) - - 0.0006 * std::sin(d) - - 0.0005 * std::sin(m + m1); - - const double latitude = 0.0895 * sin(f) - + 0.0049 * std::sin(m1 + f) - + 0.0048 * std::sin(m1 - f) - + 0.0030 * std::sin(d2 - f) - + 0.0010 * std::sin(d2 + f - m1) - + 0.0008 * std::sin(d2 - f - m1) - + 0.0006 * std::sin(d2 + f); - - const double r = 1.0 / (0.016593 - + 0.000904 * std::cos(m1) - + 0.000166 * std::cos(d2 - m1) - + 0.000137 * std::cos(d2) - + 0.000049 * std::cos(m1 * 2.0) - + 0.000015 * std::cos(d2 + m1) - + 0.000009 * std::cos(d2 - m)); - - double3 ecliptic; - ecliptic.x = r * std::cos(longitude) * std::cos(latitude); - ecliptic.y = r * std::sin(longitude) * std::cos(latitude); - ecliptic.z = r * std::sin(latitude); - - return ecliptic; -} - -double3x3 find_moon_ecliptic_rotation(double jd) -{ - const double t = (jd - 2451545.0) / 36525.0; - const double l1 = 3.8104 + 8399.7091 * t; - const double f = 1.6280 + 8433.4663 * t; - - const double az0 = f + math::pi; - const double ax = 0.026920; - const double az1 = l1 - f; - - double3x3 rz0 = - { - cos(az0), -sin(az0), 0, - sin(az0), cos(az0), 0, - 0, 0, 1 - }; - - double3x3 rx = - { - 1, 0, 0, - 0, cos(ax), -sin(ax), - 0, sin(ax), cos(ax) - }; - - double3x3 rz1 = - { - cos(az1), -sin(az1), 0, - sin(az1), cos(az1), 0, - 0, 0, 1 - }; - - return rz0 * rx * rz1; -} - -/// @see http://www.stjarnhimlen.se/comp/ppcomp.html -/// @see http://www.geoastro.de/elevazmoon/basics/index.htm -void ecliptic_to_equatorial(double longitude, double latitude, double ecl, double* right_ascension, double* declination) -{ - double eclip_x = std::cos(longitude) * std::cos(latitude); - double eclip_y = std::sin(longitude) * std::cos(latitude); - double eclip_z = std::sin(latitude); - - double equat_x = eclip_x; - double equat_y = eclip_y * std::cos(ecl) - eclip_z * std::sin(ecl); - double equat_z = eclip_y * std::sin(ecl) + eclip_z * std::cos(ecl); - - *right_ascension = std::atan2(equat_y, equat_x); - *declination = std::atan2(equat_z, sqrt(equat_x * equat_x + equat_y * equat_y)); -} - -/// @see http://www.stjarnhimlen.se/comp/ppcomp.html -/// @see http://www.geoastro.de/elevazmoon/basics/index.htm -void equatorial_to_horizontal(double right_ascension, double declination, double lmst, double latitude, double* azimuth, double* elevation) -{ - double hour_angle = lmst - right_ascension; - - double x = std::cos(hour_angle) * std::cos(declination); - double y = std::sin(hour_angle) * std::cos(declination); - double z = std::sin(declination); - - double horiz_x = x * std::cos(math::half_pi - latitude) - z * std::sin(math::half_pi - latitude); - double horiz_y = y; - double horiz_z = x * std::sin(math::half_pi - latitude) + z * std::cos(math::half_pi - latitude); - - *azimuth = math::wrap_radians(std::atan2(horiz_y, horiz_x) + math::pi); - *elevation = math::wrap_radians(std::atan2(horiz_z, std::sqrt(horiz_x * horiz_x + horiz_y * horiz_y))); -} - - -double3x3 horizontal_to_right_handed() -{ - return double3x3 - { - 0.0, 0.0, 1.0, - 1.0, 0.0, 0.0, - 0.0, -1.0, 0.0 - }; -} - -/** - * Calculates the Greenwich mean sidereal time (GMST) from a Julian day. - * - * @param jd Julian day. - * @return GMST, in radians. - */ -static double jd_to_gmst(double jd) -{ - return math::wrap_radians(4.894961212 + 6.300388098 * (jd - 2451545.0)); -} - weather_system::weather_system(entt::registry& registry): entity_system(registry), ambient_light(nullptr), @@ -239,42 +57,38 @@ void weather_system::update(double t, double dt) const float latitude = location[0]; const float longitude = location[1]; - // Calculate local mean sidereal time (LMST) + // Calculate local time double time_correction = longitude / (math::two_pi / 24.0); double local_jd = jd + time_correction / 24.0 - 0.5; double local_time = (local_jd - std::floor(local_jd)) * 24.0; - double hour = local_time; - double gmst = jd_to_gmst(jd); - double lmst = gmst + longitude; - - // Obliquity of the ecliptic - double ecl = math::radians(23.4393 - 3.563e-7 * (jd - 2451545.0)); // Solar distance in AU //double sr = ... // Apparent radius in degrees //double sradius = 0.2666 / sr; + double lmst = ast::jd_to_lmst(jd, longitude); + double ecl = ast::approx_ecliptic_obliquity(jd); double3x3 ecliptic_to_horizontal = ast::ecliptic_to_horizontal(ecl, latitude, lmst); + - double3 sun_ecliptic = calculate_sun_ecliptic(jd); + double3 sun_ecliptic = ast::approx_sun_ecliptic(jd); double3 sun_horizontal = ecliptic_to_horizontal * sun_ecliptic; + sun_horizontal.z -= 4.25875e-5; // Subtract one earth radius (in AU), for positon of observer double3 sun_spherical = ast::rectangular_to_spherical(sun_horizontal); - double3 sun_positiond = horizontal_to_right_handed() * sun_horizontal; + double3 sun_positiond = ast::horizontal_to_right_handed * sun_horizontal; float2 sun_az_el = {static_cast(sun_spherical.z) - math::pi, static_cast(sun_spherical.y)}; float3 sun_position = math::normalize(float3{static_cast(sun_positiond.x), static_cast(sun_positiond.y), static_cast(sun_positiond.z)}); - double3 moon_ecliptic = calculate_moon_ecliptic(jd); + double3 moon_ecliptic = ast::approx_moon_ecliptic(jd); double3 moon_horizontal = ecliptic_to_horizontal * moon_ecliptic; moon_horizontal.z -= 1.0; // Subtract one earth radius, for position of observer double3 moon_spherical = ast::rectangular_to_spherical(moon_horizontal); - double3 moon_positiond = horizontal_to_right_handed() * moon_horizontal; + double3 moon_positiond = ast::horizontal_to_right_handed * moon_horizontal; float2 moon_az_el = {static_cast(moon_spherical.z) - math::pi, static_cast(moon_spherical.y)}; float3 moon_position = math::normalize(float3{static_cast(moon_positiond.x), static_cast(moon_positiond.y), static_cast(moon_positiond.z)}); - - //std::cout << "new moon: " << math::degrees(moon_az_el[0]) << ", " << math::degrees(moon_az_el[1]) << std::endl; - - double3x3 moon_rotation_matrix = horizontal_to_right_handed() * ecliptic_to_horizontal; + + double3x3 moon_rotation_matrix = ast::horizontal_to_right_handed * ecliptic_to_horizontal; math::quaternion moon_rotationd = math::normalize(math::quaternion_cast(moon_rotation_matrix) * math::angle_axis(math::half_pi, double3{0, 1, 0}) * math::angle_axis(-math::half_pi, double3{0, 0, -1})); math::quaternion moon_rotation = { @@ -300,9 +114,6 @@ void weather_system::update(double t, double dt) moon_light->set_rotation(moon_az_el_rotation); } - std::size_t hour_index = static_cast(hour); - float lerp_factor = hour - std::floor(hour); - float sun_gradient_position = static_cast(std::max(0.0, ((sun_az_el[1] + math::half_pi) / math::pi))); float moon_gradient_position = static_cast(std::max(0.0, ((moon_az_el[1] + math::half_pi) / math::pi))); float sky_gradient_position = sun_gradient_position; @@ -325,7 +136,7 @@ void weather_system::update(double t, double dt) sky_pass->set_horizon_color(horizon_color); sky_pass->set_zenith_color(zenith_color); - sky_pass->set_time_of_day(static_cast(hour * 60.0 * 60.0)); + sky_pass->set_time_of_day(static_cast(local_time * 60.0 * 60.0)); sky_pass->set_observer_location(location[0], location[1], location[2]); sky_pass->set_sun_coordinates(sun_position, sun_az_el); sky_pass->set_moon_coordinates(moon_position, moon_az_el); @@ -399,10 +210,9 @@ void weather_system::set_material_pass(::material_pass* pass) material_pass = pass; } -void weather_system::set_time(int year, int month, int day, int hour, int minute, int second, double tc) +void weather_system::set_time(int year, int month, int day, int hour, int minute, double second, double tc) { - double time = ((static_cast(hour) - tc) + ((static_cast(minute) + static_cast(second) / 60.0) / 60.0)) / 24.0; - jd = julian_day(year, month, day, time); + jd = ast::ut_to_jd(year, month, day, hour, minute, second) - tc / 24.0; } void weather_system::set_time_scale(float scale) diff --git a/src/game/systems/weather-system.hpp b/src/game/systems/weather-system.hpp index dd54cd9..d352389 100644 --- a/src/game/systems/weather-system.hpp +++ b/src/game/systems/weather-system.hpp @@ -53,7 +53,7 @@ public: void set_material_pass(::material_pass* pass); /// @param tc Timezone correction, in hours - void set_time(int year, int month, int day, int hour, int minute, int second, double tc); + void set_time(int year, int month, int day, int hour, int minute, double second, double tc); void set_time_scale(float scale); void set_sky_palette(const ::image* image);