Fit the source position curves

1. Fit each Cartesian coordinate

The saved rows tell us where an object is at selected times. To reconstruct a position between rows, we will describe each coordinate with a short curve. Cartesian coordinates are the perpendicular x, y and z ruler readings from the preceding lessons. We fit them separately, producing one set of coefficients for each coordinate in each complete interval.

First, each supplied TDB Julian date becomes seconds from JD 2451545.0 TDB. The regeneration script divides that time span into 32-day intervals and samples every two days. A complete, aligned interval has 17 rows including both endpoints. It fits degree 12, which requires 13 coefficients per coordinate, using positions in kilometers. The downloaded velocities are not used in this fit.

Only complete intervals become records. For a 70-day span, two 32-day intervals cover 64 days, leaving six days unused. There is no shorter third record, and the kernel ends at the fitted end rather than the last input date. Generally, the script counts complete intervals by rounding down, with a small numerical allowance:

N = floor((last time − first time) / interval + 10⁻⁹)

Here floor means round down, and 10⁻⁹ is the allowance. The script rejects N < 1. Within an interval it includes samples up to 0.001 second beyond either boundary and requires at least degree + 1 samples. That count is a prerequisite for fitting, rather than a guarantee that the equations are well determined or the result accurate.

To use the same polynomial shapes in every interval, we put time on a common scale. Consider an invented interval from day 0 to day 32 after the reference. Its middle is day 16, and its half-length is 16 days. Both are 16 × 86400 = 1382400 seconds. Day 24 is 2073600 seconds: subtracting the midpoint leaves 691200 seconds, and dividing by the half-length gives 0.5.

s = (2073600 − 1382400) / 1382400 = 0.5

In general, let lo and hi be the endpoints in seconds, m their midpoint and r the half-length. The normalized time s has no unit because seconds cancel in the division. It runs from −1 at lo through 0 at m to +1 at hi.

m = (lo + hi) / 2; r = (hi − lo) / 2; s = (t − m) / r

Now we need shapes to combine into a curve. Start with three: a constant 1, a line s, and the curved expression 2s² − 1. These are the first Chebyshev polynomials, called T₀, T₁ and T₂. The subscript labels the shape; it is not multiplication. Each coefficient says how much of its shape to include.

For a small degree-2 fit, use the invented samples (s, x) = (−1, 3), (0, 3), (1, 7), where x is in kilometers. At either endpoint, 2s² − 1 is 1; at the midpoint, it is −1. Substituting the three rows into c₀ + c₁s + c₂(2s² − 1) gives:

c₀ − c₁ + c₂ = 3; c₀ − c₂ = 3; c₀ + c₁ + c₂ = 7

Subtract the first equation from the last: 2c₁ = 4, so c₁ = 2. Adding those endpoint equations gives c₀ + c₂ = 5. Together with the midpoint equation c₀ − c₂ = 3, this gives c₀ = 4 and c₂ = 1. Each coefficient is in kilometers, because the polynomial shapes have no unit.

p(s) = 4 + 2s + (2s² − 1) km

At day 24, our normalized time was 0.5. The curve gives p(0.5) = 4 + 1 + (0.5 − 1) = 4.5 km. This example fits and evaluates a degree-2 teaching curve. Production uses 13 shapes through T₁₂, with different coefficients for x, y and z.

The higher shapes can be built from ones we already have. Multiplying T₁ by 2s and subtracting T₀ gives T₂ = 2s × s − 1 = 2s² − 1. Repeating with T₂ and T₁ gives T₃ = 2s(2s² − 1) − s = 4s³ − 3s. This repeated rule is a recurrenceBuilding the next shape from earlier onesClose A recurrence is a rule that constructs the next item using preceding items. Here it builds the next Chebyshev polynomial from the two before it, starting with T₀ = 1 and T₁ = s. National Institute of Standards and Technology Digital Library of Mathematical Functions: recurrence relation and Chebyshev coefficients. Starting from T₀ = 1 and T₁ = s, use it for n = 1, 2, …, 11 to reach T₁₂.

Tₙ₊₁(s) = 2sTₙ(s) − Tₙ₋₁(s)

p(s) = c₀T₀(s) + c₁T₁(s) + … + c₁₂T₁₂(s)

Real intervals have more rows than coefficients, so a curve generally cannot hit every supplied value exactly. To see how we choose among imperfect fits, try a constant c for values 2 and 4 km. The differences, called residuals, are c − 2 and c − 4. Squaring and adding them gives (c − 2)² + (c − 4)² = 2(c − 3)² + 2. The squared part cannot be negative, so c = 3 km is best, with residuals +1 and −1 km.

This is the least-squares idea: minimize the sum of squared residuals. NumPy’s chebfit applies it to all the coefficients together, using equal weights here and a solver based on singular value decompositionThe numerical method used by the fitterClose NumPy’s documented least-squares solver uses this method to solve for the coefficients together. The lesson explains what the solver minimizes; it does not implement the decomposition. Having enough rows alone does not ensure that the coefficient equations are well determined. NumPy: chebfit least-squares fitting and coefficient order. It does not minimize the largest error, a different objective called a minimax fit. The source reports the largest absolute coordinate residual among its fitting rows; that is neither a three-dimensional vector distance nor an error bound between rows.

2. Write Type-2 SPK records

The coefficients now need a home that another reader can open. We write an SPK file, a format for storing the data needed to reconstruct celestial positions. One segment describes one target relative to one center, in one frame, over a time span. Inside it, Type 2 stores position coefficients for successive intervals; velocity can be derived by differentiating the position curve. It stores no separate velocity coefficients.

Begin with the toy degree-2 x curve from the previous step and let the y and z curves be zero. For its 32-day interval starting at ET 0, the midpoint and half-length are both 1382400 seconds. The record puts those two time values first, then all x coefficients, all y coefficients and all z coefficients:

[1382400, 1382400, 4, 2, 1, 0, 0, 0, 0, 0, 0]

Within each coordinate, the coefficients run from c₀ through c_d. The first values are named MID and RADIUS. RADIUS means the interval’s half-length in seconds, rather than an orbital distance; coefficients are in kilometers. A degree-2 record has two time values and three coefficients for each of three coordinates, making 11 numbers.

Generally, degree d needs d + 1 coefficients per coordinate. Each stored double-precision number occupies eight bytes. Degree 12 therefore needs 41 numbers, or 328 bytes, per record.

RSIZE = 2 + 3(d + 1)

After all N records, the writer appends four directory numbers: INIT, the first start in seconds; INTLEN, the interval length in seconds; RSIZE; and N. Our toy directory is [0, 2764800, 11, 1]. The record and directory together have 15 words, so their data occupies 120 bytes. In general:

W = N × RSIZE + 4; segment data bytes = 8W

Around that data sits a Double Precision Array File (DAF)The container around the coefficient dataClose DAF is the file structure that holds the segment data and the metadata needed to locate it. Its data addresses count eight-byte words starting at 1. A word address is therefore different from a byte offset. Jet Propulsion Laboratory Navigation and Ancillary Information Facility: Double Precision Array Files. This writer starts with three 1024-byte blocks: a file header, a summary block and a name block. It then writes the segment data. The header marks DAF/SPK as the file kind and LTL-IEEE as the little-endian number encoding, in which the least-significant byte comes first. Other SPK files can have more metadata blocks.

The summary records the start and fitted end in seconds and six integers: target, center, frame, type 2, first data address and last data address. Regeneration supplies Sun center 10 and frame 1, the SPICE J2000 frame identifier, for the ICRF source vectors. That frame label is metadata; writing it does not rotate the vectors. The name block holds a segment label.

Addresses count eight-byte words beginning at 1, whereas byte offsets begin at 0. Three 1024-byte blocks occupy 384 words. Data therefore begins at word 385, which has byte offset (385 − 1) × 8 = 3072. The toy’s 15 words occupy addresses 385 through 399, leaving 400 as the next free address.

byte offset = (word address − 1) × 8; last data address = 385 + W − 1

The toy file occupies 3072 + 120 = 3192 bytes. The writer does not pad its final data block to 1024 bytes. This is a serialization example, rather than a real orbit or a general SPICE compatibility certification. Validation must read the written representation back.

fitted end = INIT + N × INTLEN

3. Measure source-kernel differences

Writing the file is only the beginning of checking it. We now open that file and reconstruct positions, then compare them with Horizons positions at stated times. Comparing the written representation also exercises the storage and reading steps, instead of looking only at the coefficients before serialization.

The build requests 4000 evenly spaced comparison times from the fitted start through half a day before the fitted end. For a small example, five spots from day 0 to day 8 land at days 0, 2, 4, 6 and 8: four gaps divide the eight-day span. For n spots, the requested gap is (stop − start) / (n − 1). These times are not random, and some may coincide with fitting times. A short returned grid is rejected. Offline regeneration reads the responses from the recorded cache.

The validator opens the SPK with jplephem, selects the Sun-to-target segment and reconstructs positions at the returned dates. Let a be our reconstructed vector and b the Horizons vector. We compare them at the same time, origin, frame and units. There are two useful questions: how far apart are the positions, and how different are their directions?

Take invented vectors a = (3, 0, 0) km and b = (0, 4, 0) km. Subtracting gives (3, −4, 0) km. The right-triangle length rule gives √(9 + 16) = 5 km. The square root asks for the nonnegative number whose square is 25. With three nonzero components, we include the square of the third difference as well:

distance = |a − b| = √((aₓ − bₓ)² + (aᵧ − bᵧ)² + (a_z − b_z)²) km

The same invented vectors point along perpendicular axes, so their directions differ by a quarter turn: 90 degrees. To calculate direction differences for arbitrary vectors, we use the relationship between vector components and cosine. Imagine a point on a circle of radius 1. Its horizontal coordinate is the cosine of the angle around the circle: cos(0) = 1, cos(90°) = 0 and cos(180°) = −1.

Multiply matching vector components and add them. For our example, 3 × 0 + 0 × 4 + 0 × 0 = 0. This sum is called the dot product, written a·b. Divide it by the product of the vector lengths, 3 × 4 = 12 km², and we obtain q = 0. The kilometer-squared units cancel; q is the cosine of the direction angle.

a·b = aₓbₓ + aᵧbᵧ + a_zb_z; q = (a·b) / (|a| |b|)

For any vector, |a| means its length √(aₓ² + aᵧ² + a_z²). Once both vectors have nonzero lengths, the ratio q removes their scale and retains their direction relationship. The inverse cosineRecovering an angle from its cosineClose Inverse cosine, written arccos, reverses cosine for angles from 0 to 180 degrees. Its input must lie between −1 and 1. For example, arccos(1) = 0 and arccos(0) = π/2 radians, or 90 degrees. NumPy: inverse cosine, domain and radians, written arccos, turns that ratio back into an angle. In our example arccos(0) is a quarter turn, so it recovers the 90-degree separation.

NumPy returns that angle in radiansMeasuring an angle using a circleClose One radian is the angle subtended by a circular arc whose length equals the circle’s radius. A full turn is 2π radians or 360 degrees, with π approximately 3.14159. NumPy’s inverse cosine returns radians, so the validator must convert them to arcseconds. NumPy: inverse cosine, domain and radiansJet Propulsion Laboratory Navigation and Ancillary Information Facility: rotation mathematics and inverse-cosine sensitivity. A full turn is 2π radians, so a quarter turn is π/2 radians. The source multiplies by its rounded constant 206264.806247 arcseconds per radian, approximately (180 / π) × 3600. Our quarter turn is therefore approximately 324000 arcseconds.

angle = arccos(clamp(q, −1, 1)) × 206264.806247 arcsec

The clamp keeps q within the valid inverse-cosine range [−1, 1]. Floating-point rounding can otherwise put it just outside: for example, 1.0000000000000002 becomes 1. That safeguards the domain, but it does not remove the sensitivity of inverse cosine near nearly aligned directions. A displayed zero is not proof of exact equality.

Now compare a = (3, 0, 0) km with b = (6, 0, 0) km. Their distance is 3 km, but their dot product is 18 km² and their length product is also 18 km². Thus q = 1 and the angle is zero. Direction can agree while distance does not, which is why the report records both measures. It finds their maxima separately, with their worst epochs, and also reports means. The largest distance and angle can occur at different times.

4. Apply the source-kernel acceptance gate

The comparison produces a report; the regeneration workflow still needs a decision about whether to continue. It makes that decision using the worst sampled angular difference. There is a separate reference line in the report, so we need to distinguish what the report displays from what the workflow enforces.

Suppose the worst sampled angle is 2 arcseconds. The report compares it with 0.9375 degrees, which is 0.9375 × 3600 = 3375 arcseconds. Dividing 3375 by 2 gives a margin_factor of 1687.5. That large ratio may look reassuring, but the actual acceptance cutoff is 1 arcsecond. The result fails because 2 > 1.

report reference = 0.9375 × 3600 = 3375 arcsec

At 0.5 arcsecond, the report margin is 6750 and the gate passes. At exactly 1 arcsecond, it also passes; at 1.0001 it fails. A finite value is an ordinary number, excluding NaN, positive infinity and negative infinity. The gate requires a finite worst angle as well as an angle at or below its cutoff:

pass = isFinite(worst angle) AND worst angle ≤ 1.0 arcsec

The finite check is essential. A simple NaN > 1 comparison is false, so that comparison alone would not reject an undefined result. Valid direction angles from the preceding calculation are nonnegative; the function itself does not add a negative-input check. It also enforces no separate kilometer-distance cutoff.

The report’s ratio has its own protection against a zero reported angle. It divides by the larger of the angle and 10⁻¹² arcseconds, a tiny positive denominator. This avoids division by zero without changing the separate one-arcsecond gate.

report margin = 3375 / max(worst angle, 10⁻¹²)

In tools/regenerate.py, the order is write the SPK, validate it, then apply the gate. A failure stops regeneration before the later dataset builder and final usable provenance record. After a successful source-kernel check, the next lesson follows usable coverage and state evaluation.

const passesSourceFitGate = (worstArcsec: number): boolean =>
  Number.isFinite(worstArcsec) && worstArcsec <= 1.0;