From c81e93c2eaa731fa4564c60392c95b6ba84e9deb Mon Sep 17 00:00:00 2001 From: John Bytheway <52664+jbytheway@users.noreply.github.com> Date: Mon, 21 Jun 2021 07:01:41 -0400 Subject: [PATCH] Assume solar time rather than timezone Previously the sunlight angle code required you to pass a timezone. Now, instead, it assumes the timezone is local solar time (determined by the longitude). This makes some of the calculations simpler and the API less confusing. It does mean that the test cases have drifted a little further from the 'true' Boston sunrise and sunset times, but they're still fairly close. --- src/calendar.cpp | 50 ++++++++++++++------------ src/calendar.h | 5 +-- tests/sun_test.cpp | 88 +++++++++++++++++++++++----------------------- 3 files changed, 73 insertions(+), 70 deletions(-) diff --git a/src/calendar.cpp b/src/calendar.cpp index 6ccf6e2269b35..e71c63904607c 100644 --- a/src/calendar.cpp +++ b/src/calendar.cpp @@ -180,10 +180,11 @@ static units::angle sidereal_time_at( solar_effective_time t, units::angle longi } static std::pair sun_azimuth_altitude( - solar_effective_time t, lat_long location, time_duration timezone ) + solar_effective_time t, lat_long location ) { units::angle right_ascension; units::angle declination; + time_duration timezone = angle_to_time( location.longitude ); std::tie( right_ascension, declination ) = sun_ra_declination( t, timezone ); const units::angle sidereal_time = sidereal_time_at( t, location.longitude, timezone ); @@ -221,27 +222,26 @@ static std::pair sun_azimuth_altitude( return std::make_pair( azimuth, altitude ); } -std::pair sun_azimuth_altitude( - time_point t, lat_long location, time_duration timezone ) +std::pair sun_azimuth_altitude( time_point t, lat_long location ) { - return sun_azimuth_altitude( solar_effective_time( t ), location, timezone ); + return sun_azimuth_altitude( solar_effective_time( t ), location ); } -static units::angle sun_altitude( time_point t, lat_long location, time_duration timezone ) +static units::angle sun_altitude( time_point t, lat_long location ) { - return sun_azimuth_altitude( t, location, timezone ).second; + return sun_azimuth_altitude( t, location ).second; } static units::angle sun_altitude( time_point t ) { - return sun_altitude( t, location_boston, timezone_boston ); + return sun_altitude( t, location_boston ); } cata::optional sunlight_angle( const time_point &t, lat_long location ) { units::angle azimuth; units::angle altitude; - std::tie( azimuth, altitude ) = sun_azimuth_altitude( t, location, timezone_boston ); + std::tie( azimuth, altitude ) = sun_azimuth_altitude( t, location ); if( altitude <= 0_degrees ) { // Sun below horizon return cata::nullopt; @@ -259,17 +259,22 @@ cata::optional sunlight_angle( const time_point &t ) static time_point solar_noon_near( const time_point &t ) { - constexpr time_duration longitude_hours = angle_to_time( location_boston.longitude ); const time_point prior_midnight = t - time_past_midnight( t ); - return prior_midnight + 12_hours - longitude_hours + timezone_boston; + return prior_midnight + 12_hours; + // If we were using a timezone rather than local solar time this would be + // calculated as follows: + //constexpr time_duration longitude_hours = angle_to_time( location_boston.longitude ); + //return prior_midnight + 12_hours - longitude_hours + timezone; } -static units::angle offset_to_sun_altitude( const units::angle altitude, - const solar_effective_time approx_time, const bool evening ) +static units::angle offset_to_sun_altitude( + const units::angle altitude, const units::angle longitude, + const solar_effective_time approx_time, const bool evening ) { units::angle ra; units::angle declination; - std::tie( ra, declination ) = sun_ra_declination( approx_time, timezone_boston ); + time_duration timezone = angle_to_time( longitude ); + std::tie( ra, declination ) = sun_ra_declination( approx_time, timezone ); double cos_hour_angle = ( sin( altitude ) - sin( location_boston.latitude ) * sin( declination ) ) / cos( location_boston.latitude ) / cos( declination ); @@ -284,16 +289,16 @@ static units::angle offset_to_sun_altitude( const units::angle altitude, } const units::angle target_sidereal_time = hour_angle + ra; const units::angle sidereal_time_at_approx_time = - sidereal_time_at( approx_time, location_boston.longitude, timezone_boston ); + sidereal_time_at( approx_time, location_boston.longitude, timezone ); return normalize( target_sidereal_time - sidereal_time_at_approx_time ); } -static time_point sun_at_altitude( const units::angle altitude, const time_point t, - const bool evening ) +static time_point sun_at_altitude( const units::angle altitude, const units::angle longitude, + const time_point t, const bool evening ) { const time_point solar_noon = solar_noon_near( t ); units::angle initial_offset = - offset_to_sun_altitude( altitude, solar_effective_time( solar_noon ), evening ); + offset_to_sun_altitude( altitude, longitude, solar_effective_time( solar_noon ), evening ); if( !evening ) { initial_offset -= 360_degrees; } @@ -302,7 +307,8 @@ static time_point sun_at_altitude( const units::angle altitude, const time_point // Now we should have the correct time to within a few minutes; iterate to // get a more precise estimate units::angle correction_offset = - offset_to_sun_altitude( altitude, solar_effective_time( initial_approximation ), evening ); + offset_to_sun_altitude( altitude, longitude, solar_effective_time( initial_approximation ), + evening ); if( correction_offset > 180_degrees ) { correction_offset -= 360_degrees; } @@ -312,22 +318,22 @@ static time_point sun_at_altitude( const units::angle altitude, const time_point time_point sunrise( const time_point &p ) { - return sun_at_altitude( 0_degrees, p, false ); + return sun_at_altitude( 0_degrees, location_boston.longitude, p, false ); } time_point sunset( const time_point &p ) { - return sun_at_altitude( 0_degrees, p, true ); + return sun_at_altitude( 0_degrees, location_boston.longitude, p, true ); } time_point night_time( const time_point &p ) { - return sun_at_altitude( min_sun_angle_for_day, p, true ); + return sun_at_altitude( min_sun_angle_for_day, location_boston.longitude, p, true ); } time_point daylight_time( const time_point &p ) { - return sun_at_altitude( min_sun_angle_for_day, p, false ); + return sun_at_altitude( min_sun_angle_for_day, location_boston.longitude, p, false ); } bool is_night( const time_point &p ) diff --git a/src/calendar.h b/src/calendar.h index 96e1e635f4f7a..8a9cde2fab3ba 100644 --- a/src/calendar.h +++ b/src/calendar.h @@ -593,8 +593,7 @@ float sun_moon_light_at( const time_point &p ); /** How much light is provided at the solar noon nearest to given time */ double sun_moon_light_at_noon_near( const time_point &p ); -std::pair sun_azimuth_altitude( - time_point, lat_long, time_duration timezone ); +std::pair sun_azimuth_altitude( time_point, lat_long ); /** Returns the offset by which a ray of sunlight would move when shifting down * one z-level, or nullopt if the sun is below the horizon. @@ -604,8 +603,6 @@ std::pair sun_azimuth_altitude( cata::optional sunlight_angle( const time_point &, lat_long ); cata::optional sunlight_angle( const time_point & ); -constexpr time_duration timezone_boston = -5_hours; - enum class weekdays : int { SUNDAY = 0, MONDAY, diff --git a/tests/sun_test.cpp b/tests/sun_test.cpp index e44a92d8f1dfe..197c94f072137 100644 --- a/tests/sun_test.cpp +++ b/tests/sun_test.cpp @@ -241,26 +241,26 @@ TEST_CASE( "sunrise and sunset", "[sun][sunrise][sunset][equinox][solstice]" ) SECTION( "spring equinox is day 1 of spring" ) { // Actual sunrise and sunset on March 21st 2001 are 0545 and 1757 - CHECK( "Year 1, Spring, day 1 5:47:36 AM" == to_string( sunrise( spring ) ) ); - CHECK( "Year 1, Spring, day 1 5:54:55 PM" == to_string( sunset( spring ) ) ); + CHECK( "Year 1, Spring, day 1 6:03:22 AM" == to_string( sunrise( spring ) ) ); + CHECK( "Year 1, Spring, day 1 6:10:41 PM" == to_string( sunset( spring ) ) ); } SECTION( "summer solstice is day 1 of summer" ) { // Actual sunrise and sunset on June 21st 2001 are 0407 and 1924 - CHECK( "Year 1, Summer, day 1 4:13:02 AM" == to_string( sunrise( summer ) ) ); - CHECK( "Year 1, Summer, day 1 7:19:24 PM" == to_string( sunset( summer ) ) ); + CHECK( "Year 1, Summer, day 1 4:28:48 AM" == to_string( sunrise( summer ) ) ); + CHECK( "Year 1, Summer, day 1 7:35:10 PM" == to_string( sunset( summer ) ) ); } SECTION( "autumn equinox is day 1 of autumn" ) { // Actual sunrise and sunset on September 22nd 2001 are 0531 and 1741 - CHECK( "Year 1, Autumn, day 1 5:35:13 AM" == to_string( sunrise( autumn ) ) ); - CHECK( "Year 1, Autumn, day 1 5:38:27 PM" == to_string( sunset( autumn ) ) ); + CHECK( "Year 1, Autumn, day 1 5:50:59 AM" == to_string( sunrise( autumn ) ) ); + CHECK( "Year 1, Autumn, day 1 5:54:13 PM" == to_string( sunset( autumn ) ) ); } SECTION( "winter solstice is day 1 of winter" ) { // Actual sunrise and sunset on December 21st 2001 are 0710 and 1614 - CHECK( "Year 1, Winter, day 1 7:15:44 AM" == to_string( sunrise( winter ) ) ); - CHECK( "Year 1, Winter, day 1 4:09:37 PM" == to_string( sunset( winter ) ) ); + CHECK( "Year 1, Winter, day 1 7:31:30 AM" == to_string( sunrise( winter ) ) ); + CHECK( "Year 1, Winter, day 1 4:25:23 PM" == to_string( sunset( winter ) ) ); } SECTION( "spring sunrise gets earlier" ) { @@ -412,7 +412,7 @@ TEST_CASE( "sunrise_sunset_consistency", "[sun]" ) units::angle azimuth; units::angle altitude; std::tie( azimuth, altitude ) = - sun_azimuth_altitude( this_sunrise, location_boston, timezone_boston ); + sun_azimuth_altitude( this_sunrise, location_boston ); CHECK( to_degrees( altitude ) == Approx( 0 ).margin( 0.01 ) ); } { @@ -421,7 +421,7 @@ TEST_CASE( "sunrise_sunset_consistency", "[sun]" ) units::angle azimuth; units::angle altitude; std::tie( azimuth, altitude ) = - sun_azimuth_altitude( this_sunset, location_boston, timezone_boston ); + sun_azimuth_altitude( this_sunset, location_boston ); CHECK( to_degrees( altitude ) == Approx( 0 ).margin( 0.01 ) ); } { @@ -430,7 +430,7 @@ TEST_CASE( "sunrise_sunset_consistency", "[sun]" ) units::angle azimuth; units::angle altitude; std::tie( azimuth, altitude ) = - sun_azimuth_altitude( this_daylight, location_boston, timezone_boston ); + sun_azimuth_altitude( this_daylight, location_boston ); CHECK( to_degrees( altitude ) == Approx( -12 ).margin( 0.01 ) ); } } @@ -448,7 +448,7 @@ static PointSet sun_positions_regular( time_point start, time_point end, time_du CAPTURE( to_minutes( t - start ) ); units::angle azimuth; units::angle altitude; - std::tie( azimuth, altitude ) = sun_azimuth_altitude( t, location_boston, timezone_boston ); + std::tie( azimuth, altitude ) = sun_azimuth_altitude( t, location_boston ); if( altitude < 0_degrees ) { continue; } @@ -487,7 +487,7 @@ static PointSet sun_throughout_year( time_point day_start ) static void check_sun_plot( const std::vector &points, const std::string &reference ) { - static constexpr std::array symbols = { { '#', '*', '.' } }; + static constexpr std::array symbols = { { '#', '*', '-' } }; REQUIRE( points.size() <= symbols.size() ); std::ostringstream os; @@ -537,30 +537,30 @@ TEST_CASE( "movement_of_sun_through_day", "[sun]" ) " \n" " \n" " \n" -" ############### \n" +" ################ \n" " #### #### \n" " #### ## \n" " ## ## \n" -" ## ## \n" +" ## ### \n" " ## ## \n" " ## ## \n" -" ## ****** ## \n" +" ## ******* ## \n" " ## **** ***** ## \n" " # *** *** ## \n" -" # *** ** # \n" +" # *** *** # \n" " ## ** ** ## \n" " ## ** ** ## \n" -" # ** ** # \n" -" ## * * # \n" -" ## * .... * # \n" // NOLINT -" # ** ..... ..... * # \n" // NOLINT -" ## ** ... ... ** ## \n" // NOLINT -" # ** ... .. ** # \n" // NOLINT -" # * .. .. * # \n" -" ## * .. .. * ## \n" -" # ** .. .. ** ## \n" -" # ** .. .. ** # \n" -" ## * .. .. * # \n" +" # ** * # \n" +" ## * * ## \n" +" ## * ---- * # \n" +" # ** ----- ----- * # \n" +" ## ** --- --- * ## \n" +" ## ** --- -- ** ## \n" +" # * -- -- ** # \n" +" ## * -- -- * ## \n" +" # ** -- -- ** # \n" +" # ** -- -- ** # \n" +" ## * -- -- * # \n" " Azimuth\n"; // *INDENT-ON* check_sun_plot( { summer_points, equinox_points, winter_points }, reference ); @@ -580,22 +580,22 @@ TEST_CASE( "movement_of_noon_through_year", "[sun]" ) " \n" " \n" " \n" -" ###### \n" -" ## # \n" -" # ## \n" -" # ## \n" -" ### \n" -" ### \n" -" ## # \n" -" ## ## \n" -" ## # \n" -" ## ## \n" -" ## # \n" -" # # \n" -" # # \n" -" # ## \n" -" ## ## \n" -" ####### \n" +" ###### \n" +" ## ## \n" +" ## ## \n" +" ## ## \n" +" ### \n" +" ### \n" +" # ## \n" +" # ## \n" +" ## ## \n" +" # ## \n" +" ## # \n" +" # # \n" +" # # \n" +" ## # \n" +" ## ## \n" +" ####### \n" " \n" " \n" " \n"