From 9b3eff2ee85ca6627342cfbbe1ac3ba988dc377f Mon Sep 17 00:00:00 2001 From: Dan Goodliffe Date: Fri, 20 Sep 2024 20:27:51 +0100 Subject: Move getSunPos to Environment --- game/environment.cpp | 67 +++++++++++++++++++++++++++++++++++++++++++++ game/environment.h | 2 ++ test/test-environment.cpp | 70 ++--------------------------------------------- 3 files changed, 71 insertions(+), 68 deletions(-) diff --git a/game/environment.cpp b/game/environment.cpp index fd2bfd4..665c11b 100644 --- a/game/environment.cpp +++ b/game/environment.cpp @@ -16,3 +16,70 @@ Environment::render(const SceneRenderer & renderer, const SceneProvider & scene) renderer.setAmbientLight({0.5F, 0.5F, 0.5F}); renderer.setDirectionalLight({0.6F, 0.6F, 0.6F}, {-1, 1, -1}, scene); } + +// Based on the C++ code published at https://www.psa.es/sdg/sunpos.htm +// Linked from https://www.pveducation.org/pvcdrom/properties-of-sunlight/suns-position-to-high-accuracy +Direction2D +Environment::getSunPos(const Direction2D position, const time_t time) +{ + auto & longitude = position.x; + auto & latitude = position.y; + using std::acos; + using std::asin; + using std::atan2; + using std::cos; + using std::floor; + using std::sin; + using std::tan; + static const auto JD2451545 = "2000-01-01T12:00:00"_time_t; + + // Calculate difference in days between the current Julian Day + // and JD 2451545.0, which is noon 1 January 2000 Universal Time + // Calculate time of the day in UT decimal hours + const auto dDecimalHours = static_cast(time % 86400) / 3600.F; + const auto dElapsedJulianDays = static_cast(time - JD2451545) / 86400.F; + + // Calculate ecliptic coordinates (ecliptic longitude and obliquity of the + // ecliptic in radians but without limiting the angle to be less than 2*Pi + // (i.e., the result may be greater than 2*Pi) + const auto dOmega = 2.1429F - 0.0010394594F * dElapsedJulianDays; + const auto dMeanLongitude = 4.8950630F + 0.017202791698F * dElapsedJulianDays; // Radians + const auto dMeanAnomaly = 6.2400600F + 0.0172019699F * dElapsedJulianDays; + const auto dEclipticLongitude = dMeanLongitude + 0.03341607F * sin(dMeanAnomaly) + + 0.00034894F * sin(2 * dMeanAnomaly) - 0.0001134F - 0.0000203F * sin(dOmega); + const auto dEclipticObliquity = 0.4090928F - 6.2140e-9F * dElapsedJulianDays + 0.0000396F * cos(dOmega); + + // Calculate celestial coordinates ( right ascension and declination ) in radians + // but without limiting the angle to be less than 2*Pi (i.e., the result may be + // greater than 2*Pi) + const auto dSin_EclipticLongitude = sin(dEclipticLongitude); + const auto dY = cos(dEclipticObliquity) * dSin_EclipticLongitude; + const auto dX = cos(dEclipticLongitude); + auto dRightAscension = atan2(dY, dX); + if (dRightAscension < 0) { + dRightAscension = dRightAscension + two_pi; + } + const auto dDeclination = asin(sin(dEclipticObliquity) * dSin_EclipticLongitude); + + // Calculate local coordinates ( azimuth and zenith angle ) in degrees + const auto dGreenwichMeanSiderealTime = 6.6974243242F + 0.0657098283F * dElapsedJulianDays + dDecimalHours; + const auto dLocalMeanSiderealTime + = (dGreenwichMeanSiderealTime * 15.0F + (longitude / degreesToRads)) * degreesToRads; + const auto dHourAngle = dLocalMeanSiderealTime - dRightAscension; + const auto dLatitudeInRadians = latitude; + const auto dCos_Latitude = cos(dLatitudeInRadians); + const auto dSin_Latitude = sin(dLatitudeInRadians); + const auto dCos_HourAngle = cos(dHourAngle); + Direction2D udtSunCoordinates; + udtSunCoordinates.y + = (acos(dCos_Latitude * dCos_HourAngle * cos(dDeclination) + sin(dDeclination) * dSin_Latitude)); + udtSunCoordinates.x = atan2(-sin(dHourAngle), tan(dDeclination) * dCos_Latitude - dSin_Latitude * dCos_HourAngle); + if (udtSunCoordinates.x < 0) { + udtSunCoordinates.x = udtSunCoordinates.x + two_pi; + } + // Parallax Correction + const auto dParallax = (earthMeanRadius / astronomicalUnit) * sin(udtSunCoordinates.y); + udtSunCoordinates.y = half_pi - (udtSunCoordinates.y + dParallax); + + return udtSunCoordinates; +} diff --git a/game/environment.h b/game/environment.h index 62eedb8..a6f3036 100644 --- a/game/environment.h +++ b/game/environment.h @@ -1,5 +1,6 @@ #pragma once +#include "config/types.h" #include "worldobject.h" class SceneRenderer; @@ -10,6 +11,7 @@ public: Environment(); void tick(TickDuration elapsed) override; void render(const SceneRenderer &, const SceneProvider &) const; + static Direction2D getSunPos(const Direction2D position, const time_t time); private: time_t worldTime; diff --git a/test/test-environment.cpp b/test/test-environment.cpp index 0e9e37a..b6e0e4f 100644 --- a/test/test-environment.cpp +++ b/test/test-environment.cpp @@ -6,75 +6,9 @@ #include #include +#include #include -// Based on the C++ code published at https://www.psa.es/sdg/sunpos.htm -// Linked from https://www.pveducation.org/pvcdrom/properties-of-sunlight/suns-position-to-high-accuracy -Direction2D -getSunPos(const Direction2D position, const time_t time) -{ - auto & longitude = position.x; - auto & latitude = position.y; - using std::acos; - using std::asin; - using std::atan2; - using std::cos; - using std::floor; - using std::sin; - using std::tan; - static const auto JD2451545 = "2000-01-01T12:00:00"_time_t; - - // Calculate difference in days between the current Julian Day - // and JD 2451545.0, which is noon 1 January 2000 Universal Time - // Calculate time of the day in UT decimal hours - const auto dDecimalHours = static_cast(time % 86400) / 3600.F; - const auto dElapsedJulianDays = static_cast(time - JD2451545) / 86400.F; - - // Calculate ecliptic coordinates (ecliptic longitude and obliquity of the - // ecliptic in radians but without limiting the angle to be less than 2*Pi - // (i.e., the result may be greater than 2*Pi) - const auto dOmega = 2.1429F - 0.0010394594F * dElapsedJulianDays; - const auto dMeanLongitude = 4.8950630F + 0.017202791698F * dElapsedJulianDays; // Radians - const auto dMeanAnomaly = 6.2400600F + 0.0172019699F * dElapsedJulianDays; - const auto dEclipticLongitude = dMeanLongitude + 0.03341607F * sin(dMeanAnomaly) - + 0.00034894F * sin(2 * dMeanAnomaly) - 0.0001134F - 0.0000203F * sin(dOmega); - const auto dEclipticObliquity = 0.4090928F - 6.2140e-9F * dElapsedJulianDays + 0.0000396F * cos(dOmega); - - // Calculate celestial coordinates ( right ascension and declination ) in radians - // but without limiting the angle to be less than 2*Pi (i.e., the result may be - // greater than 2*Pi) - const auto dSin_EclipticLongitude = sin(dEclipticLongitude); - const auto dY = cos(dEclipticObliquity) * dSin_EclipticLongitude; - const auto dX = cos(dEclipticLongitude); - auto dRightAscension = atan2(dY, dX); - if (dRightAscension < 0) { - dRightAscension = dRightAscension + two_pi; - } - const auto dDeclination = asin(sin(dEclipticObliquity) * dSin_EclipticLongitude); - - // Calculate local coordinates ( azimuth and zenith angle ) in degrees - const auto dGreenwichMeanSiderealTime = 6.6974243242F + 0.0657098283F * dElapsedJulianDays + dDecimalHours; - const auto dLocalMeanSiderealTime - = (dGreenwichMeanSiderealTime * 15.0F + (longitude / degreesToRads)) * degreesToRads; - const auto dHourAngle = dLocalMeanSiderealTime - dRightAscension; - const auto dLatitudeInRadians = latitude; - const auto dCos_Latitude = cos(dLatitudeInRadians); - const auto dSin_Latitude = sin(dLatitudeInRadians); - const auto dCos_HourAngle = cos(dHourAngle); - Direction2D udtSunCoordinates; - udtSunCoordinates.y - = (acos(dCos_Latitude * dCos_HourAngle * cos(dDeclination) + sin(dDeclination) * dSin_Latitude)); - udtSunCoordinates.x = atan2(-sin(dHourAngle), tan(dDeclination) * dCos_Latitude - dSin_Latitude * dCos_HourAngle); - if (udtSunCoordinates.x < 0) { - udtSunCoordinates.x = udtSunCoordinates.x + two_pi; - } - // Parallax Correction - const auto dParallax = (earthMeanRadius / astronomicalUnit) * sin(udtSunCoordinates.y); - udtSunCoordinates.y = half_pi - (udtSunCoordinates.y + dParallax); - - return udtSunCoordinates; -} - using sunPosTestData = std::tuple; constexpr Direction2D Doncaster = {-1.1, 53.5}; constexpr Direction2D NewYork = {74.0, 40.7}; @@ -95,7 +29,7 @@ BOOST_DATA_TEST_CASE(sun_position, }), position, timeOfYear, expSunPos) { - const auto sunPos = getSunPos(position * degreesToRads, timeOfYear) / degreesToRads; + const auto sunPos = Environment::getSunPos(position * degreesToRads, timeOfYear) / degreesToRads; BOOST_CHECK_CLOSE(sunPos.x, expSunPos.x, 1.F); BOOST_CHECK_CLOSE(sunPos.y, expSunPos.y, 1.F); } -- cgit v1.2.3