Add GreenwichSiderealTime and LocalSiderealTime

lib/astronoby.rb

@@ -18,8 +18,9 @@
 require "astronoby/observer"
 require "astronoby/precession"
 require "astronoby/refraction"
+require "astronoby/time/greenwich_sidereal_time"
+require "astronoby/time/local_sidereal_time"
 require "astronoby/util/astrodynamics"
-require "astronoby/util/time"
 require "astronoby/util/trigonometry"
 require "astronoby/true_obliquity"
 require "astronoby/version"

lib/astronoby/body.rb

@@ -19,14 +19,13 @@ def rising_time(latitude:, longitude:, date:, apparent: true)
       return nil unless RISING_SETTING_HOUR_ANGLE_RATIO_RANGE.cover?(ratio)
 
       hour_angle = Angle.acos(ratio)
-      local_sidereal_time =
-        @equatorial_coordinates.right_ascension.hours - hour_angle.hours
-
-      Util::Time.lst_to_ut(
+      local_sidereal_time = LocalSiderealTime.new(
         date: date,
-        longitude: longitude,
-        lst: local_sidereal_time
+        time: @equatorial_coordinates.right_ascension.hours - hour_angle.hours,
+        longitude: longitude
       )
+
+      local_sidereal_time.to_gst.to_utc
     end
 
     # Source:
@@ -51,14 +50,13 @@ def setting_time(latitude:, longitude:, date:, apparent: true)
       return nil unless RISING_SETTING_HOUR_ANGLE_RATIO_RANGE.cover?(ratio)
 
       hour_angle = Angle.acos(ratio)
-      local_sidereal_time =
-        @equatorial_coordinates.right_ascension.hours + hour_angle.hours
-
-      Util::Time.lst_to_ut(
+      local_sidereal_time = LocalSiderealTime.new(
         date: date,
-        longitude: longitude,
-        lst: local_sidereal_time
+        time: @equatorial_coordinates.right_ascension.hours + hour_angle.hours,
+        longitude: longitude
       )
+
+      local_sidereal_time.to_gst.to_utc
     end
 
     # Source:

lib/astronoby/coordinates/equatorial.rb

@@ -18,12 +18,11 @@ def initialize(
       end
 
       def compute_hour_angle(time:, longitude:)
-        lst = Util::Time.local_sidereal_time(
-          time: time,
-          longitude: longitude
-        )
+        lst = GreenwichSiderealTime
+          .from_utc(time.utc)
+          .to_lst(longitude: longitude)
 
-        ha = (lst - @right_ascension.hours)
+        ha = (lst.time - @right_ascension.hours)
         ha += 24 if ha.negative?
 
         Angle.as_hours(ha)

lib/astronoby/coordinates/horizontal.rb

@@ -36,10 +36,10 @@ def to_equatorial(time:)
         end
 
         hour_angle_hours = Angle.as_degrees(hour_angle_degrees).hours
-        right_ascension_decimal = Util::Time.local_sidereal_time(
-          time: time,
-          longitude: @longitude
-        ) - hour_angle_hours
+        lst = GreenwichSiderealTime
+          .from_utc(time.utc)
+          .to_lst(longitude: @longitude)
+        right_ascension_decimal = lst.time - hour_angle_hours
         right_ascension_decimal += 24 if right_ascension_decimal.negative?
         right_ascension = Angle.as_hours(right_ascension_decimal)
 

lib/astronoby/time/greenwich_sidereal_time.rb

@@ -0,0 +1,86 @@
+# frozen_string_literal: true
+
+module Astronoby
+  class GreenwichSiderealTime
+    JULIAN_CENTURIES_EXPONENTS = [
+      BigDecimal("6.697374558"),
+      BigDecimal("2400.051336"),
+      BigDecimal("0.000025862")
+    ].freeze
+
+    attr_reader :date, :time
+
+    # Source:
+    #  Title: Practical Astronomy with your Calculator or Spreadsheet
+    #  Authors: Peter Duffett-Smith and Jonathan Zwart
+    #  Edition: Cambridge University Press
+    #  Chapter: 12 - Conversion of UT to Greenwich sidereal time (GST)
+    def self.from_utc(utc)
+      date = utc.to_date
+      julian_day = utc.to_date.ajd
+      t = (julian_day - Epoch::J2000) / Epoch::DAYS_PER_JULIAN_CENTURY
+      t0 = (
+        (JULIAN_CENTURIES_EXPONENTS[0] +
+          (JULIAN_CENTURIES_EXPONENTS[1] * t) +
+          (JULIAN_CENTURIES_EXPONENTS[2] * t * t)) % 24
+      ).abs
+
+      ut_in_hours = utc.hour +
+        utc.min / 60.0 +
+        (utc.sec + utc.subsec) / 3600.0
+
+      gmst = BigDecimal("1.002737909") * ut_in_hours + t0
+      gmst += 24 if gmst.negative?
+      gmst -= 24 if gmst > 24
+
+      new(date: date, time: gmst)
+    end
+
+    def initialize(date:, time:)
+      @date = date
+      @time = time
+    end
+
+    # Source:
+    #  Title: Practical Astronomy with your Calculator or Spreadsheet
+    #  Authors: Peter Duffett-Smith and Jonathan Zwart
+    #  Edition: Cambridge University Press
+    #  Chapter: 13 - Conversion of GST to UT
+    def to_utc
+      date = @date
+      julian_day = @date.ajd
+      t = (julian_day - Epoch::J2000) / Epoch::DAYS_PER_JULIAN_CENTURY
+
+      t0 = (
+        (JULIAN_CENTURIES_EXPONENTS[0] +
+          (JULIAN_CENTURIES_EXPONENTS[1] * t) +
+          (JULIAN_CENTURIES_EXPONENTS[2] * t * t)) % 24
+      ).abs
+
+      a = @time - t0
+      a += 24 if a.negative?
+      a -= 24 if a > 24
+
+      utc = BigDecimal("0.9972695663") * a
+
+      decimal_hour_to_time(date, utc)
+    end
+
+    def to_lst(longitude:)
+      LocalSiderealTime.from_gst(gst: self, longitude: longitude)
+    end
+
+    private
+
+    def decimal_hour_to_time(date, decimal)
+      absolute_hour = decimal.abs
+      hour = absolute_hour.floor
+      decimal_minute = 60 * (absolute_hour - hour)
+      absolute_decimal_minute = (60 * (absolute_hour - hour)).abs
+      minute = decimal_minute.floor
+      second = 60 * (absolute_decimal_minute - absolute_decimal_minute.floor)
+
+      ::Time.utc(date.year, date.month, date.day, hour, minute, second)
+    end
+  end
+end

lib/astronoby/time/local_sidereal_time.rb

@@ -0,0 +1,41 @@
+# frozen_string_literal: true
+
+module Astronoby
+  class LocalSiderealTime
+    attr_reader :date, :time, :longitude
+
+    # Source:
+    #  Title: Practical Astronomy with your Calculator or Spreadsheet
+    #  Authors: Peter Duffett-Smith and Jonathan Zwart
+    #  Edition: Cambridge University Press
+    #  Chapter: 14 - Local sidereal time (LST)
+    def self.from_gst(gst:, longitude:)
+      date = gst.date
+      time = gst.time + longitude.hours
+      time += 24 if time.negative?
+      time -= 24 if time > 24
+
+      new(date: date, time: time, longitude: longitude)
+    end
+
+    def initialize(date:, time:, longitude:)
+      @date = date
+      @time = time
+      @longitude = longitude
+    end
+
+    # Source:
+    #  Title: Practical Astronomy with your Calculator or Spreadsheet
+    #  Authors: Peter Duffett-Smith and Jonathan Zwart
+    #  Edition: Cambridge University Press
+    #  Chapter: 15 - Converting LST to GST
+    def to_gst
+      date = @date
+      time = @time - @longitude.hours
+      time += 24 if time.negative?
+      time -= 24 if time > 24
+
+      GreenwichSiderealTime.new(date: date, time: time)
+    end
+  end
+end

spec/astronoby/body_spec.rb

@@ -23,7 +23,7 @@
         apparent: false
       )
 
-      expect(rising_time).to eq(Time.utc(2016, 1, 21, 20, 40, 46))
+      expect(rising_time&.round).to eq(Time.utc(2016, 1, 21, 20, 40, 46))
     end
 
     # Source:
@@ -55,12 +55,13 @@
     #  Edition: MIT Press
     #  Chapter: 5 - Stars in the Nighttime Sky
     it "returns the body's rising time" do
+      offset = -4
       coordinates = Astronoby::Coordinates::Horizontal.new(
         azimuth: Astronoby::Angle.as_degrees(90),
         altitude: Astronoby::Angle.as_degrees(45),
         latitude: Astronoby::Angle.as_degrees(38.25),
         longitude: Astronoby::Angle.as_degrees(-78.3)
-      ).to_equatorial(time: Time.new(2015, 6, 6, 21, 0, 0, "-04:00"))
+      ).to_equatorial(time: Time.new(2015, 6, 6, 21, 0, 0, offset))
       body = described_class.new(coordinates)
 
       rising_time = body.rising_time(
@@ -70,7 +71,8 @@
         apparent: false
       )
 
-      expect(rising_time).to eq(Time.utc(2015, 6, 6, 20, 57, 48))
+      expect(rising_time&.getlocal(offset)&.round)
+        .to eq(Time.new(2015, 6, 6, 16, 57, 48, offset))
     end
 
     # Source:
@@ -91,7 +93,7 @@
         date: Date.new(2010, 8, 24)
       )
 
-      expect(rising_time).to eq(Time.utc(2010, 8, 24, 14, 16, 18))
+      expect(rising_time&.round).to eq(Time.utc(2010, 8, 24, 14, 16, 18))
     end
   end
 
@@ -140,7 +142,7 @@
         apparent: false
       )
 
-      expect(setting_time).to eq(Time.utc(2016, 1, 21, 9, 29, 50))
+      expect(setting_time&.round).to eq(Time.utc(2016, 1, 21, 9, 29, 50))
     end
 
     # Source:
@@ -187,7 +189,7 @@
         apparent: false
       )
 
-      expect(setting_time).to eq(Time.utc(2015, 6, 6, 11, 59, 51))
+      expect(setting_time&.round).to eq(Time.utc(2015, 6, 6, 11, 59, 51))
     end
 
     # Source:
@@ -208,7 +210,7 @@
         date: Date.new(2010, 8, 24)
       )
 
-      expect(setting_time).to eq(Time.utc(2010, 8, 24, 4, 10, 1))
+      expect(setting_time&.round).to eq(Time.utc(2010, 8, 24, 4, 10, 1))
     end
   end