Handle Sun's sunrise and sunset azimuth (#39)

Knowing when the Sun rises or sets is interesting, but it can also be
interesting to know _where_ on the horizon these events are going to
happen.

The location of the Sun on the horizon's axis is the azimuth.

This change adds support of the computing the azimuth angle of the Sun's
at sunrise and sunset.

`Sun#rising_azimuth` and `Sun#setting_azimuth` can return `nil` as the
Sun sometimes never rises or never sets at some places on Earth, like
locations close to the poles.

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

sun.setting_azimuth(observer: observer)&.str(:dms)
```

UPGRADING.md

@@ -58,12 +58,20 @@ Returns the apparent Sun's longitude (`Angle`) at its perigee.
 
 ### `Sun#rising_time` method added (#35)
 
-Returns the UTC `Time` of the sunrise.
+Returns the UTC `Time` of the sunrise.`
+
+### `Sun#rising_azimuth` method added (#39)
+
+Returns the Sun's azimuth (`Angle`) at sunrise.
 
 ### `Sun#setting_time` method added (#35)
 
 Returns the UTC `Time` of the sunset.
 
+### `Sun#setting_azimuth` method added (#39)
+
+Returns the Sun's azimuth (`Angle`) at sunset.
+
 ### Added comparison methods to `Angle` (#21)
 
 With the inclusion of `Comparable`, comparison methods such as `#==`, `#<`,

lib/astronoby/bodies/sun.rb

@@ -56,6 +56,16 @@ def rising_time(observer:)
       ).to_gst.to_utc
     end
 
+    # @param observer [Astronoby::Observer] Observer of the event
+    # @return [Astronoby::Angle, nil] Azimuth of sunrise
+    def rising_azimuth(observer:)
+      equatorial_coordinates = ecliptic_coordinates.to_equatorial(epoch: @epoch)
+      Body.new(equatorial_coordinates).rising_azimuth(
+        latitude: observer.latitude,
+        vertical_shift: vertical_shift
+      )
+    end
+
     # @param observer [Astronoby::Observer] Observer of the event
     # @return [Time] Time of sunset
     def setting_time(observer:)
@@ -72,6 +82,16 @@ def setting_time(observer:)
       ).to_gst.to_utc
     end
 
+    # @param observer [Astronoby::Observer] Observer of the event
+    # @return [Astronoby::Angle, nil] Azimuth of sunset
+    def setting_azimuth(observer:)
+      equatorial_coordinates = ecliptic_coordinates.to_equatorial(epoch: @epoch)
+      Body.new(equatorial_coordinates).setting_azimuth(
+        latitude: observer.latitude,
+        vertical_shift: vertical_shift
+      )
+    end
+
     # @return [Numeric] Earth-Sun distance in meters
     def earth_distance
       SEMI_MAJOR_AXIS_IN_METERS / distance_angular_size_factor
@@ -172,5 +192,11 @@ def event_local_sidereal_time_for_date(date, observer, event)
         .to_lst(longitude: observer.longitude)
         .time
     end
+
+    def vertical_shift
+      Astronoby::Body::DEFAULT_REFRACTION_VERTICAL_SHIFT +
+        Astronoby::GeocentricParallax.angle(distance: earth_distance) +
+        Astronoby::Angle.as_degrees(angular_size.degrees / 2)
+    end
   end
 end

lib/astronoby/body.rb

@@ -14,78 +14,106 @@ def initialize(equatorial_coordinates)
     #  Authors: Peter Duffett-Smith and Jonathan Zwart
     #  Edition: Cambridge University Press
     #  Chapter: 33 - Rising and setting
+
+    # @param latitude [Astronoby::Angle] Latitude of the observer
+    # @param longitude [Astronoby::Angle] Longitude of the observer
+    # @param date [Date] Date of the event
+    # @param apparent [Boolean] Compute apparent or true data
+    # @param vertical_shift [Astronoby::Angle] Vertical shift correction angle
+    # @return [Time, nil] Sunrise time
     def rising_time(
       latitude:,
       longitude:,
       date:,
       apparent: true,
       vertical_shift: nil
     )
-      ratio = ratio(latitude, apparent, vertical_shift)
-      return nil unless RISING_SETTING_HOUR_ANGLE_RATIO_RANGE.cover?(ratio)
+      time_ratio = time_ratio(latitude, apparent, vertical_shift)
+      return nil unless RISING_SETTING_HOUR_ANGLE_RATIO_RANGE.cover?(time_ratio)
 
-      hour_angle = Angle.acos(ratio)
+      hour_angle = Angle.acos(time_ratio)
       local_sidereal_time = LocalSiderealTime.new(
         date: date,
-        time: @equatorial_coordinates.right_ascension.hours - hour_angle.hours,
+        time: right_ascension.hours - hour_angle.hours,
         longitude: longitude
       )
 
       local_sidereal_time.to_gst.to_utc
     end
 
     # Source:
-    #  Title: Celestial Calculations
-    #  Author: J. L. Lawrence
-    #  Edition: MIT Press
-    #  Chapter: 5 - Stars in the Nighttime Sky
-    def rising_azimuth(latitude:)
-      ar = azimuth_component(latitude)
-      return nil if ar >= 1
-
-      Angle.acos(ar)
+    #  Title: Practical Astronomy with your Calculator or Spreadsheet
+    #  Authors: Peter Duffett-Smith and Jonathan Zwart
+    #  Edition: Cambridge University Press
+    #  Chapter: 33 - Rising and setting
+
+    # @param latitude [Astronoby::Angle] Latitude of the observer
+    # @param apparent [Boolean] Compute apparent or true data
+    # @param vertical_shift [Astronoby::Angle] Vertical shift correction angle
+    # @return [Astronoby::Angle, nil] Sunrise azimuth
+    def rising_azimuth(latitude:, apparent: true, vertical_shift: nil)
+      time_ratio = time_ratio(latitude, apparent, vertical_shift)
+      return nil unless RISING_SETTING_HOUR_ANGLE_RATIO_RANGE.cover?(time_ratio)
+
+      azimuth_ratio = azimuth_ratio(latitude, apparent, vertical_shift)
+
+      Angle.acos(azimuth_ratio)
     end
 
     # Source:
     #  Title: Practical Astronomy with your Calculator or Spreadsheet
     #  Authors: Peter Duffett-Smith and Jonathan Zwart
     #  Edition: Cambridge University Press
     #  Chapter: 33 - Rising and setting
+
+    # @param latitude [Astronoby::Angle] Latitude of the observer
+    # @param longitude [Astronoby::Angle] Longitude of the observer
+    # @param date [Date] Date of the event
+    # @param apparent [Boolean] Compute apparent or true data
+    # @param vertical_shift [Astronoby::Angle] Vertical shift correction angle
+    # @return [Time, nil] Sunset time
     def setting_time(
       latitude:,
       longitude:,
       date:,
       apparent: true,
       vertical_shift: nil
     )
-      ratio = ratio(latitude, apparent, vertical_shift)
-      return nil unless RISING_SETTING_HOUR_ANGLE_RATIO_RANGE.cover?(ratio)
+      time_ratio = time_ratio(latitude, apparent, vertical_shift)
+      return nil unless RISING_SETTING_HOUR_ANGLE_RATIO_RANGE.cover?(time_ratio)
 
-      hour_angle = Angle.acos(ratio)
+      hour_angle = Angle.acos(time_ratio)
       local_sidereal_time = LocalSiderealTime.new(
         date: date,
-        time: @equatorial_coordinates.right_ascension.hours + hour_angle.hours,
+        time: right_ascension.hours + hour_angle.hours,
         longitude: longitude
       )
 
       local_sidereal_time.to_gst.to_utc
     end
 
     # Source:
-    #  Title: Celestial Calculations
-    #  Author: J. L. Lawrence
-    #  Edition: MIT Press
-    #  Chapter: 5 - Stars in the Nighttime Sky
-    def setting_azimuth(latitude:)
-      rising_az = rising_azimuth(latitude: latitude)
-      return nil if rising_az.nil?
-
-      Angle.as_degrees(360 - rising_az.degrees)
+    #  Title: Practical Astronomy with your Calculator or Spreadsheet
+    #  Authors: Peter Duffett-Smith and Jonathan Zwart
+    #  Edition: Cambridge University Press
+    #  Chapter: 33 - Rising and setting
+
+    # @param latitude [Astronoby::Angle] Latitude of the observer
+    # @param apparent [Boolean] Compute apparent or true data
+    # @param vertical_shift [Astronoby::Angle] Vertical shift correction angle
+    # @return [Astronoby::Angle, nil] Sunset azimuth
+    def setting_azimuth(latitude:, apparent: true, vertical_shift: nil)
+      time_ratio = time_ratio(latitude, apparent, vertical_shift)
+      return nil unless RISING_SETTING_HOUR_ANGLE_RATIO_RANGE.cover?(time_ratio)
+
+      azimuth_ratio = azimuth_ratio(latitude, apparent, vertical_shift)
+
+      Angle.as_degrees(360 - Angle.acos(azimuth_ratio).degrees)
     end
 
     private
 
-    def ratio(latitude, apparent, vertical_shift)
+    def time_ratio(latitude, apparent, vertical_shift)
       shift = if vertical_shift
         vertical_shift
       elsif apparent
@@ -94,13 +122,34 @@ def ratio(latitude, apparent, vertical_shift)
         Angle.zero
       end
 
-      -(shift.sin + latitude.sin * @equatorial_coordinates.declination.sin)./(
-        latitude.cos * @equatorial_coordinates.declination.cos
-      )
+      term1 = shift.sin + latitude.sin * declination.sin
+      term2 = latitude.cos * declination.cos
+
+      -term1 / term2
+    end
+
+    def azimuth_ratio(latitude, apparent, vertical_shift)
+      shift = if vertical_shift
+        vertical_shift
+      elsif apparent
+        DEFAULT_REFRACTION_VERTICAL_SHIFT
+      else
+        Angle.zero
+      end
+
+      (declination.sin + shift.sin * latitude.cos) / (shift.cos * latitude.cos)
     end
 
     def azimuth_component(latitude)
-      @equatorial_coordinates.declination.sin / latitude.cos
+      declination.sin / latitude.cos
+    end
+
+    def right_ascension
+      @equatorial_coordinates.right_ascension
+    end
+
+    def declination
+      @equatorial_coordinates.declination
     end
   end
 end

spec/astronoby/bodies/sun_spec.rb

@@ -378,6 +378,67 @@
     end
   end
 
+  describe "#rising_azimuth" do
+    it "returns an Angle" do
+      date = Date.new
+      epoch = Astronoby::Epoch.from_time(date)
+      observer = Astronoby::Observer.new(
+        latitude: Astronoby::Angle.zero,
+        longitude: Astronoby::Angle.zero
+      )
+      sun = described_class.new(epoch: epoch)
+
+      rising_azimuth = sun.rising_azimuth(observer: observer)
+
+      expect(rising_azimuth).to be_a(Astronoby::Angle)
+    end
+
+    it "returns the Sun's rising azimuth on 2015-02-05" do
+      date = Date.new(2015, 2, 5)
+      epoch = Astronoby::Epoch.from_time(date)
+      observer = Astronoby::Observer.new(
+        latitude: Astronoby::Angle.as_degrees(38),
+        longitude: Astronoby::Angle.as_degrees(-78)
+      )
+      sun = described_class.new(epoch: epoch)
+
+      rising_azimuth = sun.rising_azimuth(observer: observer)
+
+      expect(rising_azimuth&.str(:dms)).to eq "+109° 41′ 24.0917″"
+      # Time from IMCCE: +109° 53′
+    end
+
+    it "returns the Sun's rising azimuth on 1986-03-10" do
+      date = Date.new(1986, 3, 10)
+      epoch = Astronoby::Epoch.from_time(date)
+      observer = Astronoby::Observer.new(
+        latitude: Astronoby::Angle.as_degrees(42.37),
+        longitude: Astronoby::Angle.as_degrees(-71.05)
+      )
+      sun = described_class.new(epoch: epoch)
+
+      rising_azimuth = sun.rising_azimuth(observer: observer)
+
+      expect(rising_azimuth&.str(:dms)).to eq "+94° 59′ 15.7852″"
+      # Time from IMCCE: +95° 02′
+    end
+
+    it "returns the Sun's rising azimuth on 1991-03-14" do
+      date = Date.new(1991, 3, 14)
+      epoch = Astronoby::Epoch.from_time(date)
+      observer = Astronoby::Observer.new(
+        latitude: Astronoby::Angle.as_degrees(48.8566),
+        longitude: Astronoby::Angle.as_degrees(2.3522)
+      )
+      sun = described_class.new(epoch: epoch)
+
+      rising_azimuth = sun.rising_azimuth(observer: observer)
+
+      expect(rising_azimuth&.str(:dms)).to eq "+93° 26′ 26.8564″"
+      # Time from IMCCE: +93° 26′
+    end
+  end
+
   describe "#setting_time" do
     it "returns a time" do
       date = Date.new
@@ -450,4 +511,65 @@
       # Time from IMCCE: 1991-03-14T17:52:00
     end
   end
+
+  describe "#setting_azimuth" do
+    it "returns an Angle" do
+      date = Date.new
+      epoch = Astronoby::Epoch.from_time(date)
+      observer = Astronoby::Observer.new(
+        latitude: Astronoby::Angle.zero,
+        longitude: Astronoby::Angle.zero
+      )
+      sun = described_class.new(epoch: epoch)
+
+      setting_azimuth = sun.setting_azimuth(observer: observer)
+
+      expect(setting_azimuth).to be_a(Astronoby::Angle)
+    end
+
+    it "returns the Sun's setting azimuth on 2015-02-05" do
+      date = Date.new(2015, 2, 5)
+      epoch = Astronoby::Epoch.from_time(date)
+      observer = Astronoby::Observer.new(
+        latitude: Astronoby::Angle.as_degrees(38),
+        longitude: Astronoby::Angle.as_degrees(-78)
+      )
+      sun = described_class.new(epoch: epoch)
+
+      setting_azimuth = sun.setting_azimuth(observer: observer)
+
+      expect(setting_azimuth&.str(:dms)).to eq "+250° 18′ 35.9082″"
+      # Time from IMCCE: +250° 18′
+    end
+
+    it "returns the Sun's setting azimuth on 1986-03-10" do
+      date = Date.new(1986, 3, 10)
+      epoch = Astronoby::Epoch.from_time(date)
+      observer = Astronoby::Observer.new(
+        latitude: Astronoby::Angle.as_degrees(42.37),
+        longitude: Astronoby::Angle.as_degrees(-71.05)
+      )
+      sun = described_class.new(epoch: epoch)
+
+      setting_azimuth = sun.setting_azimuth(observer: observer)
+
+      expect(setting_azimuth&.str(:dms)).to eq "+265° 0′ 44.2147″"
+      # Time from IMCCE: +265° 14′
+    end
+
+    it "returns the Sun's setting azimuth on 1991-03-14" do
+      date = Date.new(1991, 3, 14)
+      epoch = Astronoby::Epoch.from_time(date)
+      observer = Astronoby::Observer.new(
+        latitude: Astronoby::Angle.as_degrees(48.8566),
+        longitude: Astronoby::Angle.as_degrees(2.3522)
+      )
+      sun = described_class.new(epoch: epoch)
+
+      setting_azimuth = sun.setting_azimuth(observer: observer)
+
+      expect(setting_azimuth&.str(:dms)).to eq "+266° 33′ 33.1435″"
+      # Time from IMCCE: +266° 52′
+    end
+  end
 end