Introduce MeanOfDate reference frame

lib/astronoby.rb

@@ -44,6 +44,7 @@
 require "astronoby/reference_frame"
 require "astronoby/reference_frames/geometric"
 require "astronoby/reference_frames/astrometric"
+require "astronoby/reference_frames/mean_of_date"
 require "astronoby/refraction"
 require "astronoby/time/greenwich_sidereal_time"
 require "astronoby/time/local_sidereal_time"

lib/astronoby/bodies/earth.rb

@@ -28,5 +28,23 @@ def compute_astrometric(ephem)
         target_body: self.class
       )
     end
+
+    def compute_mean_of_date(ephem)
+      MeanOfDate.new(
+        position: Vector[
+          Distance.zero,
+          Distance.zero,
+          Distance.zero
+        ],
+        velocity: Vector[
+          Velocity.zero,
+          Velocity.zero,
+          Velocity.zero
+        ],
+        instant: @instant,
+        center_identifier: EARTH,
+        target_body: self.class
+      )
+    end
   end
 end

lib/astronoby/bodies/solar_system_body.rb

@@ -17,7 +17,7 @@ class SolarSystemBody
     URANUS_BARYCENTER = 7
     NEPTUNE_BARYCENTER = 8
 
-    attr_reader :geometric, :astrometric, :instant
+    attr_reader :geometric, :astrometric, :mean_of_date, :instant
 
     def self.geometric(ephem:, instant:)
       compute_geometric(ephem: ephem, instant: instant)
@@ -67,8 +67,10 @@ def self.ephemeris_segments
 
     def initialize(ephem:, instant:)
       @instant = instant
+      # TODO: Compute only values that depend on ephem and lazy load the rest
       @geometric = compute_geometric(ephem)
       @astrometric = compute_astrometric(ephem)
+      @mean_of_date = compute_mean_of_date(ephem)
     end
 
     private
@@ -85,5 +87,14 @@ def compute_astrometric(ephem)
         target_body: self
       )
     end
+
+    def compute_mean_of_date(ephem)
+      MeanOfDate.build_from_geometric(
+        ephem: ephem,
+        instant: @instant,
+        target_geometric: @geometric,
+        target_body: self
+      )
+    end
   end
 end

lib/astronoby/precession.rb

@@ -2,7 +2,11 @@
 
 module Astronoby
   class Precession
-    def initialize(instant)
+    def self.matrix_for(instant)
+      new(instant: instant).matrix
+    end
+
+    def initialize(instant:)
       @instant = instant
     end
 

lib/astronoby/reference_frames/astrometric.rb

@@ -15,7 +15,7 @@ def self.build_from_geometric(
           target: target_geometric,
           ephem: ephem
         )
-      Astrometric.new(
+      new(
         position: corrected_position - earth_geometric.position,
         velocity: corrected_velocity - earth_geometric.velocity,
         instant: instant,

lib/astronoby/reference_frames/mean_of_date.rb

@@ -0,0 +1,38 @@
+# frozen_string_literal: true
+
+module Astronoby
+  class MeanOfDate < ReferenceFrame
+    def self.build_from_geometric(
+      ephem:,
+      instant:,
+      target_geometric:,
+      target_body:
+    )
+      earth_geometric = Earth.geometric(ephem: ephem, instant: instant)
+      position = target_geometric.position - earth_geometric.position
+      velocity = target_geometric.velocity - earth_geometric.velocity
+      precession_matrix = Precession.matrix_for(instant)
+      corrected_position = Vector
+        .elements(precession_matrix * position.map(&:m))
+        .map { Distance.from_meters(_1) }
+      corrected_velocity = Vector[*precession_matrix * velocity
+        .map(&:aupd)].map { Velocity.from_astronomical_units_per_day(_1) }
+
+      new(
+        position: corrected_position,
+        velocity: corrected_velocity,
+        instant: instant,
+        center_identifier: SolarSystemBody::EARTH,
+        target_body: target_body
+      )
+    end
+
+    def ecliptic
+      @ecliptic ||= begin
+        return Coordinates::Ecliptic.zero if distance.zero?
+
+        equatorial.to_ecliptic(epoch: @instant.tdb)
+      end
+    end
+  end
+end

spec/astronoby/bodies/earth_spec.rb

@@ -137,4 +137,53 @@
         .to eq([0, 0, 0])
     end
   end
+
+  describe "#mean_of_date" do
+    it "returns a MeanOfDate position" do
+      time = Time.utc(2025, 2, 7, 12)
+      instant = Astronoby::Instant.from_time(time)
+      state = double(
+        position: Ephem::Core::Vector[1, 2, 3],
+        velocity: Ephem::Core::Vector[4, 5, 6]
+      )
+      segment = double(compute_and_differentiate: state)
+      ephem = double(:[] => segment)
+      planet = described_class.new(instant: instant, ephem: ephem)
+
+      mean_of_date = planet.mean_of_date
+
+      expect(mean_of_date).to be_a(Astronoby::MeanOfDate)
+      expect(mean_of_date.equatorial)
+        .to be_a(Astronoby::Coordinates::Equatorial)
+      expect(mean_of_date.distance).to be_a(Astronoby::Distance)
+    end
+
+    it "computes the correct position" do
+      time = Time.utc(2025, 4, 1)
+      instant = Astronoby::Instant.from_time(time)
+      ephem = test_ephem
+      planet = described_class.new(instant: instant, ephem: ephem)
+
+      mean_of_date = planet.mean_of_date
+
+      expect(mean_of_date.equatorial.right_ascension)
+        .to eq(Astronoby::Angle.zero)
+      expect(mean_of_date.equatorial.declination).to eq(Astronoby::Angle.zero)
+      expect(mean_of_date.ecliptic.latitude).to eq(Astronoby::Angle.zero)
+      expect(mean_of_date.ecliptic.longitude).to eq(Astronoby::Angle.zero)
+      expect(mean_of_date.distance).to eq(Astronoby::Distance.zero)
+    end
+
+    it "computes the correct velocity" do
+      time = Time.utc(2025, 3, 1)
+      instant = Astronoby::Instant.from_time(time)
+      ephem = test_ephem
+      planet = described_class.new(instant: instant, ephem: ephem)
+
+      mean_of_date = planet.mean_of_date
+
+      expect(mean_of_date.velocity.to_a.map(&:mps).map { _1.round(5) })
+        .to eq([0, 0, 0])
+    end
+  end
 end

spec/astronoby/bodies/jupiter_spec.rb

@@ -157,4 +157,67 @@
       # Skyfield: -32111.8691  21668.91368 9711.35459
     end
   end
+
+  describe "#mean_of_date" do
+    it "returns a MeanOfDate position" do
+      time = Time.utc(2025, 2, 7, 12)
+      instant = Astronoby::Instant.from_time(time)
+      state = double(
+        position: Ephem::Core::Vector[1, 2, 3],
+        velocity: Ephem::Core::Vector[4, 5, 6]
+      )
+      segment = double(compute_and_differentiate: state)
+      ephem = double(:[] => segment)
+      planet = described_class.new(instant: instant, ephem: ephem)
+
+      mean_of_date = planet.mean_of_date
+
+      expect(mean_of_date).to be_a(Astronoby::MeanOfDate)
+      expect(mean_of_date.equatorial).to be_a(Astronoby::Coordinates::Equatorial)
+      expect(mean_of_date.ecliptic).to be_a(Astronoby::Coordinates::Ecliptic)
+      expect(mean_of_date.distance).to be_a(Astronoby::Distance)
+    end
+
+    it "computes the correct position" do
+      time = Time.utc(2025, 5, 1)
+      instant = Astronoby::Instant.from_time(time)
+      ephem = test_ephem
+      planet = described_class.new(instant: instant, ephem: ephem)
+
+      mean_of_date = planet.mean_of_date
+
+      expect(mean_of_date.equatorial.right_ascension.str(:hms))
+        .to eq("5h 22m 31.431s")
+      # IMCCE:  5h 22m 31.4319s
+
+      expect(mean_of_date.equatorial.declination.str(:dms))
+        .to eq("+22° 55′ 5.8692″")
+      # IMCCE:  +22° 55′ 5.858″
+
+      expect(mean_of_date.ecliptic.latitude.str(:dms))
+        .to eq("-0° 14′ 15.6927″")
+      # IMCCE:  -0° 14′ 15.705″
+
+      expect(mean_of_date.ecliptic.longitude.str(:dms))
+        .to eq("+81° 22′ 34.8729″")
+      # IMCCE:  +81° 22′ 34.883″
+
+      expect(mean_of_date.distance.au)
+        .to eq(5.847671713145707)
+      # IMCCE: 5.847671760041
+    end
+
+    it "computes the correct velocity" do
+      time = Time.utc(2025, 5, 1)
+      instant = Astronoby::Instant.from_time(time)
+      ephem = test_ephem
+      planet = described_class.new(instant: instant, ephem: ephem)
+
+      mean_of_date = planet.mean_of_date
+
+      expect(mean_of_date.velocity.to_a.map(&:mps).map { _1.round(5) })
+        .to eq([-32257.90536, 21486.0083, 9631.88754])
+      # IMCCE:  -32257.90781  21486.0049  9631.88872
+    end
+  end
 end

spec/astronoby/bodies/mars_spec.rb

@@ -157,4 +157,67 @@
       # Skyfield: -16789.69861 9770.68797 4114.24472
     end
   end
+
+  describe "#mean_of_date" do
+    it "returns a MeanOfDate position" do
+      time = Time.utc(2025, 2, 7, 12)
+      instant = Astronoby::Instant.from_time(time)
+      state = double(
+        position: Ephem::Core::Vector[1, 2, 3],
+        velocity: Ephem::Core::Vector[4, 5, 6]
+      )
+      segment = double(compute_and_differentiate: state)
+      ephem = double(:[] => segment)
+      planet = described_class.new(instant: instant, ephem: ephem)
+
+      mean_of_date = planet.mean_of_date
+
+      expect(mean_of_date).to be_a(Astronoby::MeanOfDate)
+      expect(mean_of_date.equatorial).to be_a(Astronoby::Coordinates::Equatorial)
+      expect(mean_of_date.ecliptic).to be_a(Astronoby::Coordinates::Ecliptic)
+      expect(mean_of_date.distance).to be_a(Astronoby::Distance)
+    end
+
+    it "computes the correct position" do
+      time = Time.utc(2025, 4, 1)
+      instant = Astronoby::Instant.from_time(time)
+      ephem = test_ephem
+      planet = described_class.new(instant: instant, ephem: ephem)
+
+      mean_of_date = planet.mean_of_date
+
+      expect(mean_of_date.equatorial.right_ascension.str(:hms))
+        .to eq("7h 43m 51.3343s")
+      # IMCCE:  7h 43m 51.3351s
+
+      expect(mean_of_date.equatorial.declination.str(:dms))
+        .to eq("+23° 59′ 54.0557″")
+      # IMCCE:  +23° 59′ 54.056″
+
+      expect(mean_of_date.ecliptic.latitude.str(:dms))
+        .to eq("+2° 39′ 52.0958″")
+      # IMCCE:  +2° 39′ 52.098″
+
+      expect(mean_of_date.ecliptic.longitude.str(:dms))
+        .to eq("+113° 36′ 9.7732″")
+      # IMCCE:  +113° 36′ 9.784″
+
+      expect(mean_of_date.distance.au)
+        .to eq(1.1389410337239916)
+      # IMCCE: 1.1389410695
+    end
+
+    it "computes the correct velocity" do
+      time = Time.utc(2025, 4, 1)
+      instant = Astronoby::Instant.from_time(time)
+      ephem = test_ephem
+      planet = described_class.new(instant: instant, ephem: ephem)
+
+      mean_of_date = planet.mean_of_date
+
+      expect(mean_of_date.velocity.to_a.map(&:mps).map { _1.round(5) })
+        .to eq([-16853.58684, 9675.16449, 4072.69891])
+      # IMCCE:  -16853.58852  9675.16272  4072.69963
+    end
+  end
 end

spec/astronoby/bodies/mercury_spec.rb

@@ -157,4 +157,67 @@
       # Skyfield: 38473.18754 -32529.64572 -18788.09165
     end
   end
+
+  describe "#mean_of_date" do
+    it "returns a MeanOfDate position" do
+      time = Time.utc(2025, 2, 7, 12)
+      instant = Astronoby::Instant.from_time(time)
+      state = double(
+        position: Ephem::Core::Vector[1, 2, 3],
+        velocity: Ephem::Core::Vector[4, 5, 6]
+      )
+      segment = double(compute_and_differentiate: state)
+      ephem = double(:[] => segment)
+      planet = described_class.new(instant: instant, ephem: ephem)
+
+      mean_of_date = planet.mean_of_date
+
+      expect(mean_of_date).to be_a(Astronoby::MeanOfDate)
+      expect(mean_of_date.equatorial).to be_a(Astronoby::Coordinates::Equatorial)
+      expect(mean_of_date.ecliptic).to be_a(Astronoby::Coordinates::Ecliptic)
+      expect(mean_of_date.distance).to be_a(Astronoby::Distance)
+    end
+
+    it "computes the correct position" do
+      time = Time.utc(2025, 1, 1)
+      instant = Astronoby::Instant.from_time(time)
+      ephem = test_ephem
+      planet = described_class.new(instant: instant, ephem: ephem)
+
+      mean_of_date = planet.mean_of_date
+
+      expect(mean_of_date.equatorial.right_ascension.str(:hms))
+        .to eq("17h 16m 19.4051s")
+      # IMCCE:  17h 16m 19.4064s
+
+      expect(mean_of_date.equatorial.declination.str(:dms))
+        .to eq("-21° 56′ 26.4339″")
+      # IMCCE:  -21° 56′ 26.426″
+
+      expect(mean_of_date.ecliptic.latitude.str(:dms))
+        .to eq("+1° 6′ 45.242″")
+      # IMCCE:  +1° 6′ 45.252″
+
+      expect(mean_of_date.ecliptic.longitude.str(:dms))
+        .to eq("+259° 52′ 42.4276″")
+      # IMCCE:  +259° 52′ 42.445″
+
+      expect(mean_of_date.distance.au)
+        .to eq(1.1480561241626033)
+      # IMCCE: 1.148056160668
+    end
+
+    it "computes the correct velocity" do
+      time = Time.utc(2025, 1, 1)
+      instant = Astronoby::Instant.from_time(time)
+      ephem = test_ephem
+      planet = described_class.new(instant: instant, ephem: ephem)
+
+      mean_of_date = planet.mean_of_date
+
+      expect(mean_of_date.velocity.to_a.map(&:mps).map { _1.round(5) })
+        .to eq([38717.64165, -32306.63286, -18692.31286])
+      # IMCCE:  38717.64791  -32306.62736  -18692.31383
+    end
+  end
 end