diff --git a/CHANGELOG.md b/CHANGELOG.md index a96ca2f..3bf941d 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/README.md b/README.md index 0e28a44..5c6bdad 100644 --- a/README.md +++ b/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] diff --git a/benchmarks/run.rb b/benchmarks/run.rb index 4483a5c..779f597 100644 --- a/benchmarks/run.rb +++ b/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 # ========================================================================= diff --git a/lib/ephem.rb b/lib/ephem.rb index f31c178..35f7ea5 100644 --- a/lib/ephem.rb +++ b/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" diff --git a/lib/ephem/cli.rb b/lib/ephem/cli.rb index 986e376..79ba96c 100644 --- a/lib/ephem/cli.rb +++ b/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 diff --git a/lib/ephem/computation/chebyshev_polynomial.rb b/lib/ephem/computation/chebyshev_polynomial.rb index c6973fa..3ad3101 100644 --- a/lib/ephem/computation/chebyshev_polynomial.rb +++ b/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>] 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, Array)] [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 diff --git a/lib/ephem/core/orientation.rb b/lib/ephem/core/orientation.rb new file mode 100644 index 0000000..25f4350 --- /dev/null +++ b/lib/ephem/core/orientation.rb @@ -0,0 +1,118 @@ +# frozen_string_literal: true + +module Ephem + module Core + # The orientation of a body, expressed as the three Euler angles that rotate + # a reference frame (e.g. J2000/ICRF) into the body-fixed frame, optionally + # with their time derivatives. + # + # For binary PCK orientation kernels the angles are the classical 3-1-3 + # (Z-X-Z) sequence: +phi+ and +theta+ orient the pole and +psi+ is the + # rotation about it (the prime meridian). Angles are in radians and rates, + # when present, in radians per day — matching ephem's per-day rate + # convention for SPK velocities (divide by 86400 for radians per second). + class Orientation + # @return [Numeric] first Euler angle (radians) + attr_reader :phi + + # @return [Numeric] second Euler angle (radians) + attr_reader :theta + + # @return [Numeric] third Euler angle (radians) + attr_reader :psi + + # @return [Array, nil] [phi, theta, psi] rates (radians/day), + # or nil when the orientation carries no rates + attr_reader :rates + + # @param phi [Numeric] first Euler angle (radians) + # @param theta [Numeric] second Euler angle (radians) + # @param psi [Numeric] third Euler angle (radians) + # @param rates [Array, nil] optional [phi, theta, psi] rates + # (radians/day) + # @raise [Ephem::InvalidInputError] if any angle or rate is not numeric + def initialize(phi, theta, psi, rates: nil) + unless phi.is_a?(Numeric) && theta.is_a?(Numeric) && psi.is_a?(Numeric) + raise InvalidInputError, "Orientation angles must be numeric" + end + + unless rates.nil? || valid_rates?(rates) + raise InvalidInputError, "Orientation rates must be three numerics" + end + + @phi = phi + @theta = theta + @psi = psi + @rates = rates&.freeze + freeze + end + + def self.[](phi, theta, psi, rates: nil) + new(phi, theta, psi, rates: rates) + end + + # @return [Boolean] whether this orientation carries rates + def rates? + !@rates.nil? + end + + # @return [Array] the three Euler angles [phi, theta, psi] + def to_a + [phi, theta, psi] + end + + # The rotation matrix that maps the reference frame into the body-fixed + # frame, built from the 3-1-3 (Z-X-Z) Euler angles: + # +M = Rz(psi) * Rx(theta) * Rz(phi)+. Rates are ignored. + # + # @return [Array>] a 3x3 rotation matrix + def to_matrix + Rotation.multiply( + Rotation.about_z(psi), + Rotation.about_x(theta), + Rotation.about_z(phi) + ) + end + + # @param index [Integer] 0 for phi, 1 for theta, 2 for psi + # @return [Numeric] the angle at the given index + # @raise [Ephem::IndexError] if index is not 0, 1, or 2 + def [](index) + case index + when 0 then phi + when 1 then theta + when 2 then psi + else raise IndexError, "Invalid index: #{index}" + end + end + + def inspect + body = "phi: #{phi}, theta: #{theta}, psi: #{psi}" + body += ", rates: #{rates}" if rates? + "Orientation[#{body}]" + end + alias_method :to_s, :inspect + + def hash + [phi, theta, psi, rates, self.class].hash + end + + def ==(other) + unless other.is_a?(self.class) + raise InvalidInputError, "Can only compare with another Orientation" + end + + to_a == other.to_a && rates == other.rates + end + alias_method :eql?, :== + + private + + def valid_rates?(rates) + rates.is_a?(Array) && + rates.size == 3 && + rates.all?(Numeric) + end + end + end +end diff --git a/lib/ephem/core/rotation.rb b/lib/ephem/core/rotation.rb new file mode 100644 index 0000000..32c24a3 --- /dev/null +++ b/lib/ephem/core/rotation.rb @@ -0,0 +1,82 @@ +# frozen_string_literal: true + +module Ephem + module Core + # Builds and applies 3x3 rotation matrices. Kernel-agnostic: callers choose + # the axis sequence and order that matches their frame convention. + # + # The elementary rotations use the coordinate-frame (passive) convention: + # they express a fixed vector in a frame rotated by +angle about the axis + module Rotation + # @param angle [Numeric] rotation angle in radians + # @return [Array>] rotation about the X axis + def self.about_x(angle) + cosine = Math.cos(angle) + sine = Math.sin(angle) + [ + [1.0, 0.0, 0.0], + [0.0, cosine, sine], + [0.0, -sine, cosine] + ] + end + + # @param angle [Numeric] rotation angle in radians + # @return [Array>] rotation about the Y axis + def self.about_y(angle) + cosine = Math.cos(angle) + sine = Math.sin(angle) + [ + [cosine, 0.0, -sine], + [0.0, 1.0, 0.0], + [sine, 0.0, cosine] + ] + end + + # @param angle [Numeric] rotation angle in radians + # @return [Array>] rotation about the Z axis + def self.about_z(angle) + cosine = Math.cos(angle) + sine = Math.sin(angle) + [ + [cosine, sine, 0.0], + [-sine, cosine, 0.0], + [0.0, 0.0, 1.0] + ] + end + + # Product of rotation matrices in the given order, as standard matrix + # multiplication: +multiply(a, b, c)+ returns +a * b * c+. + # + # @param matrices [Array>>] one or more 3x3 matrices + # @return [Array>] the combined rotation matrix + def self.multiply(*matrices) + matrices.reduce { |product, matrix| multiply_pair(product, matrix) } + end + + # Applies a rotation matrix to a vector. + # + # @param matrix [Array>] a 3x3 rotation matrix + # @param vector [Core::Vector, Array] the vector to rotate + # @return [Core::Vector] the rotated vector + def self.apply(matrix, vector) + x, y, z = vector.to_a + Vector.new( + matrix[0][0] * x + matrix[0][1] * y + matrix[0][2] * z, + matrix[1][0] * x + matrix[1][1] * y + matrix[1][2] * z, + matrix[2][0] * x + matrix[2][1] * y + matrix[2][2] * z + ) + end + + def self.multiply_pair(left, right) + Array.new(3) do |row| + Array.new(3) do |column| + left[row][0] * right[0][column] + + left[row][1] * right[1][column] + + left[row][2] * right[2][column] + end + end + end + private_class_method :multiply_pair + end + end +end diff --git a/lib/ephem/download.rb b/lib/ephem/download.rb index 32f0513..7af24d9 100644 --- a/lib/ephem/download.rb +++ b/lib/ephem/download.rb @@ -11,6 +11,8 @@ module Ephem class Download JPL_BASE_URL = "https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/" IMCCE_BASE_URL = "https://ftp.imcce.fr/pub/ephem/planets/" + NAIF_PCK_BASE_URL = + "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/" JPL_KERNELS = %w[ de102.bsp @@ -59,6 +61,11 @@ class Download de441.bsp ].freeze + NAIF_PCK_KERNELS = %w[ + moon_pa_de421_1900-2050.bpc + moon_pa_de440_200625.bpc + ].freeze + IMCCE_KERNELS = { "inpop10b.bsp" => "inpop10b_TDB_m100_p100_spice.bsp", "inpop10b_large.bsp" => "inpop10b_TDB_m1000_p1000_spice.bsp", @@ -90,7 +97,8 @@ class Download "inpop21a_large.bsp" => "inpop21a/inpop21a_TDB_m1000_p1000_spice.tar.gz" }.freeze - SUPPORTED_KERNELS = (JPL_KERNELS + IMCCE_KERNELS.keys).freeze + SUPPORTED_KERNELS = + (JPL_KERNELS + NAIF_PCK_KERNELS + IMCCE_KERNELS.keys).freeze def self.call(name:, target:) new(name, target).call @@ -104,12 +112,22 @@ def initialize(name, target_path) def call FileUtils.mkdir_p(@target_path.dirname) - jpl_kernel? ? download_jpl : download_imcce + download_kernel true end private + def download_kernel + if jpl_kernel? + download_jpl + elsif pck_kernel? + download_naif_pck + else + download_imcce + end + end + def validate_requested_kernel! unless SUPPORTED_KERNELS.include?(@name) raise UnsupportedError, @@ -121,8 +139,20 @@ def jpl_kernel? JPL_KERNELS.include?(@name) end + def pck_kernel? + NAIF_PCK_KERNELS.include?(@name) + end + def download_jpl - uri = URI.join(JPL_BASE_URL, @name) + download_direct(JPL_BASE_URL) + end + + def download_naif_pck + download_direct(NAIF_PCK_BASE_URL) + end + + def download_direct(base_url) + uri = URI.join(base_url, @name) @target_path.open("wb") do |file| stream_http_to_file(uri, file) end diff --git a/lib/ephem/excerpt.rb b/lib/ephem/excerpt.rb index 09e8971..b4bfe4a 100644 --- a/lib/ephem/excerpt.rb +++ b/lib/ephem/excerpt.rb @@ -1,9 +1,9 @@ # frozen_string_literal: true module Ephem - # The Excerpt class creates SPK file excerpts with reduced time spans and - # target bodies. This is useful for creating smaller files that focus only on - # the data needed for specific applications. + # The Excerpt class creates DAF excerpts (SPK or binary PCK) with reduced time + # spans and bodies. This is useful for creating smaller files that focus only + # on the data needed for specific applications. # # @example Create an excerpt with specific time range and bodies # spk = Ephem::SPK.open("de421.bsp") @@ -19,10 +19,10 @@ class Excerpt J2000_EPOCH = Core::Constants::Time::J2000_EPOCH RECORD_SIZE = 1024 - # @param spk [Ephem::SPK] The SPK object to create an excerpt from - def initialize(spk) - @spk = spk - @daf = spk.daf + # @param kernel [Ephem::SPK, Ephem::PCK] The kernel to excerpt from + def initialize(kernel) + @kernel = kernel + @daf = kernel.daf @binary_reader = @daf.binary_reader end @@ -31,11 +31,12 @@ def initialize(spk) # @param output_path [String] Path where the excerpt will be written # @param start_jd [Float] Start time as Julian Date # @param end_jd [Float] End time as Julian Date - # @param target_ids [Array, nil] Optional list of target IDs to - # include + # @param target_ids [Array, nil] Optional list of target/body IDs + # to include # @param debug [Boolean] Whether to print debug information # - # @return [Ephem::SPK] A new SPK instance for the excerpt file + # @return [Ephem::SPK, Ephem::PCK] A new instance for the excerpt file, + # matching the source kernel kind def extract(output_path:, start_jd:, end_jd:, target_ids: nil, debug: false) start_seconds = seconds_since_j2000(start_jd) end_seconds = seconds_since_j2000(end_jd) @@ -46,11 +47,15 @@ def extract(output_path:, start_jd:, end_jd:, target_ids: nil, debug: false) process_segments(writer, start_seconds, end_seconds, target_ids, debug) output_file.close - SPK.open(output_path) + reopen(output_path) end private + def reopen(path) + (@daf.file_type == :pck) ? PCK.open(path) : SPK.open(path) + end + def seconds_since_j2000(jd) (jd - J2000_EPOCH) * S_PER_DAY end diff --git a/lib/ephem/io/daf.rb b/lib/ephem/io/daf.rb index ff684e1..cbbce85 100644 --- a/lib/ephem/io/daf.rb +++ b/lib/ephem/io/daf.rb @@ -15,6 +15,11 @@ def comments @comments ||= load_comments end + # @return [Symbol, nil] + def file_type + @file_type ||= detect_file_type + end + def summaries(&block) return enum_for(:summaries) unless block_given? @@ -36,6 +41,18 @@ def close private + def detect_file_type + case @record_data.locator_identifier + when "DAF/SPK" then :spk + when "DAF/PCK" then :pck + else + case @record_data.integer_count + when 6 then :spk + when 5 then :pck + end + end + end + def setup_file_format file_record = @binary_reader.read_record(1) endianness_info = EndiannessManager.new(file_record).detect_endianness diff --git a/lib/ephem/pck.rb b/lib/ephem/pck.rb new file mode 100644 index 0000000..4e6da56 --- /dev/null +++ b/lib/ephem/pck.rb @@ -0,0 +1,118 @@ +# frozen_string_literal: true + +module Ephem + # Reads a binary PCK (+DAF/PCK+) orientation kernel: the orientation of one or + # more body frames over time, expressed as Euler angles. + class PCK + # Index of the data type within a PCK summary descriptor + # ([start, end, body, reference_frame, data_type, start_i, end_i]). + DATA_TYPE_IDENTIFIER = 4 + + attr_reader :daf, :segments + + # @param daf [Ephem::IO::DAF] DAF containing PCK data + # @raise [ArgumentError] if the DAF is nil + def initialize(daf:) + raise ArgumentError, "DAF cannot be nil" if daf.nil? + + @daf = daf + @segments = load_segments + @bodies = build_bodies + end + + # Opens a binary PCK file. + # + # @param path [String] Path to the PCK (+.bpc+) file + # @return [PCK] + # @raise [ArgumentError] if the file is not a binary PCK or cannot be read + def self.open(path) + daf = IO::DAF.new(File.open(path, "rb")) + unless daf.file_type == :pck + raise ArgumentError, "#{path} is not a binary PCK (DAF/PCK) file" + end + + new(daf: daf) + rescue Errno::EACCES => e + raise ArgumentError, "File permission denied: #{path} (#{e.message})" + rescue + daf&.close + raise + end + + # @return [void] + def close + @daf&.close + @segments&.each(&:clear_data) + end + + # Retrieves the orientation source for a body frame. + # + # @param body [Integer] NAIF frame ID of the oriented body + # @return [Segments::OrientationSegment, Segments::OrientationGroup] a + # single segment, or a group routing each query to the covering segment + # when the body spans several time intervals + # @raise [KeyError] if no segment is found for the given body + def [](body) + @bodies.fetch(body) do + raise KeyError, "No orientation segment found for body: #{body}" + end + end + + # @return [String] the comments stored in the PCK file + def comments + @daf.comments + end + + # @yieldparam segment [Segments::OrientationSegment] + # @return [Enumerator] if no block is given + def each_segment(&block) + return enum_for(:each_segment) unless block_given? + + @segments.each(&block) + end + + # @return [String] a description of the PCK file and its segments + def to_s + <<~DESCRIPTION + PCK file with #{@segments.size} segments: + #{@segments.join("\n")} + DESCRIPTION + end + + def excerpt(output_path:, start_jd:, end_jd:, target_ids: nil, debug: false) + Excerpt + .new(self) + .extract( + output_path: output_path, + start_jd: start_jd, + end_jd: end_jd, + target_ids: target_ids, + debug: debug + ) + end + + private + + def load_segments + @daf.summaries.map do |source, descriptor| + build_segment(source: source, descriptor: descriptor) + end + end + + def build_bodies + @segments.group_by(&:body).transform_values do |segments| + Segments::OrientationGroup.wrap(segments) + end + end + + def build_segment(source:, descriptor:) + data_type = descriptor[DATA_TYPE_IDENTIFIER] + segment_class = Segments::Registry.lookup(:pck, data_type) + unless segment_class + raise UnsupportedError, "Unsupported PCK data type: #{data_type}" + end + + segment_class.new(daf: @daf, source: source, descriptor: descriptor) + end + end +end diff --git a/lib/ephem/segments/base_segment.rb b/lib/ephem/segments/base_segment.rb index a94442d..f47d3e8 100644 --- a/lib/ephem/segments/base_segment.rb +++ b/lib/ephem/segments/base_segment.rb @@ -38,6 +38,10 @@ class BaseSegment attr_reader :center # @return [String] the source of the segment attr_reader :source + # @return [Float] start of coverage as a Julian Date + attr_reader :start_jd + # @return [Float] end of coverage as a Julian Date + attr_reader :end_jd # Initialize a new segment # @@ -55,14 +59,7 @@ class BaseSegment def initialize(daf:, source:, descriptor:) @daf = daf @source = source - @start_second, - @end_second, - @target, - @center, - @frame, - @data_type, - @start_i, - @end_i = descriptor + parse_descriptor(descriptor) @start_jd = compute_julian_date(@start_second) @end_jd = compute_julian_date(@end_second) end @@ -97,6 +94,15 @@ def describe(verbose: false) DESCRIPTION end + # Whether the given time falls within this segment's coverage. + # + # @param tdb [Numeric] time in TDB Julian Date + # @param tdb2 [Numeric] optional fractional part of the TDB date + # @return [Boolean] + def covers?(tdb, tdb2 = 0.0) + (tdb + tdb2).between?(@start_jd, @end_jd) + end + def compute(_tdb, _tdb2 = 0.0) raise NotImplementedError, "#{self.class} has not implemented compute() for data type #{@data_type}" @@ -113,6 +119,17 @@ def clear_data private + def parse_descriptor(descriptor) + @start_second, + @end_second, + @target, + @center, + @frame, + @data_type, + @start_i, + @end_i = descriptor + end + def compute_julian_date(seconds) Time::J2000_EPOCH + seconds / Time::SECONDS_PER_DAY end diff --git a/lib/ephem/segments/chebyshev_type2.rb b/lib/ephem/segments/chebyshev_type2.rb new file mode 100644 index 0000000..79d89c1 --- /dev/null +++ b/lib/ephem/segments/chebyshev_type2.rb @@ -0,0 +1,142 @@ +# frozen_string_literal: true + +module Ephem + module Segments + # Shared evaluation machinery for DAF "type 2" Chebyshev segments. + module ChebyshevType2 + include Core::Constants + + # @return [void] + def clear_data + @data_lock.synchronize do + @data_loaded = false + @midpoints = nil + @radii = nil + @coefficients = nil + end + end + + private + + def load_data + @data_lock.synchronize do + return if @data_loaded + + process_coefficient_data(load_coefficient_data) + + @data_loaded = true + end + end + + def load_coefficient_data + metadata = @daf.read_array(@end_i - 3, @end_i) + _start_index, _end_index, record_size, segment_count = metadata + + coefficient_count = ((record_size - 2) / component_count).to_i + coefficients_raw = @daf.map_array(@start_i, @end_i - 4) + + [ + coefficients_raw, + record_size.to_i, + segment_count.to_i, + coefficient_count + ] + end + + def process_coefficient_data(data) + coefficients_raw, record_size, segment_count, coefficient_count = data + + coefficients = coefficients_raw.each_slice(record_size).to_a + + @midpoints = coefficients.map { |row| row[0] } + @radii = coefficients.map { |row| row[1] } + n_terms = coefficient_count + n_components = component_count + + @coefficients = Array.new(segment_count) do |i| + row = coefficients[i][2..] + Array.new(n_terms) do |k| + Array.new(n_components) do |j| + row[k + j * n_terms] + end + end + end + end + + def convert_to_seconds(tdb, tdb2) + case tdb + when Array + tdb.map { |t| time_to_seconds(t, tdb2) } + else + time_to_seconds(tdb, tdb2) + end + end + + def time_to_seconds(time, offset) + (time - Time::J2000_EPOCH + offset) * Time::SECONDS_PER_DAY + end + + def generate_position(tdb_seconds) + interval = find_interval(tdb_seconds) + normalized_time = compute_normalized_time(tdb_seconds, interval) + coeffs = @coefficients[interval] + Computation::ChebyshevPolynomial.evaluate(coeffs, normalized_time) + end + + def generate_single(tdb_seconds) + interval = find_interval(tdb_seconds) + normalized_time = compute_normalized_time(tdb_seconds, interval) + + coeffs = @coefficients[interval] # already [n_terms][3] + Computation::ChebyshevPolynomial.evaluate_with_derivative( + coeffs, + normalized_time, + @radii[interval] + ) + end + + def generate_multiple(tdb_seconds) + tdb_seconds.map { |time| generate_single(time) } + end + + def find_interval(tdb_seconds) + left = 0 + right = @midpoints.size - 1 + + if @last_interval && time_in_interval?(tdb_seconds, @last_interval) + return @last_interval + end + + while left <= right + mid = (left + right) / 2 + min_time = @midpoints[mid] - @radii[mid] + max_time = @midpoints[mid] + @radii[mid] + + if tdb_seconds < min_time + right = mid - 1 + elsif tdb_seconds > max_time + left = mid + 1 + else + @last_interval = mid + return mid + end + end + + raise OutOfRangeError.new( + "Time #{tdb_seconds} is outside the coverage of this segment", + tdb_seconds + ) + end + + def time_in_interval?(time, interval) + min_time = @midpoints[interval] - @radii[interval] + max_time = @midpoints[interval] + @radii[interval] + time.between?(min_time, max_time) + end + + def compute_normalized_time(time_seconds, interval) + (time_seconds - @midpoints[interval]) / @radii[interval] + end + end + end +end diff --git a/lib/ephem/segments/orientation_group.rb b/lib/ephem/segments/orientation_group.rb new file mode 100644 index 0000000..bb47475 --- /dev/null +++ b/lib/ephem/segments/orientation_group.rb @@ -0,0 +1,59 @@ +# frozen_string_literal: true + +module Ephem + module Segments + # The orientation segments for one PCK body frame. Routes each query to the + # segment covering the requested time. Returned by PCK#[]. + # + # @see Ephem::Segments::OrientationSegment + class OrientationGroup < SegmentGroup + include OrientationSource + + # @return [Integer] NAIF frame ID of the oriented body frame + def body + @segments.first.body + end + + # @return [Integer] NAIF ID of the inertial reference frame + def reference_frame + @segments.first.reference_frame + end + + # Euler angles at the given time, without rates. See + # {OrientationSegment#angles_at}. + # + # @param tdb [Numeric, Array] Time(s) in TDB Julian Date + # @param tdb2 [Numeric] Optional fractional part of TDB date + # @return [Core::Orientation, Array] + def angles_at(tdb, tdb2 = 0.0) + query(tdb, tdb2) do |segment, time, fraction| + segment.angles_at(time, fraction) + end + end + + # Euler angles and their rates at the given time. See + # {OrientationSegment#orientation_at}. + # + # @param tdb [Numeric, Array] Time(s) in TDB Julian Date + # @param tdb2 [Numeric] Optional fractional part of TDB date + # @return [Core::Orientation, Array] + def orientation_at(tdb, tdb2 = 0.0) + query(tdb, tdb2) do |segment, time, fraction| + segment.orientation_at(time, fraction) + end + end + + # The reference-frame to body-fixed rotation matrix at the given time. + # See {OrientationSegment#matrix_at}. + # + # @param tdb [Numeric, Array] Time(s) in TDB Julian Date + # @param tdb2 [Numeric] Optional fractional part of TDB date + # @return [Array>, Array>>] + def matrix_at(tdb, tdb2 = 0.0) + query(tdb, tdb2) do |segment, time, fraction| + segment.matrix_at(time, fraction) + end + end + end + end +end diff --git a/lib/ephem/segments/orientation_segment.rb b/lib/ephem/segments/orientation_segment.rb new file mode 100644 index 0000000..9321956 --- /dev/null +++ b/lib/ephem/segments/orientation_segment.rb @@ -0,0 +1,118 @@ +# frozen_string_literal: true + +module Ephem + module Segments + # Binary PCK orientation segment (data type 2): the orientation of a body + # frame relative to an inertial reference frame, stored as three Euler + # angles in Chebyshev coefficients. + class OrientationSegment < BaseSegment + include ChebyshevType2 + include OrientationSource + + COMPONENT_COUNT = 3 # phi, theta, psi + + def initialize(daf:, source:, descriptor:) + super + @data_loaded = false + @data_lock = Mutex.new + end + + # @return [Integer] NAIF frame ID of the oriented body frame + alias_method :body, :target + # @return [Integer] NAIF ID of the inertial reference frame + alias_method :reference_frame, :center + + # Euler angles at the given time, without rates. + # + # @param tdb [Numeric, Array] Time(s) in TDB Julian Date + # @param tdb2 [Numeric] Optional fractional part of TDB date + # @return [Core::Orientation, Array] angles in radians + # @raise [Ephem::OutOfRangeError] if time is outside segment coverage + def angles_at(tdb, tdb2 = 0.0) + load_data + tdb_seconds = convert_to_seconds(tdb, tdb2) + + case tdb_seconds + when Numeric + to_orientation(generate_position(tdb_seconds)) + else + tdb_seconds.map do |seconds| + to_orientation(generate_position(seconds)) + end + end + end + + # Euler angles and their rates at the given time. + # + # @param tdb [Numeric, Array] Time(s) in TDB Julian Date + # @param tdb2 [Numeric] Optional fractional part of TDB date + # @return [Core::Orientation, Array] angles (radians) + # carrying rates (radians/day) + # @raise [Ephem::OutOfRangeError] if time is outside segment coverage + def orientation_at(tdb, tdb2 = 0.0) + load_data + tdb_seconds = convert_to_seconds(tdb, tdb2) + + case tdb_seconds + when Numeric + to_orientation(*generate_single(tdb_seconds)) + else + generate_multiple(tdb_seconds).map do |angles, rates| + to_orientation(angles, rates) + end + end + end + + # The reference-frame to body-fixed rotation matrix at the given time, + # built from the 3-1-3 Euler angles. See {Core::Orientation#to_matrix}. + # + # @param tdb [Numeric, Array] Time(s) in TDB Julian Date + # @param tdb2 [Numeric] Optional fractional part of TDB date + # @return [Array>, Array>>] a 3x3 matrix, + # or one per time for an array input + # @raise [Ephem::OutOfRangeError] if time is outside segment coverage + def matrix_at(tdb, tdb2 = 0.0) + angles = angles_at(tdb, tdb2) + angles.is_a?(Array) ? angles.map(&:to_matrix) : angles.to_matrix + end + + def describe(verbose: false) + start_date = format_date(*julian_to_gregorian(@start_jd)) + end_date = format_date(*julian_to_gregorian(@end_jd)) + + description = + "#{start_date}..#{end_date} Type #{@data_type} orientation of " \ + "frame #{body} relative to frame #{reference_frame}" + return description unless verbose + + <<~DESCRIPTION.chomp + #{description} + source=#{@source} + DESCRIPTION + end + + private + + def parse_descriptor(descriptor) + @start_second, + @end_second, + @target, + @center, + @data_type, + @start_i, + @end_i = descriptor + @frame = @center + end + + def component_count + COMPONENT_COUNT + end + + def to_orientation(angles, rates = nil) + Core::Orientation.new(angles[0], angles[1], angles[2], rates: rates) + end + + Registry.register(:pck, 2, self) + end + end +end diff --git a/lib/ephem/segments/orientation_source.rb b/lib/ephem/segments/orientation_source.rb new file mode 100644 index 0000000..613a1e5 --- /dev/null +++ b/lib/ephem/segments/orientation_source.rb @@ -0,0 +1,17 @@ +# frozen_string_literal: true + +module Ephem + module Segments + module OrientationSource + def compute(*) + raise NotImplementedError, + "Use #angles_at or #orientation_at for orientation kernels" + end + + def compute_and_differentiate(*) + raise NotImplementedError, + "Use #orientation_at for orientation kernels" + end + end + end +end diff --git a/lib/ephem/segments/position_group.rb b/lib/ephem/segments/position_group.rb new file mode 100644 index 0000000..fd7d29d --- /dev/null +++ b/lib/ephem/segments/position_group.rb @@ -0,0 +1,46 @@ +# frozen_string_literal: true + +module Ephem + module Segments + # The position segments for one SPK center/target pair. Routes each query to + # the segment covering the requested time. Returned by SPK#[]. + # + # @see Ephem::Segments::Segment + class PositionGroup < SegmentGroup + # @return [Integer] the center body ID + def center + @segments.first.center + end + + # @return [Integer] the target body ID + def target + @segments.first.target + end + + # Position at the given time. See {Segment#compute}. + # + # @param tdb [Numeric, Array] Time(s) in TDB Julian Date + # @param tdb2 [Numeric] Optional fractional part of TDB date + # @return [Core::Vector, Array] + def compute(tdb, tdb2 = 0.0) + query(tdb, tdb2) do |segment, time, fraction| + segment.compute(time, fraction) + end + end + alias_method :position_at, :compute + + # Position and velocity at the given time. See + # {Segment#compute_and_differentiate}. + # + # @param tdb [Numeric, Array] Time(s) in TDB Julian Date + # @param tdb2 [Numeric] Optional fractional part of TDB date + # @return [Core::State, Array] + def compute_and_differentiate(tdb, tdb2 = 0.0) + query(tdb, tdb2) do |segment, time, fraction| + segment.compute_and_differentiate(time, fraction) + end + end + alias_method :state_at, :compute_and_differentiate + end + end +end diff --git a/lib/ephem/segments/registry.rb b/lib/ephem/segments/registry.rb index ff21304..21d6e29 100644 --- a/lib/ephem/segments/registry.rb +++ b/lib/ephem/segments/registry.rb @@ -2,9 +2,15 @@ module Ephem module Segments - class Registry - def self.register(type, klass) - SPK::SEGMENT_CLASSES[type] = klass + module Registry + TABLES = {spk: {}, pck: {}}.freeze + + def self.register(kind, type, klass) + TABLES.fetch(kind)[type] = klass + end + + def self.lookup(kind, type, default = nil) + TABLES.fetch(kind).fetch(type, default) end end end diff --git a/lib/ephem/segments/segment.rb b/lib/ephem/segments/segment.rb index ff1db3b..c629d7b 100644 --- a/lib/ephem/segments/segment.rb +++ b/lib/ephem/segments/segment.rb @@ -2,18 +2,9 @@ module Ephem module Segments - # Manages data segments within SPICE kernel (SPK) files, providing methods - # to compute positions and velocities of celestial bodies using Chebyshev - # polynomial approximations. - # - # Each segment contains data for a specific celestial body (target) relative - # to another body (center) within a specific time range. The data is stored - # as Chebyshev polynomial coefficients that can be evaluated to obtain - # position and velocity vectors. - # - # The class provides thread-safe data loading and caching mechanisms to - # optimize performance while ensuring data consistency in multithreaded - # environments. + # SPK trajectory segment: position (type 2) and position/velocity (type 3) + # of a target body relative to a center body, stored as Chebyshev + # coefficients. # # @example Computing position at a specific time # segment = Ephem::Segments::Segment.new( @@ -25,29 +16,18 @@ module Segments # # @example Computing position and velocity # state = segment.compute_and_differentiate(time) # returns State - # position = state.position # Vector - # velocity = state.velocity # Vector # # @see Ephem::Core::Vector # @see Ephem::Core::State - # @see Ephem::Computation::ChebyshevPolynomial + # @see Ephem::Segments::ChebyshevType2 class Segment < BaseSegment + include ChebyshevType2 + COMPONENT_COUNTS = { 2 => 3, # Type 2: position (x, y, z) 3 => 6 # Type 3: position (x, y, z) and velocity (vx, vy, vz) }.freeze - # @param daf [Ephem::IO::DAF] DAF file object containing the segment data - # @param source [String] Name of the source SPK file - # @param descriptor [Array] Array containing segment metadata: - # - start_second [Float] Start time in seconds from J2000 - # - end_second [Float] End time in seconds from J2000 - # - target [Integer] NAIF ID of target body - # - center [Integer] NAIF ID of center body - # - frame [Integer] Reference frame ID - # - data_type [Integer] Type of data (2 for position, 3 for pos/vel) - # - start_i [Integer] Start index in DAF array - # - end_i [Integer] End index in DAF array def initialize(daf:, source:, descriptor:) super @data_loaded = false @@ -57,17 +37,10 @@ def initialize(daf:, source:, descriptor:) # Computes the position of the target body relative to the center body at # the specified time. # - # Uses Chebyshev polynomial approximation to interpolate the position from - # stored coefficients. The computation is thread-safe and uses cached data - # when available. - # # @param tdb [Numeric, Array] Time(s) in TDB Julian Date # @param tdb2 [Numeric] Optional fractional part of TDB date # @return [Ephem::Core::Vector] Position vector in kilometers # @raise [Ephem::OutOfRangeError] if time is outside segment coverage - # - # @example Computing Earth's position relative to Solar System Barycenter - # position = segment.compute(2451545.0) # J2000 epoch def compute(tdb, tdb2 = 0.0) load_data tdb_seconds = convert_to_seconds(tdb, tdb2) @@ -87,21 +60,12 @@ def compute(tdb, tdb2 = 0.0) # Computes both position and velocity vectors at the specified time. # - # Uses Chebyshev polynomial approximation and its derivative to compute - # both position and velocity. The computation is thread-safe and uses - # cached data when available. - # # @param tdb [Numeric, Array] Time(s) in TDB Julian Date # @param tdb2 [Numeric] Optional fractional part of TDB date # @return [Ephem::Core::State, Array] State object(s) - # containing position and velocity vectors. Returns array if input is - # array. + # containing position (km) and velocity (km/day) vectors. Returns an + # array if the input is an array. # @raise [Ephem::OutOfRangeError] if time is outside segment coverage - # - # @example Computing Earth's state vector - # state = segment.compute_and_differentiate(2451545.0) - # position = state.position # in kilometers - # velocity = state.velocity # in kilometers/second def compute_and_differentiate(tdb, tdb2 = 0.0) load_data tdb_seconds = convert_to_seconds(tdb, tdb2) @@ -124,186 +88,16 @@ def compute_and_differentiate(tdb, tdb2 = 0.0) end alias_method :state_at, :compute_and_differentiate - # Clears cached coefficient data, forcing reload on next computation. - # - # This method is thread-safe and can be used to free memory or force - # fresh data loading if needed. - # - # @return [void] - def clear_data - @data_lock.synchronize do - @data_loaded = false - @midpoints = nil - @radii = nil - @coefficients = nil - end - end - private - def load_data - # Synchronize access to data loading using a mutex lock - # to prevent race conditions in multithreaded environments - @data_lock.synchronize do - return if @data_loaded - - component_count = determine_component_count - coefficients_data = load_coefficient_data - process_coefficient_data(coefficients_data, component_count) - - @data_loaded = true - end - end - - def determine_component_count + def component_count COMPONENT_COUNTS.fetch(@data_type) do raise "Unsupported data type: #{@data_type}" end end - def load_coefficient_data - # Read metadata from the end of the segment - # start_index: index of first coefficient in segment - # end_index: index of last coefficient in segment - # record_size: total size of each record (coefficients + 2) - # segment_count: number of records in the segment - metadata = @daf.read_array(@end_i - 3, @end_i) - _start_index, _end_index, record_size, segment_count = metadata - - coefficient_count = ((record_size - 2) / determine_component_count).to_i - coefficients_raw = @daf.map_array(@start_i, @end_i - 4) - - [ - coefficients_raw, - record_size.to_i, - segment_count.to_i, - coefficient_count - ] - end - - def process_coefficient_data(data, component_count) - coefficients_raw, record_size, segment_count, coefficient_count = data - - coefficients = coefficients_raw.each_slice(record_size).to_a - - @midpoints = coefficients.map { |row| row[0] } - @radii = coefficients.map { |row| row[1] } - n_terms = coefficient_count - n_components = component_count - - @coefficients = Array.new(segment_count) do |i| - row = coefficients[i][2..] - Array.new(n_terms) do |k| - Array.new(n_components) do |j| - row[k + j * n_terms] - end - end - end - end - - def convert_to_seconds(tdb, tdb2) - case tdb - when Array - tdb.map { |t| time_to_seconds(t, tdb2) } - else - time_to_seconds(tdb, tdb2) - end - end - - def time_to_seconds(time, offset) - (time - Time::J2000_EPOCH + offset) * Time::SECONDS_PER_DAY - end - - def generate(tdb, tdb2) - load_data - tdb_seconds = convert_to_seconds(tdb, tdb2) - - case tdb_seconds - when Numeric - generate_single(tdb_seconds) - else - generate_multiple(tdb_seconds) - end - end - - def generate_position(tdb_seconds) - interval = find_interval(tdb_seconds) - normalized_time = compute_normalized_time(tdb_seconds, interval) - coeffs = @coefficients[interval] - Computation::ChebyshevPolynomial.evaluate(coeffs, normalized_time) - end - - def generate_single(tdb_seconds) - interval = find_interval(tdb_seconds) - normalized_time = compute_normalized_time(tdb_seconds, interval) - - coeffs = @coefficients[interval] # already [n_terms][3] - position = Computation::ChebyshevPolynomial.evaluate( - coeffs, - normalized_time - ) - velocity = Computation::ChebyshevPolynomial.evaluate_derivative( - coeffs, - normalized_time, - @radii[interval] - ) - [position, velocity] - end - - def generate_multiple(tdb_seconds) - positions = [] - velocities = [] - - tdb_seconds.each do |time| - pos, vel = generate_single(time) - positions << pos - velocities << vel - end - - [positions, velocities] - end - - def find_interval(tdb_seconds) - left = 0 - right = @midpoints.size - 1 - - if @last_interval && time_in_interval?(tdb_seconds, @last_interval) - return @last_interval - end - - while left <= right - mid = (left + right) / 2 - min_time = @midpoints[mid] - @radii[mid] - max_time = @midpoints[mid] + @radii[mid] - - if tdb_seconds < min_time - right = mid - 1 - elsif tdb_seconds > max_time - left = mid + 1 - else - @last_interval = mid - return mid - end - end - - raise OutOfRangeError.new( - "Time #{tdb_seconds} is outside the coverage of this segment", - tdb_seconds - ) - end - - def time_in_interval?(time, interval) - min_time = @midpoints[interval] - @radii[interval] - max_time = @midpoints[interval] + @radii[interval] - time.between?(min_time, max_time) - end - - def compute_normalized_time(time_seconds, interval) - (time_seconds - @midpoints[interval]) / @radii[interval] - end - - Registry.register(2, self) - Registry.register(3, self) + Registry.register(:spk, 2, self) + Registry.register(:spk, 3, self) end end end diff --git a/lib/ephem/segments/segment_group.rb b/lib/ephem/segments/segment_group.rb new file mode 100644 index 0000000..98a38fa --- /dev/null +++ b/lib/ephem/segments/segment_group.rb @@ -0,0 +1,84 @@ +# frozen_string_literal: true + +module Ephem + module Segments + # Several segments that share the same key (an SPK center/target pair, or a + # PCK body) but cover different, contiguous time intervals. Each query is + # routed to the segment covering the requested time, so a body that a kernel + # splits across several intervals behaves as a single, continuous source. + # + # SPK and PCK only build a group when a key actually has more than one + # segment; the common single-segment case returns the bare segment, so this + # routing never sits in the hot path for it. + # + # Subclasses ({PositionGroup}, {OrientationGroup}) add the query methods + # appropriate to the segments they hold. + class SegmentGroup + # Wraps segments that share a key. A single segment is returned as-is, so + # the common case carries no routing overhead; only a key spanning several + # time intervals becomes a group. + # + # @param segments [Array] segments sharing the same key + # @return [BaseSegment, SegmentGroup] + def self.wrap(segments) + segments.one? ? segments.first : new(segments) + end + + # @return [Array] the underlying segments + attr_reader :segments + + # @param segments [Array] segments sharing the same key + def initialize(segments) + @segments = segments + end + + # Clears cached data for every segment in the group. + # + # @return [void] + def clear_data + @segments.each(&:clear_data) + end + + def to_s + @segments.join("\n") + end + + private + + # Routes a query to the covering segment(s) and assembles the result. For + # a scalar time the block is called once with that segment and time; for + # an array, times are grouped by covering segment so each is queried in a + # single batched call, then results are reassembled in input order. + def query(tdb, tdb2) + if tdb.is_a?(Array) + query_many(tdb, tdb2) { |segment, times| yield segment, times, tdb2 } + else + yield segment_for(tdb, tdb2), tdb, tdb2 + end + end + + def query_many(times, tdb2) + results = Array.new(times.size) + indices_by_segment = times.each_index.group_by do |index| + segment_for(times[index], tdb2) + end + + indices_by_segment.each do |segment, indices| + segment_results = yield(segment, indices.map { |index| times[index] }) + indices.each_with_index do |original_index, position| + results[original_index] = segment_results[position] + end + end + + results + end + + def segment_for(tdb, tdb2) + @segments.find { |segment| segment.covers?(tdb, tdb2) } || + raise(OutOfRangeError.new( + "Time #{tdb} is outside the coverage of this group", tdb + )) + end + end + end +end diff --git a/lib/ephem/spk.rb b/lib/ephem/spk.rb index 8a89d46..72cdb86 100644 --- a/lib/ephem/spk.rb +++ b/lib/ephem/spk.rb @@ -25,7 +25,6 @@ class SPK DE_FILENAME = "NIO2SPK" DATA_TYPE_IDENTIFIER = 5 - SEGMENT_CLASSES = {} attr_reader :daf, :segments, :pairs @@ -49,9 +48,16 @@ def initialize(daf:) # @raise [ArgumentError] If the file cannot be accessed due to permissions def self.open(path) daf = IO::DAF.new(File.open(path, "rb")) + if daf.file_type == :pck + raise ArgumentError, "#{path} is a binary PCK file, use Ephem::PCK.open" + end + new(daf: daf) rescue Errno::EACCES => e raise ArgumentError, "File permission denied: #{path} (#{e.message})" + rescue + daf&.close + raise end # Closes the SPK file and cleans up resources. @@ -77,8 +83,9 @@ def to_s # # @param center [Integer] NAIF ID of the center body # @param target [Integer] NAIF ID of the target body - # @return [Segments::BaseSegment] The segment containing data for the - # specified bodies + # @return [Segments::Segment, Segments::PositionGroup] the position source + # for the pair: a single segment, or a group routing each query to the + # covering segment when the pair spans several time intervals # @raise [KeyError] If no segment is found for the given center-target pair def [](center, target) @pairs.fetch([center, target]) do @@ -139,14 +146,18 @@ def load_segments end def build_pairs - @segments.to_h do |segment| - [[segment.center, segment.target], segment] - end + @segments + .group_by { |segment| [segment.center, segment.target] } + .transform_values { |segments| Segments::PositionGroup.wrap(segments) } end def build_segment(source:, descriptor:) data_type = descriptor[DATA_TYPE_IDENTIFIER] - segment_class = SEGMENT_CLASSES.fetch(data_type, Segments::BaseSegment) + segment_class = Segments::Registry.lookup( + :spk, + data_type, + Segments::BaseSegment + ) segment_class.new(daf: @daf, source: source, descriptor: descriptor) end end diff --git a/lib/ephem/tasks/validate_accuracy.rb b/lib/ephem/tasks/validate_accuracy.rb index d137864..af25702 100644 --- a/lib/ephem/tasks/validate_accuracy.rb +++ b/lib/ephem/tasks/validate_accuracy.rb @@ -128,7 +128,7 @@ def max_errors_output pos = @max_errors.slice(:dx, :dy, :dz) vel = @max_errors.slice(:dvx, :dvy, :dvz) "Max position errors (km): #{pos.map { |k, v| "#{k}=#{v}" }.join(", ")}\n" \ - "Max velocity errors (km/s): #{vel.map { |k, v| "#{k}=#{v}" }.join(", ")}" + "Max velocity errors (km/day): #{vel.map { |k, v| "#{k}=#{v}" }.join(", ")}" end end end diff --git a/spec/ephem/cli_spec.rb b/spec/ephem/cli_spec.rb index 3a5cb48..64d46eb 100644 --- a/spec/ephem/cli_spec.rb +++ b/spec/ephem/cli_spec.rb @@ -43,7 +43,7 @@ expect(output).to include("Ruby Ephem") expect(output).to include("excerpt") expect(output).to include("--targets") - expect(output).to include("Example:") + expect(output).to include("Examples:") end end @@ -185,6 +185,63 @@ FileUtils.remove_entry(temp_dir) if Dir.exist?(temp_dir) end end + + it "excerpts a binary PCK kernel, detecting the format" do + temp_dir = Dir.mktmpdir("ephem_test_") + input_file = File.join(temp_dir, "input.bpc") + output_file = File.join(temp_dir, "output.bpc") + + begin + FileUtils.cp(moon_pa_de440_excerpt, input_file) + + args = [ + "--targets", + "31008", + "2001-01-01", + "2002-01-01", + input_file, + output_file + ] + output = capture_stdout { Ephem::CLI.handle_excerpt(args) } + + expect(output).to include("Excerpt created successfully") + + excerpt = Ephem::PCK.open(output_file) + expect(excerpt[31008].body).to eq(31008) + excerpt.close + ensure + FileUtils.remove_entry(temp_dir) if Dir.exist?(temp_dir) + end + end + end + end + + describe ".open_kernel" do + it "opens a binary PCK file as Ephem::PCK" do + kernel = Ephem::CLI.open_kernel(moon_pa_de440_excerpt) + + expect(kernel).to be_a(Ephem::PCK) + + kernel.close + end + + it "opens an SPK file as Ephem::SPK" do + kernel = Ephem::CLI.open_kernel(test_spk) + + expect(kernel).to be_a(Ephem::SPK) + + kernel.close + end + + it "closes the DAF when the kernel cannot be built" do + daf = instance_double(Ephem::IO::DAF, file_type: :pck) + allow(File).to receive(:open).and_return(instance_double(File)) + allow(Ephem::IO::DAF).to receive(:new).and_return(daf) + allow(Ephem::PCK).to receive(:new).and_raise(Ephem::UnsupportedError) + + expect(daf).to receive(:close) + expect { Ephem::CLI.open_kernel("kernel.bpc") } + .to raise_error(Ephem::UnsupportedError) end end diff --git a/spec/ephem/computation/chebyshev_polynomial_spec.rb b/spec/ephem/computation/chebyshev_polynomial_spec.rb index aa177b0..14b523f 100644 --- a/spec/ephem/computation/chebyshev_polynomial_spec.rb +++ b/spec/ephem/computation/chebyshev_polynomial_spec.rb @@ -191,6 +191,35 @@ end end + describe ".evaluate_with_derivative" do + it "matches #evaluate and #evaluate_derivative exactly" do + coeffs = [ + [1.0, 2.0, 3.0], + [0.5, 0.2, -1.5], + [-0.25, 1.0, 2.0], + [0.1, -0.4, 0.7] + ] + radius = 2000.0 + + [-1.0, -0.3, 0.0, 0.42, 1.0].each do |t| + position, velocity = + described_class.evaluate_with_derivative(coeffs, t, radius) + + expect(position).to eq(described_class.evaluate(coeffs, t)) + expect(velocity) + .to eq(described_class.evaluate_derivative(coeffs, t, radius)) + end + end + + it "returns a zero velocity for a constant polynomial" do + position, velocity = + described_class.evaluate_with_derivative([[1.0, 2.0, 3.0]], 0.3, 100.0) + + expect(position).to eq([1.0, 2.0, 3.0]) + expect(velocity).to eq([0.0, 0.0, 0.0]) + end + end + def expect_vector_close(result, expected, tol = 1e-10) expect(result.size).to eq(expected.size) result.zip(expected).each do |a, b| diff --git a/spec/ephem/core/orientation_spec.rb b/spec/ephem/core/orientation_spec.rb new file mode 100644 index 0000000..36b34b8 --- /dev/null +++ b/spec/ephem/core/orientation_spec.rb @@ -0,0 +1,128 @@ +# frozen_string_literal: true + +RSpec.describe Ephem::Core::Orientation do + describe ".new" do + it "exposes the three Euler angles" do + orientation = described_class.new(1.0, 2.0, 3.0) + + expect(orientation.phi).to eq(1.0) + expect(orientation.theta).to eq(2.0) + expect(orientation.psi).to eq(3.0) + end + + it "carries no rates by default" do + orientation = described_class.new(1.0, 2.0, 3.0) + + expect(orientation.rates).to be_nil + expect(orientation.rates?).to be(false) + end + + it "carries rates when given" do + orientation = described_class.new(1.0, 2.0, 3.0, rates: [0.1, 0.2, 0.3]) + + expect(orientation.rates).to eq([0.1, 0.2, 0.3]) + expect(orientation.rates?).to be(true) + end + + it "is frozen" do + expect(described_class.new(1.0, 2.0, 3.0)).to be_frozen + end + + it "raises for non-numeric angles" do + expect { described_class.new(1.0, "x", 3.0) }.to raise_error( + Ephem::InvalidInputError, "Orientation angles must be numeric" + ) + end + + it "raises for malformed rates" do + expect { described_class.new(1.0, 2.0, 3.0, rates: [0.1, 0.2]) } + .to raise_error( + Ephem::InvalidInputError, "Orientation rates must be three numerics" + ) + end + end + + describe "#[] and #to_a" do + it "indexes angles and converts to an array" do + orientation = described_class.new(1.0, 2.0, 3.0) + + expect(orientation.to_a).to eq([1.0, 2.0, 3.0]) + expect(orientation[0]).to eq(1.0) + expect(orientation[2]).to eq(3.0) + end + + it "raises for an invalid index" do + expect { described_class.new(1.0, 2.0, 3.0)[3] } + .to raise_error(Ephem::IndexError) + end + end + + describe "#to_matrix" do + it "returns the identity for zero angles" do + matrix = described_class.new(0.0, 0.0, 0.0).to_matrix + + expect(matrix).to eq( + [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]] + ) + end + + it "composes the 3-1-3 sequence Rz(psi) * Rx(theta) * Rz(phi)" do + phi = 0.3 + theta = 0.7 + psi = 1.1 + expected = Ephem::Core::Rotation.multiply( + Ephem::Core::Rotation.about_z(psi), + Ephem::Core::Rotation.about_x(theta), + Ephem::Core::Rotation.about_z(phi) + ) + + expect(described_class.new(phi, theta, psi).to_matrix).to eq(expected) + end + + it "is orthonormal" do + matrix = described_class.new(1.2, 0.4, 5.9).to_matrix + product = Ephem::Core::Rotation.multiply(matrix, matrix.transpose) + + product.each_with_index do |row, row_index| + row.each_with_index do |value, column_index| + expected = (row_index == column_index) ? 1.0 : 0.0 + expect(value).to be_within(1e-12).of(expected) + end + end + end + end + + describe "#==" do + it "compares angles and rates" do + with_rates = described_class.new(1.0, 2.0, 3.0, rates: [0.1, 0.2, 0.3]) + same_with_rates = described_class.new( + 1.0, + 2.0, + 3.0, + rates: [0.1, 0.2, 0.3] + ) + without_rates = described_class.new(1.0, 2.0, 3.0) + + expect(with_rates).to eq(same_with_rates) + expect(with_rates).not_to eq(without_rates) + end + + it "raises when compared with a non-Orientation" do + expect { described_class.new(1.0, 2.0, 3.0) == 5 } + .to raise_error(Ephem::InvalidInputError) + end + end + + describe "#inspect" do + it "renders angles, and rates when present" do + angles_only = described_class.new(1.0, 2.0, 3.0) + with_rates = described_class.new(1.0, 2.0, 3.0, rates: [0.1, 0.2, 0.3]) + + expect(angles_only.inspect) + .to eq("Orientation[phi: 1.0, theta: 2.0, psi: 3.0]") + expect(with_rates.inspect).to eq( + "Orientation[phi: 1.0, theta: 2.0, psi: 3.0, rates: [0.1, 0.2, 0.3]]" + ) + end + end +end diff --git a/spec/ephem/core/rotation_spec.rb b/spec/ephem/core/rotation_spec.rb new file mode 100644 index 0000000..c31d0a6 --- /dev/null +++ b/spec/ephem/core/rotation_spec.rb @@ -0,0 +1,104 @@ +# frozen_string_literal: true + +RSpec.describe Ephem::Core::Rotation do + describe ".about_x / .about_y / .about_z" do + it "rotates coordinates about the Z axis (passive convention)" do + rotation = described_class.about_z(half_pi) + rotated = described_class.apply( + rotation, + Ephem::Core::Vector.new(1.0, 0.0, 0.0) + ) + + expect(rotated.x).to be_within(1e-12).of(0.0) + expect(rotated.y).to be_within(1e-12).of(-1.0) + expect(rotated.z).to be_within(1e-12).of(0.0) + end + + it "rotates coordinates about the X axis (passive convention)" do + rotation = described_class.about_x(half_pi) + rotated = described_class.apply( + rotation, + Ephem::Core::Vector.new(0.0, 1.0, 0.0) + ) + + expect(rotated.to_a.map { |value| value.round(12) }) + .to eq([0.0, 0.0, -1.0]) + end + + it "rotates coordinates about the Y axis (passive convention)" do + rotation = described_class.about_y(half_pi) + rotated = described_class.apply( + rotation, + Ephem::Core::Vector.new(1.0, 0.0, 0.0) + ) + + expect(rotated.to_a.map { |value| value.round(12) }) + .to eq([0.0, 0.0, 1.0]) + end + + it "returns the identity for a zero angle" do + expect(described_class.about_z(0.0)).to eq( + [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]] + ) + end + end + + describe ".multiply" do + it "composes rotations as standard matrix multiplication" do + composed = described_class.multiply( + described_class.about_z(2.5), + described_class.about_x(0.4), + described_class.about_z(1.2) + ) + + # An orthonormal rotation: composed * composed^T == identity + product = described_class.multiply(composed, composed.transpose) + identity = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]] + product.each_with_index do |row, row_index| + row.each_with_index do |value, column_index| + expect(value) + .to be_within(1e-12).of(identity[row_index][column_index]) + end + end + end + + it "is associative across the supplied factors" do + left = described_class.about_x(0.3) + middle = described_class.about_y(0.7) + right = described_class.about_z(1.1) + + all_at_once = described_class.multiply(left, middle, right) + pairwise = described_class.multiply( + described_class.multiply(left, middle), right + ) + + all_at_once.flatten.zip(pairwise.flatten).each do |combined, stepwise| + expect(combined).to be_within(1e-12).of(stepwise) + end + end + end + + describe ".apply" do + it "returns a Core::Vector" do + identity = described_class.about_z(0.0) + result = described_class.apply( + identity, + Ephem::Core::Vector.new(4.0, 5.0, 6.0) + ) + + expect(result).to be_a(Ephem::Core::Vector) + expect(result.to_a).to eq([4.0, 5.0, 6.0]) + end + + it "accepts a plain array as the vector" do + identity = described_class.about_z(0.0) + result = described_class.apply(identity, [7.0, 8.0, 9.0]) + + expect(result.to_a).to eq([7.0, 8.0, 9.0]) + end + end + + def half_pi + Math::PI / 2 + end +end diff --git a/spec/ephem/download_spec.rb b/spec/ephem/download_spec.rb index 9b2f1bf..eead576 100644 --- a/spec/ephem/download_spec.rb +++ b/spec/ephem/download_spec.rb @@ -32,6 +32,37 @@ end end + context "when downloading a NAIF binary PCK kernel" do + it "downloads and writes the PCK kernel file from NAIF" do + name = "moon_pa_de440_200625.bpc" + target_path = "tmp/kernel.bpc" + mock_content = "binary-pck-data" + file_io = StringIO.new + pathname_double = instance_double(Pathname, dirname: "tmp") + http_response = instance_double(Net::HTTPResponse) + allow(Pathname).to receive(:new) + .with(target_path) + .and_return(pathname_double) + allow(pathname_double).to receive(:open) + .with("wb") + .and_yield(file_io) + .once + naif_uri = URI.join(described_class::NAIF_PCK_BASE_URL, name) + allow(Net::HTTP).to receive(:start) + .with(naif_uri.host, naif_uri.port, use_ssl: true) + .and_yield( + double(request: nil).tap do |http| + allow(http).to receive(:request).and_yield(http_response) + end + ) + allow(http_response).to receive(:read_body).and_yield(mock_content) + + described_class.call(name: name, target: target_path) + + expect(file_io.string).to eq(mock_content) + end + end + context "when downloading an IMCCE kernel" do it "downloads, extracts, and writes the IMCCE kernel file" do name = "inpop19a.bsp" diff --git a/spec/ephem/io/daf_spec.rb b/spec/ephem/io/daf_spec.rb new file mode 100644 index 0000000..88473d8 --- /dev/null +++ b/spec/ephem/io/daf_spec.rb @@ -0,0 +1,32 @@ +# frozen_string_literal: true + +RSpec.describe Ephem::IO::DAF do + include TestSpkHelper + + describe "#file_type" do + it "identifies a DAF/SPK kernel" do + daf = described_class.new(File.open(test_spk, "rb")) + + expect(daf.file_type).to eq(:spk) + + daf.close + end + + it "identifies a binary DAF/PCK kernel" do + daf = described_class.new(File.open(moon_pa_de440_excerpt, "rb")) + + expect(daf.file_type).to eq(:pck) + + daf.close + end + + it "falls back to the integer count for a legacy NAIF/DAF SPK" do + daf = described_class.new(File.open(de405_2000_excerpt, "rb")) + + expect(daf.record_data.locator_identifier).to eq("NAIF/DAF") + expect(daf.file_type).to eq(:spk) + + daf.close + end + end +end diff --git a/spec/ephem/pck_spec.rb b/spec/ephem/pck_spec.rb new file mode 100644 index 0000000..be732a3 --- /dev/null +++ b/spec/ephem/pck_spec.rb @@ -0,0 +1,177 @@ +# frozen_string_literal: true + +RSpec.describe Ephem::PCK do + include TestSpkHelper + + describe ".open" do + it "opens a binary PCK file" do + pck = described_class.open(moon_pa_de440_excerpt) + + expect(pck).to be_a(described_class) + expect(pck.daf.file_type).to eq(:pck) + + pck.close + end + + it "rejects an SPK file" do + expect { described_class.open(test_spk) }.to raise_error( + ArgumentError, /not a binary PCK/ + ) + end + + it "closes the DAF when the kernel cannot be built" do + daf = instance_double(Ephem::IO::DAF, file_type: :pck) + allow(File).to receive(:open).and_return(instance_double(File)) + allow(Ephem::IO::DAF).to receive(:new).and_return(daf) + allow(described_class).to receive(:new) + .and_raise(Ephem::UnsupportedError) + + expect(daf).to receive(:close) + expect { described_class.open("kernel.bpc") } + .to raise_error(Ephem::UnsupportedError) + end + end + + describe "#[]" do + it "returns the orientation source for a body frame" do + pck = described_class.open(moon_pa_de440_excerpt) + + source = pck[31008] + expect(source.body).to eq(31008) + expect(source.reference_frame).to eq(1) + + pck.close + end + + it "raises KeyError for an unknown body" do + pck = described_class.open(moon_pa_de440_excerpt) + + expect { pck[42] }.to raise_error( + KeyError, /No orientation segment found for body: 42/ + ) + + pck.close + end + end + + describe "#to_s" do + it "describes the file and its segments" do + pck = described_class.open(moon_pa_de440_excerpt) + + expect(pck.to_s).to include("PCK file with 1 segments") + expect(pck.to_s).to include("orientation of frame 31008") + + pck.close + end + end + + describe "data type support" do + it "raises UnsupportedError for a non type-2 segment" do + descriptor = [0.0, 1.0, 31008, 1, 3, 1, 100] + daf = instance_double(Ephem::IO::DAF) + allow(daf).to receive(:summaries).and_return([["src", descriptor]]) + + expect { described_class.new(daf: daf) }.to raise_error( + Ephem::UnsupportedError, /Unsupported PCK data type: 3/ + ) + end + end + + describe "a body split across multiple segments" do + it "is served by an OrientationGroup" do + pck = described_class.open(moon_pa_de440_boundary_excerpt) + + group = pck[31008] + expect(group).to be_a(Ephem::Segments::OrientationGroup) + expect(group.segments.size).to eq(2) + + pck.close + end + + it "routes each time to the segment that covers it" do + pck = described_class.open(moon_pa_de440_boundary_excerpt) + group = pck[31008] + early_segment, late_segment = group.segments + early_time = early_segment.start_jd + 10 + late_time = late_segment.start_jd + 10 + + expect(group.angles_at(early_time).to_a) + .to eq(early_segment.angles_at(early_time).to_a) + expect(group.angles_at(late_time).to_a) + .to eq(late_segment.angles_at(late_time).to_a) + + pck.close + end + + it "preserves input order for an array spanning both segments" do + pck = described_class.open(moon_pa_de440_boundary_excerpt) + group = pck[31008] + early_segment, late_segment = group.segments + early_time = early_segment.start_jd + 10 + late_time = late_segment.start_jd + 10 + + results = group.angles_at([late_time, early_time]) + + expect(results[0].to_a).to eq(late_segment.angles_at(late_time).to_a) + expect(results[1].to_a).to eq(early_segment.angles_at(early_time).to_a) + + pck.close + end + + it "rejects position queries with a helpful error" do + pck = described_class.open(moon_pa_de440_boundary_excerpt) + group = pck[31008] + + expect { group.compute(2452000.0) }.to raise_error( + NotImplementedError, /angles_at or #orientation_at/ + ) + expect { group.compute_and_differentiate(2452000.0) }.to raise_error( + NotImplementedError, /orientation_at/ + ) + + pck.close + end + end + + describe "accuracy against jplephem (MOON_PA_DE440)" do + it "matches reference Euler angles and rates" do + pck = described_class.open(moon_pa_de440_excerpt) + source = pck[31008] + + reference_orientations.each do |reference| + orientation = source.orientation_at(reference[:jd]) + + orientation.to_a.zip(reference[:angles]).each do |actual, expected| + expect(actual).to be_within(1e-9).of(expected) + end + orientation.rates.zip(reference[:rates]).each do |actual, expected| + expect(actual).to be_within(1e-9).of(expected) + end + end + + pck.close + end + end + + # Ground truth from jplephem 2.24 reading moon_pa_de440_200625.bpc. + # Rates converted to ephem's per-day convention (jplephem returns per second). + def reference_orientations + [ + { + jd: 2452000.0, + angles: [-0.0645389064927909, 0.415637068726414, 2668.90512416828], + rates: [0.000371042029495874, -0.00014568869379002, 0.229623457019946] + }, + { + jd: 2456000.0, + angles: [0.0609875385648408, 0.419110500824661, 3588.67325163054], + rates: [-0.000208821705039922, 9.83855110550238e-05, 0.230149931169316] + }, + { + jd: 2460000.0, + angles: [-0.0423518935446004, 0.387741447813132, 4508.65143781339], + rates: [4.41729372107692e-05, -7.15944475461727e-05, 0.229930175592345] + } + ] + end +end diff --git a/spec/ephem/segments/orientation_segment_spec.rb b/spec/ephem/segments/orientation_segment_spec.rb new file mode 100644 index 0000000..6d27ec2 --- /dev/null +++ b/spec/ephem/segments/orientation_segment_spec.rb @@ -0,0 +1,148 @@ +# frozen_string_literal: true + +RSpec.describe Ephem::Segments::OrientationSegment do + describe "#angles_at" do + it "returns an Orientation of the three Euler angles" do + segment = create_segment_with_data + time = Ephem::Core::Constants::Time::J2000_EPOCH + + orientation = segment.angles_at(time) + + expect(orientation).to be_a(Ephem::Core::Orientation) + expect(orientation.phi).to be_within(1e-12).of(1.0) + expect(orientation.theta).to be_within(1e-12).of(2.0) + expect(orientation.psi).to be_within(1e-12).of(3.0) + end + + it "carries no rates" do + segment = create_segment_with_data + time = Ephem::Core::Constants::Time::J2000_EPOCH + + expect(segment.angles_at(time).rates?).to be(false) + end + + it "returns an Orientation per time for an array input" do + segment = create_segment_with_data + time = Ephem::Core::Constants::Time::J2000_EPOCH + + orientations = segment.angles_at([time, time + (15.0 / 86400.0)]) + + expect(orientations.length).to eq(2) + expect(orientations).to all(be_a(Ephem::Core::Orientation)) + end + + it "raises OutOfRangeError outside the coverage" do + segment = create_segment_with_data + early = Ephem::Core::Constants::Time::J2000_EPOCH - 1.0 + + expect { segment.angles_at(early) } + .to raise_error(Ephem::OutOfRangeError) + end + end + + describe "#orientation_at" do + it "returns an Orientation carrying rates" do + segment = create_segment_with_data + time = Ephem::Core::Constants::Time::J2000_EPOCH + + orientation = segment.orientation_at(time) + + expect(orientation.rates?).to be(true) + expect(orientation.rates.length).to eq(3) + end + end + + describe "#matrix_at" do + it "returns the rotation matrix from the angles" do + segment = create_segment_with_data + time = Ephem::Core::Constants::Time::J2000_EPOCH + + expect(segment.matrix_at(time)) + .to eq(segment.angles_at(time).to_matrix) + end + + it "returns one matrix per time for an array input" do + segment = create_segment_with_data + time = Ephem::Core::Constants::Time::J2000_EPOCH + + matrices = segment.matrix_at([time, time + (15.0 / 86400.0)]) + + expect(matrices.length).to eq(2) + expect(matrices).to all(be_an(Array)) + end + end + + describe "#body and #reference_frame" do + it "expose the oriented frame and the reference frame" do + segment = create_segment_with_data + + expect(segment.body).to eq(31008) + expect(segment.reference_frame).to eq(1) + end + end + + describe "#describe" do + it "describes the orientation by frame IDs" do + segment = create_segment_with_data + + expect(segment.describe).to include( + "Type 2 orientation of frame 31008 relative to frame 1" + ) + end + end + + describe "#compute" do + it "directs callers to the orientation methods" do + segment = create_segment_with_data + + expect { segment.compute(0.0) }.to raise_error( + NotImplementedError, /angles_at or #orientation_at/ + ) + expect { segment.compute_and_differentiate(0.0) }.to raise_error( + NotImplementedError, /orientation_at/ + ) + end + end + + def create_mock_daf + daf = instance_double(Ephem::IO::DAF) + + allow(daf).to receive(:read_array) do |start_word, end_word| + if end_word >= 97 + [1, 96, 11, 1] + else + [ + 0.0, # midpoint (seconds from J2000) + 43200.0, # radius (half day in seconds) + 1.0, 0.05, 0.0, # phi coefficients + 2.0, 0.1, 0.0, # theta coefficients + 3.0, 0.15, 0.0 # psi coefficients + ][(start_word - 1)..(end_word - 1)] + end + end + + allow(daf).to receive(:map_array) do |start_word, end_word| + daf.read_array(start_word, end_word) + end + + daf + end + + def create_segment_with_data + descriptor = [ + 0.0, # start_second + 86400.0, # end_second + 31008, # body (MOON_PA_DE440 frame) + 1, # reference_frame (J2000) + 2, # data_type + 1, # start_i + 100 # end_i + ] + + described_class.new( + daf: create_mock_daf, + source: "moon_pa_de440.bpc", + descriptor: descriptor + ) + end +end diff --git a/spec/ephem/segments/segment_group_spec.rb b/spec/ephem/segments/segment_group_spec.rb new file mode 100644 index 0000000..931e73a --- /dev/null +++ b/spec/ephem/segments/segment_group_spec.rb @@ -0,0 +1,99 @@ +# frozen_string_literal: true + +RSpec.describe Ephem::Segments::PositionGroup do + describe "#compute" do + it "routes a scalar time to the covering segment" do + early = build_position_segment(covers: 0.0..99.0) + late = build_position_segment(covers: 100.0..200.0) + allow(late).to receive(:compute).with(150.0, 0.0).and_return(:late) + group = described_class.new([early, late]) + + expect(group.compute(150.0)).to eq(:late) + end + + it "routes each time of an array to its segment, in order" do + early = build_position_segment(covers: 0.0..99.0) + late = build_position_segment(covers: 100.0..200.0) + allow(early).to receive(:compute).with([10.0], 0.0).and_return([:a]) + allow(late).to receive(:compute).with([150.0], 0.0).and_return([:b]) + group = described_class.new([early, late]) + + expect(group.compute([10.0, 150.0])).to eq([:a, :b]) + end + + it "raises OutOfRangeError when no segment covers the time" do + segment = build_position_segment(covers: 0.0..1.0) + group = described_class.new([segment]) + + expect { group.compute(150.0) } + .to raise_error(Ephem::OutOfRangeError) + end + end + + describe "#center and #target" do + it "delegate to the first segment" do + segment = build_position_segment(covers: 0.0..1.0) + group = described_class.new([segment]) + + expect(group.center).to eq(0) + expect(group.target).to eq(301) + end + end + + describe "#clear_data" do + it "clears data on every segment" do + first = build_position_segment(covers: 0.0..1.0) + second = build_position_segment(covers: 0.0..1.0) + group = described_class.new([first, second]) + + expect(first).to receive(:clear_data) + expect(second).to receive(:clear_data) + + group.clear_data + end + end + + def build_position_segment(covers:) + segment = instance_double( + Ephem::Segments::Segment, + center: 0, + target: 301 + ) + allow(segment).to receive(:covers?) { |time| covers.cover?(time) } + allow(segment).to receive(:clear_data) + segment + end +end + +RSpec.describe Ephem::Segments::OrientationGroup do + describe "#angles_at" do + it "routes to the covering segment" do + early = build_orientation_segment(covers: 0.0..99.0) + late = build_orientation_segment(covers: 100.0..200.0) + allow(early).to receive(:angles_at).with(10.0, 0.0).and_return(:a) + group = described_class.new([early, late]) + + expect(group.angles_at(10.0)).to eq(:a) + end + end + + describe "#body and #reference_frame" do + it "delegate to the first segment" do + segment = build_orientation_segment(covers: 0.0..1.0) + group = described_class.new([segment]) + + expect(group.body).to eq(31008) + expect(group.reference_frame).to eq(1) + end + end + + def build_orientation_segment(covers:) + segment = instance_double( + Ephem::Segments::OrientationSegment, + body: 31008, + reference_frame: 1 + ) + allow(segment).to receive(:covers?) { |time| covers.cover?(time) } + segment + end +end diff --git a/spec/ephem/segments/segment_spec.rb b/spec/ephem/segments/segment_spec.rb index 253b3cc..3161a35 100644 --- a/spec/ephem/segments/segment_spec.rb +++ b/spec/ephem/segments/segment_spec.rb @@ -78,6 +78,20 @@ expect(result.velocity).to be_a(Ephem::Core::Vector) end end + + it "pairs each time with its own position and velocity" do + first_time = Ephem::Core::Constants::Time::J2000_EPOCH + second_time = first_time + (15.0 / 86400.0) + segment = create_segment_with_data + + array_results = + segment.compute_and_differentiate([first_time, second_time]) + first_scalar = segment.compute_and_differentiate(first_time) + second_scalar = segment.compute_and_differentiate(second_time) + + expect(array_results[0].to_arrays).to eq(first_scalar.to_arrays) + expect(array_results[1].to_arrays).to eq(second_scalar.to_arrays) + end end describe "#describe" do diff --git a/spec/ephem/spk_spec.rb b/spec/ephem/spk_spec.rb index 8f83ff3..dbfd4b6 100644 --- a/spec/ephem/spk_spec.rb +++ b/spec/ephem/spk_spec.rb @@ -18,11 +18,19 @@ allow(daf_double).to receive(:summaries) .and_return([]) + allow(daf_double).to receive(:file_type) + .and_return(:spk) spk = described_class.open("spk_file.bsp") expect(spk).to be_an_instance_of(described_class) end + + it "rejects a binary PCK file" do + expect { described_class.open(moon_pa_de440_excerpt) }.to raise_error( + ArgumentError, /binary PCK file, use Ephem::PCK/ + ) + end end describe "#initialize" do diff --git a/spec/support/data/moon_pa_de440_boundary_excerpt.bpc b/spec/support/data/moon_pa_de440_boundary_excerpt.bpc new file mode 100644 index 0000000..001edc1 Binary files /dev/null and b/spec/support/data/moon_pa_de440_boundary_excerpt.bpc differ diff --git a/spec/support/data/moon_pa_de440_excerpt.bpc b/spec/support/data/moon_pa_de440_excerpt.bpc new file mode 100644 index 0000000..e0ebbbc Binary files /dev/null and b/spec/support/data/moon_pa_de440_excerpt.bpc differ diff --git a/spec/support/test_spk_helper.rb b/spec/support/test_spk_helper.rb index 7df1cc9..05ceb23 100644 --- a/spec/support/test_spk_helper.rb +++ b/spec/support/test_spk_helper.rb @@ -16,4 +16,12 @@ def de421_2000_excerpt def inpop21a_2000_excerpt File.path("#{__dir__}/data/inpop21a_2000_excerpt.bsp") end + + def moon_pa_de440_excerpt + File.path("#{__dir__}/data/moon_pa_de440_excerpt.bpc") + end + + def moon_pa_de440_boundary_excerpt + File.path("#{__dir__}/data/moon_pa_de440_boundary_excerpt.bpc") + end end