More accurate rising, transit and setting times (#50)

This change is a complete refactoring of rising and setting times and azimuth.
It is based on Astronomical Algorithms by Jean Meeus with a more complex but more accurate method than the ones used so far.
It also introduces transit time and altitude.

The overall precision of these events have been increased, especially for the Sun. Most of the times are accurate with the IMCCE within 2 minutes precision, with results sometimes accurate to a few seconds.

More work may be needed in the future to increase the precision for distant stars, that should be in theory simpler to deal with but are less accurate for some reason (angles are off by a dozen minutes of arc).

The following methods have been added to `Sun`:
- `#rise_transit_set_times`
- `#rise_set_azimuths`
- `#transit_time`
- `#transit_altitude`
They all require an `observer` key argument. Memoization on the observer has been implemented so that each of these methods calls don't compute the values all over again.

This change includes a breaking change:
`Astronoby::Body` is dropped. It will probably be replaced by something like `DistantStar`, `Moon` and `Planet` in the future, in the same way we have `Sun`.

```rb
date = Date.new(2015, 2, 5)
epoch = Astronoby::Epoch.from_time(date)
observer = Astronoby::Observer.new(
  latitude: Astronoby::Angle.from_degrees(38),
  longitude: Astronoby::Angle.from_degrees(-78)
)
sun = Astronoby::Sun.new(epoch: epoch)

sun.rising_time(observer: observer)
# => 2015-02-05 12:12:59 UTC

sun.rising_azimuth(observer: observer).str(:dms)
# => "+109° 46′ 43.1427″"

sun.transit_time(observer: observer)
# => 2015-02-05 17:25:59 UTC

sun.transit_altitude(observer: observer).str(:dms)
# => "+36° 8′ 15.7669″"

sun.setting_time(observer: observer)
# => 2015-02-05 22:39:27 UTC

sun.setting_azimuth(observer: observer).str(:dms)
# => "+250° 23′ 33.6177″"

# Also available in more compact forms:
sun.rise_transit_set_times(observer: observer)
# => [..., ..., ...]

sun.rise_set_azimuths(observer: observer)
# => [..., ...]
```

README.md

@@ -102,16 +102,29 @@ observer = Astronoby::Observer.new(
 sun = Astronoby::Sun.new(epoch: epoch)
 
 sun.rising_time(observer: observer)
-# => 2015-02-05 12:13:26 UTC
+# => 2015-02-05 12:12:59 UTC
 
 sun.rising_azimuth(observer: observer).str(:dms)
-# => "+109° 41′ 22.2585″"
+# => "+109° 46′ 43.1427″"
+
+sun.transit_time(observer: observer)
+# => 2015-02-05 17:25:59 UTC
+
+sun.transit_altitude(observer: observer).str(:dms)
+# => "+36° 8′ 15.7669″"
 
 sun.setting_time(observer: observer)
-# => 2015-02-05 22:35:12 UTC
+# => 2015-02-05 22:39:27 UTC
 
 sun.setting_azimuth(observer: observer).str(:dms)
-# => "+250° 18′ 37.7414″"
+# => "+250° 23′ 33.6177″"
+
+# Also available in more compact forms:
+sun.rise_transit_set_times(observer: observer)
+# => [..., ..., ...]
+
+sun.rise_set_azimuths(observer: observer)
+# => [..., ...]
 ```
 
 ### Solstice and Equinox times

lib/astronoby.rb

@@ -4,7 +4,6 @@
 require "astronoby/angles/dms"
 require "astronoby/angles/hms"
 require "astronoby/epoch"
-require "astronoby/body"
 require "astronoby/bodies/sun"
 require "astronoby/coordinates/ecliptic"
 require "astronoby/coordinates/equatorial"
@@ -18,6 +17,7 @@
 require "astronoby/observer"
 require "astronoby/precession"
 require "astronoby/refraction"
+require "astronoby/rise_transit_set"
 require "astronoby/time/greenwich_sidereal_time"
 require "astronoby/time/local_sidereal_time"
 require "astronoby/util/astrodynamics"

lib/astronoby/bodies/sun.rb

@@ -75,58 +75,51 @@ def horizontal_coordinates(latitude:, longitude:)
         .to_horizontal(time: time, latitude: latitude, longitude: longitude)
     end
 
+    # @param observer [Astronoby::Observer] Observer of the events
+    # @return [Array<Time, Time, Time>, nil] Sunrise, transit and sunset times
+    def rise_transit_set_times(observer:)
+      rise_transit_set(observer).times
+    end
+
+    # @param observer [Astronoby::Observer] Observer of the events
+    # @return [Array<Astronoby::Angle, Astronoby::Angle>] Azimuths of sunrise
+    #   and sunset
+    def rise_set_azimuths(observer:)
+      rise_transit_set(observer).azimuths
+    end
+
     # @param observer [Astronoby::Observer] Observer of the event
     # @return [Time] Time of sunrise
     def rising_time(observer:)
-      event_date = Epoch.to_utc(@epoch).to_date
-      lst1 = event_local_sidereal_time_for_date(event_date, observer, :rising)
-      next_day = event_date.next_day(1)
-      lst2 = event_local_sidereal_time_for_date(next_day, observer, :rising)
-      time = (INTERPOLATION_FACTOR * lst1) / (INTERPOLATION_FACTOR + lst1 - lst2)
-
-      LocalSiderealTime.new(
-        date: event_date,
-        time: time,
-        longitude: observer.longitude
-      ).to_gst.to_utc
+      rise_transit_set(observer).times[0]
     end
 
     # @param observer [Astronoby::Observer] Observer of the event
     # @return [Astronoby::Angle, nil] Azimuth of sunrise
     def rising_azimuth(observer:)
-      equatorial_coordinates = apparent_ecliptic_coordinates
-        .to_apparent_equatorial(epoch: @epoch)
-      Body.new(equatorial_coordinates).rising_azimuth(
-        latitude: observer.latitude,
-        vertical_shift: vertical_shift
-      )
+      rise_transit_set(observer).azimuths[0]
     end
 
     # @param observer [Astronoby::Observer] Observer of the event
     # @return [Time] Time of sunset
     def setting_time(observer:)
-      event_date = Epoch.to_utc(@epoch).to_date
-      lst1 = event_local_sidereal_time_for_date(event_date, observer, :setting)
-      next_day = event_date.next_day(1)
-      lst2 = event_local_sidereal_time_for_date(next_day, observer, :setting)
-      time = (INTERPOLATION_FACTOR * lst1) / (INTERPOLATION_FACTOR + lst1 - lst2)
-
-      LocalSiderealTime.new(
-        date: event_date,
-        time: time,
-        longitude: observer.longitude
-      ).to_gst.to_utc
+      rise_transit_set(observer).times[2]
     end
 
     # @param observer [Astronoby::Observer] Observer of the event
     # @return [Astronoby::Angle, nil] Azimuth of sunset
     def setting_azimuth(observer:)
-      equatorial_coordinates = apparent_ecliptic_coordinates
-        .to_apparent_equatorial(epoch: @epoch)
-      Body.new(equatorial_coordinates).setting_azimuth(
-        latitude: observer.latitude,
-        vertical_shift: vertical_shift
-      )
+      rise_transit_set(observer).azimuths[1]
+    end
+
+    def transit_time(observer:)
+      rise_transit_set(observer).times[1]
+    end
+
+    # @param observer [Astronoby::Observer] Observer of the event
+    # @return [Astronoby::Angle] Altitude at transit
+    def transit_altitude(observer:)
+      rise_transit_set(observer).transit_altitude
     end
 
     # @return [Numeric] Earth-Sun distance in meters
@@ -206,43 +199,39 @@ def distance_angular_size_factor
       term1 / term2
     end
 
-    def event_local_sidereal_time_for_date(date, observer, event)
-      midnight_utc = Time.utc(date.year, date.month, date.day)
-      epoch = Epoch.from_time(midnight_utc)
-      sun_at_midnight = self.class.new(epoch: epoch)
-      shift = Body::DEFAULT_REFRACTION_VERTICAL_SHIFT +
-        GeocentricParallax.angle(distance: sun_at_midnight.earth_distance) +
-        Angle.from_degrees(sun_at_midnight.angular_size.degrees / 2)
-      ecliptic_coordinates = sun_at_midnight.apparent_ecliptic_coordinates
-      equatorial_coordinates = ecliptic_coordinates
-        .to_apparent_equatorial(epoch: epoch)
-
-      event_time = if event == :rising
-        Body.new(equatorial_coordinates).rising_time(
-          latitude: observer.latitude,
-          longitude: observer.longitude,
-          date: midnight_utc.to_date,
-          vertical_shift: shift
-        )
-      else
-        Body.new(equatorial_coordinates).setting_time(
-          latitude: observer.latitude,
-          longitude: observer.longitude,
-          date: midnight_utc.to_date,
-          vertical_shift: shift
+    def current_date
+      Epoch.to_utc(@epoch).to_date
+    end
+
+    def rise_transit_set(observer)
+      @rise_transit_set ||= {}
+      return @rise_transit_set[observer] if @rise_transit_set.key?(observer)
+
+      @rise_transit_set[observer] = begin
+        date = Epoch.to_utc(@epoch).to_date
+        yesterday_epoch = Epoch.from_time(date.prev_day)
+        tomorrow_epoch = Epoch.from_time(date.next_day)
+
+        coordinates_of_the_previous_day = self.class
+          .new(epoch: yesterday_epoch)
+          .apparent_ecliptic_coordinates
+          .to_apparent_equatorial(epoch: yesterday_epoch)
+        coordinates_of_the_day =
+          apparent_ecliptic_coordinates.to_apparent_equatorial(epoch: @epoch)
+        coordinates_of_the_next_day = self.class
+          .new(epoch: tomorrow_epoch)
+          .apparent_ecliptic_coordinates
+          .to_apparent_equatorial(epoch: tomorrow_epoch)
+
+        RiseTransitSet.new(
+          observer: observer,
+          date: date,
+          coordinates_of_the_previous_day: coordinates_of_the_previous_day,
+          coordinates_of_the_day: coordinates_of_the_day,
+          coordinates_of_the_next_day: coordinates_of_the_next_day,
+          additional_altitude: Angle.from_degrees(angular_size.degrees / 2)
         )
       end
-
-      GreenwichSiderealTime
-        .from_utc(event_time.utc)
-        .to_lst(longitude: observer.longitude)
-        .time
-    end
-
-    def vertical_shift
-      Astronoby::Body::DEFAULT_REFRACTION_VERTICAL_SHIFT +
-        Astronoby::GeocentricParallax.angle(distance: earth_distance) +
-        Astronoby::Angle.from_degrees(angular_size.degrees / 2)
     end
   end
 end