Add binary PCK orientation kernel support

CHANGELOG.md

@@ -1,5 +1,37 @@
 # Changelog
 
+## [0.5.0] - 2026-06-22
+
+### Features
+
+* Read binary PCK (`DAF/PCK`) orientation kernels via `Ephem::PCK`, exposing a
+  body's Euler angles and rates over time (`angles_at`, `orientation_at`), the
+  foundation for DE440-grade lunar libration ([#76])
+* Add `Ephem::Core::Orientation` (Euler angles + optional rates) and
+  `Ephem::Core::Rotation` (kernel-agnostic rotation-matrix helpers), plus
+  `Orientation#to_matrix` / `OrientationSegment#matrix_at` for the built-in
+  3-1-3 (Z-X-Z) reference→body convention
+* Excerpt and the `excerpt` CLI now support binary PCK kernels, detecting the
+  kernel kind automatically
+* Download binary PCK lunar orientation kernels from NAIF via `Ephem::Download`
+
+### Improvements
+
+* Route queries to the covering segment when a body/pair spans multiple
+  time-split segments (SPK and PCK), with no overhead for single-segment keys
+* Share the type-2 Chebyshev machinery between SPK and PCK segments
+* Evaluate position and velocity in a single Chebyshev pass
+  (`ChebyshevPolynomial.evaluate_with_derivative`), speeding up every state /
+  orientation query (`compute_and_differentiate`, `state_at`, `orientation_at`)
+  with bit-for-bit identical results
+* Fix `compute_and_differentiate` returning mismatched velocities for an array
+  of times
+* Correct the documented velocity unit to km/day (the actual, validated value)
+
+[#76]: https://github.com/rhannequin/ruby-ephem/issues/76
+
+**Full Changelog**: https://github.com/rhannequin/ruby-ephem/compare/v0.4.1...v0.5.0
+
 ## [0.4.1] - 2025-08-03
 
 ### Improvements

README.md

@@ -12,8 +12,9 @@ These files are a collection of numerical integrations of the equations of
 motion of the Solar System, used to calculate the positions of the planets,
 the Moon, and other celestial bodies with high precision.
 
-Ephem currently only support planetary ephemerides like DE421, DE430,
-INPOP19A, etc.
+Ephem supports planetary ephemerides like DE421, DE430, INPOP19A, etc. (SPK),
+as well as binary PCK orientation kernels such as the lunar orientation kernel
+`moon_pa_de440_200625.bpc`.
 
 The library is in high development mode and does not have a stable version yet.
 The API is subject to major changes at the moment, please keep that in mind if
@@ -99,6 +100,67 @@ puts "Velocity: #{state.velocity}"
 # The velocity is expressed in km/day
 ```
 
+## Orientation kernels (binary PCK)
+
+Binary PCK (`DAF/PCK`) kernels store the orientation of a body frame over time
+as Euler angles, rather than position. The most common use is the Moon's
+physical libration via a lunar orientation kernel such as
+`moon_pa_de440_200625.bpc`.
+
+```rb
+# Download and load a binary PCK orientation kernel
+Ephem::Download.call(
+  name: "moon_pa_de440_200625.bpc",
+  target: "tmp/moon_pa_de440_200625.bpc"
+)
+pck = Ephem::PCK.open("tmp/moon_pa_de440_200625.bpc")
+
+puts pck
+# => PCK file with 2 segments:
+# 1549-12-21..2426-02-16 Type 2 orientation of frame 31008 relative to frame 1
+# 2426-02-16..2650-01-25 Type 2 orientation of frame 31008 relative to frame 1
+
+# Query a body frame by its NAIF frame ID (31008 = MOON_PA_DE440).
+# The time is expressed in Julian Date (TDB).
+frame = pck[31008]
+orientation = frame.orientation_at(2451545.0)
+
+# The three classical 3-1-3 (Z-X-Z) Euler angles, in radians
+puts "phi:   #{orientation.phi}"
+puts "theta: #{orientation.theta}"
+puts "psi:   #{orientation.psi}"
+
+# Rates are in radians/day (ephem's per-day convention, like SPK velocities)
+puts "rates: #{orientation.rates.inspect}"
+
+# Use #angles_at when you do not need the rates (it skips the derivative)
+angles = frame.angles_at(2451545.0)
+
+pck.close
+```
+
+### Building the body-fixed frame
+
+The orientation maps the reference frame into the body-fixed frame through the
+3-1-3 (Z-X-Z) Euler sequence `M = Rz(psi) * Rx(theta) * Rz(phi)`. That
+convention is built in, so you can get the rotation matrix directly:
+
+```rb
+# Straight from a frame, at a given time
+matrix = frame.matrix_at(2451545.0)
+
+# ...or from an Orientation you already have
+matrix = orientation.to_matrix
+
+# Project a geocentric Earth->Moon vector into the body-fixed frame
+body_fixed = Ephem::Core::Rotation.apply(matrix, earth_to_moon)
+```
+
+`Ephem::Core::Rotation` also exposes the lower-level, kernel-agnostic building
+blocks (`about_x`, `about_y`, `about_z`, `multiply`, `apply`) if you need a
+different axis order or an additional fixed rotation (such as the PA -> ME
+offset, whose constants the consumer supplies).
+
 ## CLI
 
 The gem also provides a CLI to generate an excerpt from an original kernel file.
@@ -121,6 +183,13 @@ While DE440s originally supports 14 segments from 1849 to 2150 with a size of
 Not only the excerpt is smaller, but most importantly it is way more
 efficient to parse and to use in your application.
 
+The input kernel kind (SPK or binary PCK) is detected automatically, so the
+same command excerpts orientation kernels too:
+
+```bash
+ruby-ephem excerpt --targets 31008 2000-01-01 2030-01-01 /path/to/moon_pa_de440.bpc moon_excerpt.bpc
+```
+
 ## Accuracy
 
 Data from this library has been tested against the Python library [jplephem]

benchmarks/run.rb

@@ -12,6 +12,8 @@
 SPK_FULL = File.join(ROOT, "spec", "support", "data", "de432s.bsp")
 SPK_EXCERPT = File
   .join(ROOT, "spec", "support", "data", "de421_2000_excerpt.bsp")
+PCK_EXCERPT = File
+  .join(ROOT, "spec", "support", "data", "moon_pa_de440_excerpt.bpc")
 
 JD_J2000 = Ephem::Core::Constants::Time::J2000_EPOCH # 2451545.0
 JD_TEST = 2459000.0 # 2020-09-30, well within de432s range
@@ -293,6 +295,12 @@ def ensure_file!(path)
     )
   end
 
+  x.report("evaluate_with_derivative (pos+vel, 1 pass)") do
+    Ephem::Computation::ChebyshevPolynomial.evaluate_with_derivative(
+      test_coeffs, test_t, test_radius
+    )
+  end
+
   x.compare!
 end
 
@@ -388,6 +396,32 @@ def ensure_file!(path)
   puts "  Components per term:             #{coeffs.first&.first&.size}"
 end
 
+# =========================================================================
+# 13. PCK ORIENTATION (BINARY PCK)
+# =========================================================================
+
+separator "13. PCK Orientation (Binary PCK)"
+
+puts "angles_at (Euler angles) vs orientation_at (angles + rates)"
+puts
+
+pck = Ephem::PCK.open(PCK_EXCERPT)
+moon = pck[31008] # MOON_PA_DE440 frame
+moon.angles_at(JD_J2000) # warm up
+
+orientation_times = Array.new(SEQUENTIAL_STEPS) { |index| JD_J2000 + index }
+
+GC.start
+Benchmark.ips do |x|
+  x.report("angles_at (scalar)") { moon.angles_at(JD_J2000) }
+  x.report("orientation_at (scalar)") { moon.orientation_at(JD_J2000) }
+  x.report("angles_at (batch 1000)") { moon.angles_at(orientation_times) }
+
+  x.compare!
+end
+
+pck.close
+
 # =========================================================================
 # Cleanup
 # =========================================================================

lib/ephem.rb

@@ -6,6 +6,8 @@
 require_relative "ephem/core/calendar_calculations"
 require_relative "ephem/core/state"
 require_relative "ephem/core/vector"
+require_relative "ephem/core/orientation"
+require_relative "ephem/core/rotation"
 require_relative "ephem/error"
 require_relative "ephem/io/binary_reader"
 require_relative "ephem/io/daf"
@@ -15,9 +17,16 @@
 require_relative "ephem/io/record_parser"
 require_relative "ephem/io/summary_manager"
 require_relative "ephem/spk"
+require_relative "ephem/pck"
 require_relative "ephem/segments/base_segment"
 require_relative "ephem/segments/registry"
+require_relative "ephem/segments/chebyshev_type2"
+require_relative "ephem/segments/orientation_source"
 require_relative "ephem/segments/segment"
+require_relative "ephem/segments/orientation_segment"
+require_relative "ephem/segments/segment_group"
+require_relative "ephem/segments/position_group"
+require_relative "ephem/segments/orientation_group"
 require_relative "ephem/excerpt"
 require_relative "ephem/cli"
 require_relative "ephem/version"

lib/ephem/cli.rb

@@ -47,22 +47,22 @@ def self.show_help
         Ruby Ephem - A tool for working with JPL Ephemerides
 
         Commands:
-          excerpt - Create an excerpt of an SPK file
+          excerpt - Create an excerpt of an SPK or binary PCK kernel
           help    - Show this help message
 
         Excerpt command:
           ruby-ephem excerpt [options] START_DATE END_DATE INPUT_FILE OUTPUT_FILE
 
         Options:
-          --targets TARGET_IDS  - Comma-separated list of target IDs to include
-                                  (default: all targets)
+          --targets TARGET_IDS  - Comma-separated list of target/body IDs to
+                                  include (default: all)
 
-        Example:
+        Examples:
           ruby-ephem excerpt --targets 3,10,399 2000-01-01 2030-01-01 de440s.bsp excerpt.bsp
+          ruby-ephem excerpt --targets 31008 2000-01-01 2030-01-01 moon_pa_de440.bpc moon_excerpt.bpc
 
-        This will create an excerpt of de440s.bsp containing only the specified
-        targets (Earth-Moon barycenter, Sun, Earth) for the period from
-        2000-01-01 to 2030-01-01.
+        The input kernel kind (SPK or binary PCK) is detected automatically and
+        the excerpt is written in the same format.
       HELP
     end
 
@@ -125,9 +125,9 @@ def self.handle_excerpt(args)
           puts "Including all targets"
         end
 
-        spk = Ephem::SPK.open(input_file)
+        kernel = open_kernel(input_file)
 
-        excerpt_spk = spk.excerpt(
+        excerpt_kernel = kernel.excerpt(
           output_path: output_file,
           start_jd: start_jd,
           end_jd: end_jd,
@@ -136,8 +136,8 @@ def self.handle_excerpt(args)
         )
 
         puts "Excerpt created successfully!"
-        puts "Original segments: #{spk.segments.size}"
-        puts "Excerpt segments: #{excerpt_spk.segments.size}"
+        puts "Original segments: #{kernel.segments.size}"
+        puts "Excerpt segments: #{excerpt_kernel.segments.size}"
 
         original_size = File.size(input_file)
         excerpt_size = File.size(output_file)
@@ -148,12 +148,25 @@ def self.handle_excerpt(args)
         puts "Original: #{original_size} bytes"
         puts "Excerpt: #{excerpt_size} bytes"
 
-        spk.close
-        excerpt_spk.close
+        kernel.close
+        excerpt_kernel.close
       rescue => e
         puts "Error creating excerpt: #{e.message}"
         puts e.backtrace if options[:debug]
       end
     end
+
+    def self.open_kernel(path)
+      daf = Ephem::IO::DAF.new(File.open(path, "rb"))
+
+      if daf.file_type == :pck
+        Ephem::PCK.new(daf: daf)
+      else
+        Ephem::SPK.new(daf: daf)
+      end
+    rescue
+      daf&.close
+      raise
+    end
   end
 end

lib/ephem/computation/chebyshev_polynomial.rb

@@ -85,6 +85,65 @@ def self.evaluate_derivative(coeffs, t, radius)
         scale = Ephem::Core::Constants::Time::SECONDS_PER_DAY / (2.0 * radius)
         [d1x * scale, d1y * scale, d1z * scale]
       end
+
+      ##
+      # Evaluates a 3D Chebyshev polynomial and its time derivative in a single
+      # pass. It runs the same value and derivative recurrences as {evaluate}
+      # and {evaluate_derivative}, but fused into one loop so the coefficient
+      # fetch and loop control are shared. Results are bit-for-bit identical to
+      # calling the two methods separately.
+      #
+      # @param coeffs [Array<Array<Float>>] coefficients; shape [n_terms][3].
+      # @param t [Float] normalized independent variable, in [-1, 1].
+      # @param radius [Float] half-length of the time interval (seconds).
+      # @return [Array(Array<Float>, Array<Float>)] [position, velocity], with
+      #   velocity in units per day.
+      def self.evaluate_with_derivative(coeffs, t, radius)
+        n = coeffs.size
+        b1x = b1y = b1z = 0.0
+        b2x = b2y = b2z = 0.0
+        d1x = d1y = d1z = 0.0
+        d2x = d2y = d2z = 0.0
+
+        t2 = 2.0 * t
+        k = n - 1
+        while k > 0
+          c = coeffs[k]
+          c0 = c[0]
+          c1 = c[1]
+          c2 = c[2]
+          k2 = 2 * k
+
+          bx = t2 * b1x - b2x + c0
+          by = t2 * b1y - b2y + c1
+          bz = t2 * b1z - b2z + c2
+          dx = t2 * d1x - d2x + k2 * c0
+          dy = t2 * d1y - d2y + k2 * c1
+          dz = t2 * d1z - d2z + k2 * c2
+
+          b2x = b1x
+          b2y = b1y
+          b2z = b1z
+          b1x = bx
+          b1y = by
+          b1z = bz
+          d2x = d1x
+          d2y = d1y
+          d2z = d1z
+          d1x = dx
+          d1y = dy
+          d1z = dz
+          k -= 1
+        end
+
+        c0, c1, c2 = coeffs[0]
+        position = [t * b1x - b2x + c0, t * b1y - b2y + c1, t * b1z - b2z + c2]
+
+        scale = Ephem::Core::Constants::Time::SECONDS_PER_DAY / (2.0 * radius)
+        velocity = [d1x * scale, d1y * scale, d1z * scale]
+
+        [position, velocity]
+      end
     end
   end
 end