Introduce MeanOfDate reference frame (#138)

> The mean coordinates of the date are related to the ecliptic or the mean equator and the mean equinox of the date. These are the geometric coordinates corrected for the precessional motion of the Earth’s rotation axis.

This introduces `MeanOfDate` reference frame, from `Geometric`, with precession of the date applied. By "date", we actually mean "instant", although the precession only changes slightly from one day to the next.

In Astronoby, the `MeanOfDate` reference frame is geocentric (like `Astrometric`). Mean of date coordinates don't _have to_ be geocentric, but that's the design choice for Astronoby at the moment.

Like `Geometric` and `Astrometry`, the `MeanOfDate` reference frame offers equatorial and ecliptic coordinates, among with the Cartesian position and velocity, and the distance to the Earth's center.

```rb
instant = Astronoby::Instant.from_time(Time.utc(2025, 2, 7))
ephem = Astronoby::Ephem.load("de440s.bsp")

sun = Astronoby::Sun.new(instant: instant, ephem: ephem)

sun.mean_of_date.position.map(&:km).map(&:round)
# => Vector[110436790, -89767243, -38912451]

sun.mean_of_date.velocity.map(&:kmpd).map(&:round)
# => Vector[1747570, 1776918, 770305]

sun.mean_of_date.equatorial.right_ascension.str(:hms)
# => "21h 23m 34.6729s"

sun.mean_of_date.equatorial.declination.str(:dms)
# => "-15° 17′ 31.1747″"

sun.mean_of_date.ecliptic.latitude.str(:dms)
# => "+0° 0′ 0.4003″"

sun.mean_of_date.ecliptic.longitude.str(:dms)
# => "+318° 27′ 41.6839″"

sun.mean_of_date.distance.km.round
# => 147541931
```

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