Introduce observation events for Moon (#86)

Like the Sun, the Moon is a candidate for events defined in `Astronoby::Events::ObservationEvents`, such as rising, transit and setting times, or rising and setting azimuths, or transit altitude.

This adds `Astronoby::Moon#observation_events` method, with an `observer` as parameter.

This change also refactors how observation events are computed, introducing a loop to increase precision, as explained in _Astronomical Algorithms_.

From a few tests, the error margin with the IMCCE is around 2 minutes maximum.

```rb
time = Time.utc(2024, 5, 31)
observer = Astronoby::Observer.new(
  latitude: Astronoby::Angle.zero,
  longitude: Astronoby::Angle.zero
)
moon = Astronoby::Moon.new(time: time)

observation_events = moon.observation_events(observer: observer)

observation_events.rising_time
# => 2024-05-31 00:29:50 UTC

observation_events.transit_time
# => 2024-05-31 06:41:22 UTC

observation_events.setting_time
# => 2024-05-31 12:52:48 UTC

observation_events.rising_azimuth.str(:dms)
# => "+98° 47′ 20.2285″"

observation_events.transit_altitude.str(:dms)
# => "+83° 2′ 9.8145″"

observation_events.setting_azimuth.str(:dms)
# => "+264° 35′ 11.0134″"
```

lib/astronoby.rb

@@ -18,6 +18,7 @@
 require "astronoby/errors"
 require "astronoby/events/moon_phases"
 require "astronoby/events/observation_events"
+require "astronoby/events/rise_transit_set_iteration"
 require "astronoby/events/twilight_events"
 require "astronoby/geocentric_parallax"
 require "astronoby/mean_obliquity"

lib/astronoby/bodies/moon.rb

@@ -2,6 +2,7 @@
 
 module Astronoby
   class Moon
+    SEMIDIAMETER_VARIATION = 0.7275
     MEAN_GEOCENTRIC_DISTANCE = Astronoby::Distance.from_meters(385_000_560)
 
     # Source:
@@ -166,6 +167,44 @@ def illuminated_fraction
       @illuminated_fraction ||= (1 + phase_angle.cos) / 2
     end
 
+    # @param observer [Astronoby::Observer] Observer of the event
+    # @return [Astronoby::Events::ObservationEvents] Moon's observation events
+    def observation_events(observer:)
+      today = @time.to_date
+      leap_seconds = Util::Time.terrestrial_universal_time_delta(today)
+      yesterday = today.prev_day
+      yesterday_midnight_terrestrial_time =
+        Time.utc(yesterday.year, yesterday.month, yesterday.day) - leap_seconds
+      today_midnight_terrestrial_time =
+        Time.utc(today.year, today.month, today.day) - leap_seconds
+      tomorrow = today.next_day
+      tomorrow_midnight_terrestrial_time =
+        Time.utc(tomorrow.year, tomorrow.month, tomorrow.day) - leap_seconds
+
+      coordinates_of_the_previous_day = self.class
+        .new(time: yesterday_midnight_terrestrial_time)
+        .apparent_equatorial_coordinates
+      coordinates_of_the_day = self.class
+        .new(time: today_midnight_terrestrial_time)
+        .apparent_equatorial_coordinates
+      coordinates_of_the_next_day = self.class
+        .new(time: tomorrow_midnight_terrestrial_time)
+        .apparent_equatorial_coordinates
+      additional_altitude = -Angle.from_degrees(
+        SEMIDIAMETER_VARIATION *
+          GeocentricParallax.angle(distance: distance).degrees
+      )
+
+      Events::ObservationEvents.new(
+        observer: observer,
+        date: today,
+        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: additional_altitude
+      )
+    end
+
     private
 
     def terrestrial_time

lib/astronoby/events/observation_events.rb

@@ -6,6 +6,7 @@ class ObservationEvents
       STANDARD_ALTITUDE = Angle.from_dms(0, -34, 0)
       RISING_SETTING_HOUR_ANGLE_RATIO_RANGE = (-1..1)
       EARTH_SIDEREAL_ROTATION_RATE = 360.98564736629
+      ITERATION_PRECISION = 0.0001
 
       attr_reader :rising_time,
         :rising_azimuth,
@@ -51,46 +52,80 @@ def initialize(
       private
 
       def compute
-        @transit_time = Util::Time.decimal_hour_to_time(@date, initial_transit)
+        @initial_transit = initial_transit
+        @transit_time = Util::Time.decimal_hour_to_time(@date, @initial_transit)
         @transit_altitude = local_horizontal_altitude_transit
 
         return if h0.nil?
 
-        delta_m_rising = (local_horizontal_altitude_rising - shift).degrees./(
-          Constants::DEGREES_PER_CIRCLE *
-            declination_rising.cos *
-            @observer.latitude.cos *
-            local_hour_angle_rising.sin
+        initial_rising = rationalize_decimal_time(
+          @initial_transit - h0.degrees / Constants::DEGREES_PER_CIRCLE
         )
-        delta_m_transit = -local_hour_angle_transit.degrees / Constants::DEGREES_PER_CIRCLE
-        delta_m_setting = (local_horizontal_altitude_setting - shift).degrees./(
-          Constants::DEGREES_PER_CIRCLE *
-            declination_setting.cos *
-            @observer.latitude.cos *
-            local_hour_angle_setting.sin
+
+        initial_setting = rationalize_decimal_time(
+          @initial_transit + h0.degrees / Constants::DEGREES_PER_CIRCLE
         )
 
-        corrected_rising = rationalize_decimal_hours(
-          Constants::HOURS_PER_DAY * (initial_rising + delta_m_rising)
+        @final_rising, @final_transit, @final_setting =
+          iterate(initial_rising, @initial_transit, initial_setting)
+
+        rationalized_corrected_rising = rationalize_decimal_hours(
+          Constants::HOURS_PER_DAY * @final_rising
         )
-        corrected_transit = rationalize_decimal_hours(
-          Constants::HOURS_PER_DAY * (initial_transit + delta_m_transit)
+        rationalized_corrected_transit = rationalize_decimal_hours(
+          Constants::HOURS_PER_DAY * @final_transit
         )
-        corrected_setting = rationalize_decimal_hours(
-          Constants::HOURS_PER_DAY * (initial_setting + delta_m_setting)
+        rationalized_corrected_setting = rationalize_decimal_hours(
+          Constants::HOURS_PER_DAY * @final_setting
         )
 
-        @rising_time = Util::Time.decimal_hour_to_time(@date, corrected_rising)
+        @rising_time = Util::Time.decimal_hour_to_time(@date, rationalized_corrected_rising)
         @rising_azimuth = local_horizontal_azimuth_rising
-        @transit_time = Util::Time.decimal_hour_to_time(@date, corrected_transit)
-        @setting_time = Util::Time.decimal_hour_to_time(@date, corrected_setting)
+        @transit_time = Util::Time.decimal_hour_to_time(@date, rationalized_corrected_transit)
+        @transit_altitude = local_horizontal_altitude_transit
+        @setting_time = Util::Time.decimal_hour_to_time(@date, rationalized_corrected_setting)
         @setting_azimuth = local_horizontal_azimuth_setting
       end
 
+      def iterate(initial_rising, initial_transit, initial_setting)
+        delta = 1
+        corrected_rising = initial_rising
+        corrected_transit = initial_transit
+        corrected_setting = initial_setting
+        until delta < ITERATION_PRECISION
+          iterate = RiseTransitSetIteration.new(
+            observer: @observer,
+            date: @date,
+            coordinates_of_the_next_day: @coordinates_of_the_next_day,
+            coordinates_of_the_day: @coordinates_of_the_day,
+            coordinates_of_the_previous_day: @coordinates_of_the_previous_day,
+            shift: shift,
+            initial_rising: corrected_rising,
+            initial_transit: corrected_transit,
+            initial_setting: corrected_setting
+          ).iterate
+          delta = iterate.sum
+          corrected_rising = rationalize_decimal_time corrected_rising + iterate[0]
+          corrected_transit = rationalize_decimal_time corrected_transit + iterate[1]
+          corrected_setting = rationalize_decimal_time corrected_setting + iterate[2]
+        end
+        [corrected_rising, corrected_transit, corrected_setting]
+      end
+
       def observer_longitude
         # Longitude must be treated positively westwards from the meridian of
         # Greenwich, and negatively to the east
-        @observer_longitude ||= -@observer.longitude
+        -@observer.longitude
+      end
+
+      def initial_transit
+        rationalize_decimal_time(
+          (
+            @coordinates_of_the_day.right_ascension.degrees +
+              observer_longitude.degrees -
+              apparent_gst_at_midnight.degrees
+          ) / Constants::DEGREES_PER_CIRCLE
+        )
       end
 
       def h0
@@ -106,115 +141,46 @@ def h0
       end
 
       def apparent_gst_at_midnight
-        @apparent_gst_at_midnight ||= Angle.from_hours(
+        Angle.from_hours(
           GreenwichSiderealTime.from_utc(
             Time.utc(@date.year, @date.month, @date.day)
           ).time
         )
       end
 
-      def initial_transit
-        @initial_transit ||= rationalize_decimal_time(
-          (
-            @coordinates_of_the_day.right_ascension.degrees +
-              observer_longitude.degrees -
-              apparent_gst_at_midnight.degrees
-          ) / Constants::DEGREES_PER_CIRCLE
-        )
-      end
-
-      def initial_rising
-        @initial_rising ||=
-          rationalize_decimal_time(
-            initial_transit - h0.degrees / Constants::DEGREES_PER_CIRCLE
-          )
-      end
-
-      def initial_setting
-        @initial_setting ||=
-          rationalize_decimal_time(
-            initial_transit + h0.degrees / Constants::DEGREES_PER_CIRCLE
-          )
-      end
-
-      def gst_rising
-        @gst_rising ||= Angle.from_degrees(
-          apparent_gst_at_midnight.degrees +
-            EARTH_SIDEREAL_ROTATION_RATE * initial_rising
-        )
-      end
-
       def gst_transit
-        @gst_transit ||= Angle.from_degrees(
-          apparent_gst_at_midnight.degrees +
-            EARTH_SIDEREAL_ROTATION_RATE * initial_transit
-        )
-      end
-
-      def gst_setting
-        @gst_setting ||= Angle.from_degrees(
+        Angle.from_degrees(
           apparent_gst_at_midnight.degrees +
-            EARTH_SIDEREAL_ROTATION_RATE * initial_setting
+            EARTH_SIDEREAL_ROTATION_RATE * (@final_transit || @initial_transit)
         )
       end
 
       def leap_day_portion
-        @leap_day_portion ||= begin
-          leap_seconds = Util::Time.terrestrial_universal_time_delta(@date)
-          leap_seconds / Constants::SECONDS_PER_DAY
-        end
-      end
-
-      def local_hour_angle_rising
-        @local_hour_angle_rising ||=
-          gst_rising - observer_longitude - right_ascension_rising
+        leap_seconds = Util::Time.terrestrial_universal_time_delta(@date)
+        leap_seconds / Constants::SECONDS_PER_DAY
       end
 
       def local_hour_angle_transit
-        @local_hour_angle_transit ||=
-          gst_transit - observer_longitude - right_ascension_transit
-      end
-
-      def local_hour_angle_setting
-        @local_hour_angle_setting ||=
-          gst_setting - observer_longitude - right_ascension_setting
-      end
-
-      def local_horizontal_altitude_rising
-        @local_horizontal_altitude_rising ||= Angle.asin(
-          @observer.latitude.sin * declination_rising.sin +
-            @observer.latitude.cos * declination_rising.cos * local_hour_angle_rising.cos
-        )
+        gst_transit - observer_longitude - right_ascension_transit
       end
 
       def local_horizontal_azimuth_rising
-        @local_horizontal_azimuth_rising ||= begin
-          shift = -@standard_altitude
-          term1 = declination_rising.sin + shift.sin * @observer.latitude.cos
-          term2 = shift.cos * @observer.latitude.cos
-          angle = term1 / term2
-          Angle.acos(angle)
-        end
+        term1 = declination_rising.sin + (-shift).sin * @observer.latitude.cos
+        term2 = (-shift).cos * @observer.latitude.cos
+        angle = term1 / term2
+        Angle.acos(angle)
       end
 
       def local_horizontal_altitude_transit
-        @local_horizontal_altitude_transit ||= Angle.asin(
+        Angle.asin(
           @observer.latitude.sin * declination_transit.sin +
             @observer.latitude.cos * declination_transit.cos * local_hour_angle_transit.cos
         )
       end
 
-      def local_horizontal_altitude_setting
-        @local_horizontal_altitude_setting ||= Angle.asin(
-          @observer.latitude.sin * declination_setting.sin +
-            @observer.latitude.cos * declination_setting.cos * local_hour_angle_setting.cos
-        )
-      end
-
       def local_horizontal_azimuth_setting
-        shift = -@standard_altitude
-        term1 = declination_setting.sin + shift.sin * @observer.latitude.cos
-        term2 = shift.cos * @observer.latitude.cos
+        term1 = declination_setting.sin + (-shift).sin * @observer.latitude.cos
+        term2 = (-shift).cos * @observer.latitude.cos
         angle = term1 / term2
         Angle.from_degrees(
           Constants::DEGREES_PER_CIRCLE - Angle.acos(angle).degrees
@@ -233,19 +199,6 @@ def rationalize_decimal_hours(decimal_hours)
         decimal_hours
       end
 
-      def right_ascension_rising
-        Angle.from_degrees(
-          Util::Maths.interpolate(
-            [
-              @coordinates_of_the_previous_day.right_ascension.degrees,
-              @coordinates_of_the_day.right_ascension.degrees,
-              @coordinates_of_the_next_day.right_ascension.degrees
-            ],
-            initial_rising + leap_day_portion
-          )
-        )
-      end
-
       def right_ascension_transit
         Angle.from_degrees(
           Util::Maths.interpolate(
@@ -254,20 +207,7 @@ def right_ascension_transit
               @coordinates_of_the_day.right_ascension.degrees,
               @coordinates_of_the_next_day.right_ascension.degrees
             ],
-            initial_transit + leap_day_portion
-          )
-        )
-      end
-
-      def right_ascension_setting
-        Angle.from_degrees(
-          Util::Maths.interpolate(
-            [
-              @coordinates_of_the_previous_day.right_ascension.degrees,
-              @coordinates_of_the_day.right_ascension.degrees,
-              @coordinates_of_the_next_day.right_ascension.degrees
-            ],
-            initial_setting + leap_day_portion
+            (@final_transit || @initial_transit) + leap_day_portion
           )
         )
       end
@@ -280,7 +220,7 @@ def declination_rising
               @coordinates_of_the_day.declination.degrees,
               @coordinates_of_the_next_day.declination.degrees
             ],
-            initial_rising + leap_day_portion
+            @final_rising + leap_day_portion
           )
         )
       end
@@ -293,7 +233,7 @@ def declination_transit
               @coordinates_of_the_day.declination.degrees,
               @coordinates_of_the_next_day.declination.degrees
             ],
-            initial_transit + leap_day_portion
+            (@final_transit || @initial_transit) + leap_day_portion
           )
         )
       end
@@ -306,7 +246,7 @@ def declination_setting
               @coordinates_of_the_day.declination.degrees,
               @coordinates_of_the_next_day.declination.degrees
             ],
-            initial_setting + leap_day_portion
+            @final_setting + leap_day_portion
           )
         )
       end