Sunlight angle (CleverRaven#48090)

* Make turn_zero constexpr

* Add features to rl_vec3d

* Construct from rl_vec2d and z-coord
* Extract xy as rl_vec2d
* In-place scaling

* StringMaker support for rl_vec2d

* Add StringMaker specialization for time_point

* Add atan2( rl_vec2d ) overload

* Add units::asin returning units::angle

* Add sunlight angle calculation and tests

This is an initial implementation of sunlight_angle, which should work
OK for the first few days, but will break down before the end of the
first year.

One test plots the path of the sun through the sky at each solstice and
an equinox.  This is intended as a simple smell test that's easy for a
human to see when something has gone wrong with the sun angle algorithm.

* Add analemma unit test

* Tweak sunlight calculation away from reality

Because the game years don't match real-world years, we need to make
some adjustments to the sunlight angle calculations to allow for that.

* Use logical origin of sidereal time

We had an origin of 177, but 180 makes logical sense, so switch to that
as it's easier to understand.

* Derive sidereal time from angle_per_day

Removing one magic constant.

* mean_long, mean_anomaly in terms of angle_per_day

Remove more arbitrary constants from the calculation, and document those
which remain.

* Use game year length, not real-world one

This removes one more arbitrary constant and means that the sunlight
angle calculations ought to work for other season lengths.

* Unit test for handling season length

* Base sunrise and sunset on astronomical maths

This is the first change that allows the sun angle calculations to
actually influence the game.

We switch the sunrise and sunset times to match the sun positions as
calculated by the astronomical equations surrounding the location of the
sun.

This requires changing a few other places in the code, and updating a
lot of tests.

Heavinly inspired by Hirmuolio's work in CleverRaven#47570, but reimplemented
rather than cherry-picked.

Co-authored-by: Hirmuolio <[email protected]>

* Support eternal season for new Sun calculations

When eternal season is set, always have the Sun behave as it does on
whatever day the game starts.

* Update some comments in sun_test.cpp

* Appease clang 3.8

* Improve comments and variable names

Based on Github review.

Co-authored-by: actual-nh <[email protected]>

* Fix clang-tidy issues

* 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.

Co-authored-by: Hirmuolio <[email protected]>
Co-authored-by: actual-nh <[email protected]>

src/calendar.h

@@ -6,12 +6,23 @@
 #include <utility>
 #include <vector>
 
+#include "units_fwd.h"
+
 class JsonIn;
 class JsonOut;
+struct lat_long;
+struct rl_vec2d;
 class time_duration;
 class time_point;
 template<typename T> struct enum_traits;
 
+namespace cata
+{
+template<typename T>
+class optional;
+} // namespace cata
+
+
 /** Real world seasons */
 enum season_type {
     SPRING = 0,
@@ -109,18 +120,6 @@ extern time_point start_of_game;
 extern time_point turn;
 extern season_type initial_season;
 
-/**
- * A time point that is always before the current turn, even when the game has
- * just started. This implies `before_time_starts < calendar::turn` is always
- * true. It can be used to initialize `time_point` values that denote that last
- * time a cache was update.
- */
-extern const time_point before_time_starts;
-/**
- * Represents time point 0.
- * TODO: flesh out the documentation
- */
-extern const time_point turn_zero;
 } // namespace calendar
 
 template<typename T>
@@ -504,6 +503,25 @@ class time_point
         // TODO: implement minutes_of_hour and so on and use it.
 };
 
+namespace calendar
+{
+
+/**
+ * A time point that is always before the current turn, even when the game has
+ * just started. This implies `before_time_starts < calendar::turn` is always
+ * true. It can be used to initialize `time_point` values that denote the last
+ * time a cache was update.
+ */
+constexpr time_point before_time_starts = time_point::from_turn( -1 );
+/**
+ * Represents time point 0.
+ * TODO: flesh out the documentation
+ */
+
+constexpr time_point turn_zero = time_point::from_turn( 0 );
+
+} // namespace calendar
+
 inline time_duration time_past_midnight( const time_point &p )
 {
     return ( p - calendar::turn_zero ) % 1_days;
@@ -560,15 +578,30 @@ bool is_day( const time_point &p );
 bool is_dusk( const time_point &p );
 /** Returns true if it's currently dawn - between sunrise and twilight_duration after sunrise. */
 bool is_dawn( const time_point &p );
-/** Returns the current seasonally-adjusted maximum daylight level */
-double current_daylight_level( const time_point &p );
 /** How much light is provided in full daylight */
 double default_daylight_level();
-/** Returns the current sunlight or moonlight level through the preceding functions.
- *  By default, returns sunlight level for vision, with moonlight providing a measurable amount
- *  of light.  with vision == false, returns sunlight for solar panel purposes, and moonlight
- *  provides 0 light */
-float sunlight( const time_point &p, bool vision = true );
+/** Returns the current sunlight.
+ *  Based entirely on astronomical circumstances; does not account for e.g.
+ *  weather.
+ *  For most situations you actually want to call the below function which also
+ *  includes moonlight. */
+float sun_light_at( const time_point &p );
+/** Returns the current sunlight plus moonlight level.
+ *  Based entirely on astronomical circumstances; does not account for e.g.
+ *  weather. */
+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 );
+
+/** 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.
+ *
+ * If lat_long not provided it defaults to Boston.
+ */
+cata::optional<rl_vec2d> sunlight_angle( const time_point &, lat_long );
+cata::optional<rl_vec2d> sunlight_angle( const time_point & );
 
 enum class weekdays : int {
     SUNDAY = 0,

src/character.cpp

@@ -12013,14 +12013,17 @@ bool Character::sees_with_infrared( const Creature &critter ) const
     }
 
     map &here = get_map();
+    // Use range based on default daylight, not actual current light, since
+    // we're seeing in infra red not via light.
+    int range = sight_range( default_daylight_level() );
+
     if( is_player() || critter.is_player() ) {
         // Players should not use map::sees
         // Likewise, players should not be "looked at" with map::sees, not to break symmetry
-        return here.pl_line_of_sight( critter.pos(),
-                                      sight_range( current_daylight_level( calendar::turn ) ) );
+        return here.pl_line_of_sight( critter.pos(), range );
     }
 
-    return here.sees( pos(), critter.pos(), sight_range( current_daylight_level( calendar::turn ) ) );
+    return here.sees( pos(), critter.pos(), range );
 }
 
 bool Character::is_visible_in_range( const Creature &critter, const int range ) const

src/editmap.cpp

@@ -743,9 +743,9 @@ void editmap::update_view_with_help( const std::string &txt, const std::string &
     const std::string u_see_msg = player_character.sees( target ) ? _( "yes" ) : _( "no" );
     mvwprintw( w_info, point( 1, off++ ), _( "dist: %d u_see: %s veh: %s scent: %d" ),
                rl_dist( player_character.pos(), target ), u_see_msg, veh_msg, get_scent().get( target ) );
-    mvwprintw( w_info, point( 1, off++ ), _( "sight_range: %d, daylight_sight_range: %d," ),
+    mvwprintw( w_info, point( 1, off++ ), _( "sight_range: %d, noon_sight_range: %d," ),
                player_character.sight_range( g->light_level( player_character.posz() ) ),
-               player_character.sight_range( current_daylight_level( calendar::turn ) ) );
+               player_character.sight_range( sun_moon_light_at_noon_near( calendar::turn ) ) );
     mvwprintw( w_info, point( 1, off++ ), _( "cache{transp:%.4f seen:%.4f cam:%.4f}" ),
                map_cache.transparency_cache[target.x][target.y],
                map_cache.seen_cache[target.x][target.y],

src/game.cpp

@@ -4037,7 +4037,7 @@ float game::natural_light_level( const int zlev ) const
 
     // Sunlight/moonlight related stuff
     if( !weather.lightning_active ) {
-        ret = sunlight( calendar::turn );
+        ret = sun_moon_light_at( calendar::turn );
     } else {
         // Recent lightning strike has lit the area
         ret = default_daylight_level();
@@ -7711,8 +7711,9 @@ void game::reset_item_list_state( const catacurses::window &window, int height,
 
 void game::list_items_monsters()
 {
-    std::vector<Creature *> mons = u.get_visible_creatures( current_daylight_level( calendar::turn ) );
-    // whole reality bubble
+    // Search whole reality bubble because each function internally verifies
+    // the visibilty of the items / monsters in question.
+    std::vector<Creature *> mons = u.get_visible_creatures( 60 );
     const std::vector<map_item_stack> items = find_nearby_items( 60 );
 
     if( mons.empty() && items.empty() ) {

src/line.cpp

@@ -276,6 +276,11 @@ units::angle atan2( const point &p )
     return units::atan2( p.y, p.x );
 }
 
+units::angle atan2( const rl_vec2d &p )
+{
+    return units::atan2( p.y, p.x );
+}
+
 // This more general version of this function gives correct values for larger values.
 unsigned make_xyz( const tripoint &p )
 {
@@ -562,6 +567,11 @@ std::vector<point> squares_in_direction( const point &p1, const point &p2 )
     return adjacent_squares;
 }
 
+rl_vec2d rl_vec3d::xy() const
+{
+    return rl_vec2d( x, y );
+}
+
 float rl_vec2d::magnitude() const
 {
     return std::sqrt( x * x + y * y );
@@ -661,13 +671,17 @@ rl_vec2d rl_vec2d::operator*( const float rhs ) const
     return ret;
 }
 
+rl_vec3d &rl_vec3d::operator*=( const float rhs )
+{
+    x *= rhs;
+    y *= rhs;
+    z *= rhs;
+    return *this;
+}
+
 rl_vec3d rl_vec3d::operator*( const float rhs ) const
 {
-    rl_vec3d ret;
-    ret.x = x * rhs;
-    ret.y = y * rhs;
-    ret.z = z * rhs;
-    return ret;
+    return rl_vec3d( *this ) *= rhs;
 }
 
 // subtract
@@ -731,13 +745,17 @@ rl_vec2d rl_vec2d::operator/( const float rhs ) const
     return ret;
 }
 
+rl_vec3d &rl_vec3d::operator/=( const float rhs )
+{
+    x /= rhs;
+    y /= rhs;
+    z /= rhs;
+    return *this;
+}
+
 rl_vec3d rl_vec3d::operator/( const float rhs ) const
 {
-    rl_vec3d ret;
-    ret.x = x / rhs;
-    ret.y = y / rhs;
-    ret.z = z / rhs;
-    return ret;
+    return rl_vec3d( *this ) /= rhs;
 }
 
 void calc_ray_end( units::angle angle, const int range, const tripoint &p, tripoint &out )

src/line.h

@@ -11,6 +11,8 @@
 #include "point.h"
 #include "units_fwd.h"
 
+struct rl_vec2d;
+
 extern bool trigdist;
 
 /**
@@ -231,6 +233,7 @@ int manhattan_dist( const point &loc1, const point &loc2 );
 
 // get angle of direction represented by point
 units::angle atan2( const point & );
+units::angle atan2( const rl_vec2d & );
 
 // Get the magnitude of the slope ranging from 0.0 to 1.0
 float get_normalized_angle( const point &start, const point &end );
@@ -282,8 +285,11 @@ struct rl_vec3d {
     float z;
 
     explicit rl_vec3d( float x = 0, float y = 0, float z = 0 ) : x( x ), y( y ), z( z ) {}
+    explicit rl_vec3d( const rl_vec2d &xy, float z = 0 ) : x( xy.x ), y( xy.y ), z( z ) {}
     explicit rl_vec3d( const tripoint &p ) : x( p.x ), y( p.y ), z( p.z ) {}
 
+    rl_vec2d xy() const;
+
     float magnitude() const;
     rl_vec3d normalized() const;
     rl_vec3d rotated( float angle ) const;
@@ -293,13 +299,15 @@ struct rl_vec3d {
     tripoint as_point() const;
 
     // scale.
-    rl_vec3d operator* ( float rhs ) const;
-    rl_vec3d operator/ ( float rhs ) const;
+    rl_vec3d &operator*=( float rhs );
+    rl_vec3d &operator/=( float rhs );
+    rl_vec3d operator*( float rhs ) const;
+    rl_vec3d operator/( float rhs ) const;
     // subtract
-    rl_vec3d operator- ( const rl_vec3d &rhs ) const;
+    rl_vec3d operator-( const rl_vec3d &rhs ) const;
     // unary negation
-    rl_vec3d operator- () const;
-    rl_vec3d operator+ ( const rl_vec3d &rhs ) const;
+    rl_vec3d operator-() const;
+    rl_vec3d operator+( const rl_vec3d &rhs ) const;
 };
 
 #endif // CATA_SRC_LINE_H

src/units.h

@@ -860,6 +860,16 @@ inline units::angle atan2( double y, double x )
     return from_radians( std::atan2( y, x ) );
 }
 
+inline units::angle asin( double x )
+{
+    return from_radians( std::asin( x ) );
+}
+
+inline units::angle acos( double x )
+{
+    return from_radians( std::acos( x ) );
+}
+
 static const std::vector<std::pair<std::string, energy>> energy_units = { {
         { "mJ", 1_mJ },
         { "J", 1_J },

src/units_utility.h

@@ -31,6 +31,13 @@ units::quantity<T, U> round_to_multiple_of( units::quantity<T, U> val, units::qu
     return multiple * of;
 }
 
+struct lat_long {
+    units::angle latitude;
+    units::angle longitude;
+};
+
+constexpr lat_long location_boston{ 42.36_degrees, -71.06_degrees };
+
 /**
  * Create a units label for a weight value.
  *

src/weather.cpp

@@ -130,7 +130,7 @@ void glare( const weather_type_id &w )
 
 int incident_sunlight( const weather_type_id &wtype, const time_point &t )
 {
-    return std::max<float>( 0.0f, sunlight( t, false ) + wtype->light_modifier );
+    return std::max<float>( 0.0f, sun_light_at( t ) + wtype->light_modifier );
 }
 
 static inline void proc_weather_sum( const weather_type_id &wtype, weather_sum &data,

tests/char_sight_test.cpp

@@ -56,8 +56,8 @@ TEST_CASE( "light and fine_detail_vision_mod", "[character][sight][light][vision
         // Build map cache including lightmap
         here.build_map_cache( 0, false );
         REQUIRE( g->is_in_sunlight( dummy.pos() ) );
-        // ambient_light_at is 100.0 in full sun (this fails if lightmap cache is not built)
-        REQUIRE( here.ambient_light_at( dummy.pos() ) == Approx( 100.0f ) );
+        // ambient_light_at is ~100.0 in full sun (this fails if lightmap cache is not built)
+        REQUIRE( here.ambient_light_at( dummy.pos() ) == Approx( 100.0f ).margin( 10 ) );
 
         // 1.0 is LIGHT_AMBIENT_LIT or brighter
         CHECK( dummy.fine_detail_vision_mod() == Approx( 1.0f ) );
@@ -263,7 +263,7 @@ TEST_CASE( "ursine vision", "[character][ursine][vision]" )
             here.build_map_cache( 0, false );
             light_here = here.ambient_light_at( dummy.pos() );
             REQUIRE( g->is_in_sunlight( dummy.pos() ) );
-            REQUIRE( light_here == Approx( 100.0f ) );
+            REQUIRE( light_here == Approx( 100.0f ).margin( 10 ) );
 
             THEN( "impaired sight, with 4 tiles of range" ) {
                 dummy.recalc_sight_limits();
@@ -281,7 +281,7 @@ TEST_CASE( "ursine vision", "[character][ursine][vision]" )
                     dummy.recalc_sight_limits();
                     CHECK_FALSE( dummy.sight_impaired() );
                     CHECK( dummy.unimpaired_range() == 60 );
-                    CHECK( dummy.sight_range( light_here ) == 87 );
+                    CHECK( dummy.sight_range( light_here ) == 88 );
                 }
             }
         }