Skip to content

Commit

Permalink
Assume solar time rather than timezone
Browse files Browse the repository at this point in the history 
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.
  • Loading branch information
@jbytheway
jbytheway committed Jul 7, 2021
1 parent a19139c commit c81e93c
Show file tree
Hide file tree
Showing 3 changed files with 73 additions and 70 deletions.
50 changes: 28 additions & 22 deletions src/calendar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -180,10 +180,11 @@ static units::angle sidereal_time_at( solar_effective_time t, units::angle longi
}

static std::pair<units::angle, units::angle> 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 );

Expand Down Expand Up @@ -221,27 +222,26 @@ static std::pair<units::angle, units::angle> sun_azimuth_altitude(
return std::make_pair( azimuth, altitude );
}

std::pair<units::angle, units::angle> sun_azimuth_altitude(
time_point t, lat_long location, time_duration timezone )
std::pair<units::angle, units::angle> 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<rl_vec2d> 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;
Expand All @@ -259,17 +259,22 @@ cata::optional<rl_vec2d> 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 );
Expand All @@ -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;
}
Expand All @@ -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;
}
Expand All @@ -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 )
Expand Down
5 changes: 1 addition & 4 deletions src/calendar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<units::angle, units::angle> sun_azimuth_altitude(
time_point, lat_long, time_duration timezone );
std::pair<units::angle, units::angle> 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.
Expand All @@ -604,8 +603,6 @@ std::pair<units::angle, units::angle> sun_azimuth_altitude(
cata::optional<rl_vec2d> sunlight_angle( const time_point &, lat_long );
cata::optional<rl_vec2d> sunlight_angle( const time_point & );

constexpr time_duration timezone_boston = -5_hours;

enum class weekdays : int {
SUNDAY = 0,
MONDAY,
Expand Down
88 changes: 44 additions & 44 deletions tests/sun_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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" ) {
Expand Down Expand Up @@ -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 ) );
}
{
Expand All @@ -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 ) );
}
{
Expand All @@ -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 ) );
}
}
Expand All @@ -448,7 +448,7 @@ static PointSet sun_positions_regular( time_point start, time_point end, time_du
CAPTURE( to_minutes<int>( 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;
}
Expand Down Expand Up @@ -487,7 +487,7 @@ static PointSet sun_throughout_year( time_point day_start )

static void check_sun_plot( const std::vector<PointSet> &points, const std::string &reference )
{
static constexpr std::array<char, 3> symbols = { { '#', '*', '.' } };
static constexpr std::array<char, 3> symbols = { { '#', '*', '-' } };
REQUIRE( points.size() <= symbols.size() );

std::ostringstream os;
Expand Down Expand Up @@ -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 );
Expand All @@ -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"
Expand Down

0 comments on commit c81e93c

Please sign in to comment.