Compress the angular samples

1. Choose the production series

Think of a graph with time along the horizontal axis and one measured angle along the vertical axis. Longitude and declination need separate graphs, even when they describe the same body. The builder calls each body-and-quantity combination a series. Choosing those series tells it which curves the chart runtime will need.

Sun, Moon, Mercury, Venus, Mars, Jupiter, Saturn, Uranus, Neptune, Pluto, Chiron and Ceres each receive longitude and declination. That makes 12 × 2 = 24 series. NodeOmega adds longitude only, bringing the total to 25. Derived Earth, the paired chart nodes and chart angles are assembled later; they do not add fitted series.

number of series = 12 × 2 + 1 = 25

Each graph is divided into time segments, with a polynomial describing each segment. A degree-10 polynomial needs 11 coefficients; a degree-8 polynomial needs 9. Suppose the Sun covers a hypothetical complete 64-day span in two 32-day segments. Its initial longitude and declination settings therefore require 22 + 18 = 40 coefficients. This counts only the Sun, assumes the initial fit passes, and does not estimate the full production file.

coefficients per segment = degree + 1

The initial longitude/declination settings below give segment length in days followed by degree: Sun 32/10 and 32/8; Moon and Mercury both 8/10 and 8/10; Venus and Mars 16/10 and 16/8; Jupiter and Saturn 32/10 and 32/8; Uranus, Neptune and Pluto 32/8 and 32/6; Chiron and Ceres 32/12 and 32/10. NodeOmega longitude begins at 8/10. These are repository starting choices. The retry ladder can change them if the measured differences exceed the configured limits.

The sampling functions return degrees. Body longitude follows the light-time and frame calculations already taught, while node longitude follows the same-epoch Moon state calculation. Declination remains signed. During generation, the builder samples each body’s pair together and caches results using the exact floating-point epoch bits. These kernels and mutable caches belong to preparation; the runtime needs only the finished coefficients.

2. Sample and unwrap angles

A longitude can cross zero smoothly even though its written numbers jump. If a body moves from 359° to 1°, it has advanced by 2°, not reversed by 358°. A polynomial should follow that small forward movement. Before fitting, we therefore allow longitude values to continue beyond 360° or below zero.

For raw samples 358°, 359°, 1°, 3°, the continuous sequence becomes 358°, 359°, 361°, 363°. This is called unwrappingLongitude unwrappingClose Replace jumps at zero with short signed changes between samples, retaining a continuous angle that can extend outside 0–360°. This is a representation choice for fitting, not a change of physical direction. National Institute of Standards and Technology: Digital Library of Mathematical Functions — Chebyshev recurrence. Keep the first value, then subtract each preceding raw sample from the next raw sample. Subtract 360° while that difference exceeds +180°, or add 360° while it is below −180°. Add the adjusted difference to the preceding unwrapped value. An exact +180° or −180° stays as given. Declination uses its raw signed values without this operation.

The same rule handles backward motion. Raw values 2°, 359°, 357° have short changes −3° and −2°, giving 2°, −1°, −3°. Each segment begins its own unwrap, so neighboring fitted values may differ by a whole turn while describing the same direction.

We also need sampling times that suit the polynomial shapes. Recall the normalized interval from −1 to +1. On an invented interval from day 10 to day 18, day 16 is six days into an eight-day span. Multiplication by two and subtraction of one put it at x = 2 × 6/8 − 1 = 0.5. At the start x is −1, at the midpoint it is 0, and at the end it is +1.

x = 2(t − lo)/(hi − lo) − 1

Here lo and hi are the interval endpoints, and t is the requested time. Their units cancel, so x is dimensionless. A hand example may consistently use days; production uses ephemeris seconds.

The builder takes three samples for every coefficient: degree 2 has three coefficients and therefore nine fitting times. Degree 0 has N = 3 samples, with normalized positions approximately −0.866025, 0, +0.866025. They cluster toward the ends, but none is an endpoint. To obtain the positions, imagine equally spaced angles across half a circle and take their horizontal coordinates, using cosine. The minus sign orders the times from early to late.

N = 3(d + 1); θᵢ = π(i + 0.5)/N

xᵢ = −cos(θᵢ); tᵢ = lo + (hi − lo)(xᵢ + 1)/2

Here d is the polynomial degree and N is the number of samples. The index i runs from 0 through N − 1; θᵢ is measured in radians, with π radians equal to half a turn. The half-index offset keeps θ away from both 0 and π, excluding the endpoints. These are the fitting rows. A separate grid will later check the approximation, including its endpoints. The threefold oversampling and seam rules are choices in this builder; the polynomial convention is documented by the National Institute of Standards and Technology.

const toUnit = (t: number, lo: number, hi: number): number =>
  2 * (t - lo) / (hi - lo) - 1;
const fitTimes = (lo: number, hi: number, degree: number): number[] =>
  Array.from({ length: 3 * (degree + 1) }, (_, i) => {
    const n = 3 * (degree + 1);
    const x = -Math.cos(Math.PI * (i + 0.5) / n);
    return lo + (hi - lo) * (x + 1) / 2;
  });

3. Build and solve the fit

A fit balances several supplied values with a smaller list of coefficients. Start with an invented straight-line example: at x = −1, 0, +1, the supplied angles are 2°, 3°, 6°. No single straight line passes through all three, so we choose coefficients that make the sum of squared differences as small as possible. These evenly spaced rows simplify the arithmetic; production uses the cosine grid.

For degree 1, the familiar line is p(x) = c₀ + c₁x. Each sample contributes a row [1, x]. Multiply its two entries by c₀ and c₁ and add to get the prediction. Collecting those rows into a rectangular table gives a matrix A; the supplied angles form the list y.

To combine the rows, multiply their entries in pairs and add across all three samples. The first entries give 1² + 1² + 1² = 3. Multiplying first by second entries gives −1 + 0 + 1 = 0, and squaring the second entries gives (−1)² + 0² + 1² = 2. These sums fill the table M = [[3, 0], [0, 2]].

Combine the supplied angles in the same way: 1 × 2 + 1 × 3 + 1 × 6 = 11, and (−1) × 2 + 0 × 3 + 1 × 6 = 4. This gives the right-hand list b = [11, 4]. The resulting normal equationsNormal equationsClose The system M c = b formed by adding products of sample-row entries and sample values. In matrix notation it is AᵀA c = Aᵀy. It solves linear least squares in exact arithmetic, but can amplify floating-point rounding. LAPACK Users’ Guide: linear least squares are 3c₀ = 11 and 2c₁ = 4, so the stored coefficient list is [11/3, 2].

The predictions at the three rows are 5/3°, 11/3°, 17/3°. Prediction minus supplied value gives −1/3°, +2/3°, −1/3°, whose squared sum is 2/3 square degrees. At x = 0.5 the same fitted line gives 11/3 + 2 × 0.5 = 14/3 degrees, or about 4.666667°.

For higher degrees, replace the row [1, x] with the Chebyshev shapes already introduced. Begin with T₀(x) = 1 and T₁(x) = x. Repeating Tₖ(x) = 2xTₖ₋₁(x) − Tₖ₋₂(x) gives T₂(x) = 2x² − 1 and the subsequent shapes. The subscript identifies the shape. Each coefficient still has units of degrees; the shapes themselves are dimensionless.

p(x) = sum of cₖTₖ(x), k = 0 through d

To form entry Mᵢⱼ, multiply shape i by shape j at every sample and add the products. To form bᵢ, multiply shape i by each sampled angle and add. The sums run over all rows. Solving the resulting equations gives the coefficient list c that minimizes the squared sample residuals in exact arithmetic.

Mᵢⱼ = sum of Tᵢ(x)Tⱼ(x); bᵢ = sum of Tᵢ(x)y

solve M c = b

The local solver reduces this table by subtracting multiples of one row from another. For equations 2a + b = 5 and 4a + 3b = 11, first swap rows to use 4 as the first pivot. Subtract half that row from the second: −0.5b = −0.5. Solving backward gives b = 1 and then a = 2.

This procedure is Gaussian elimination with partial pivotingElimination and partial pivotingClose Elimination subtracts row multiples to remove unknowns from lower equations. Partial pivoting first selects the largest absolute remaining entry in the current column to use as the divisor. National Institute of Standards and Technology: Digital Library of Mathematical Functions — Gaussian elimination and partial pivoting. At each column, the solver chooses the remaining row with the largest absolute entry, swaps it into the pivot position, eliminates lower entries, and then substitutes backward. The builder uses its own small solver rather than LAPACK or the earlier NumPy fitter. Production generally uses higher degrees and does not store the sample rows.

4. Assemble complete segments

Once a segment is fitted, we can store its coefficients and continue with the next segment. Before joining the results, however, we must distinguish the interval we requested from the interval the complete segments actually cover.

Take an invented span from day 10 to day 80, with 32-day segments. The 70-day span contains two complete segments: days 10–42 and 42–74. The remaining six days do not form another complete segment. The fitted end is day 74, even though the supplied high endpoint was day 80.

n = max(1, floor((hi − lo)/L + 10⁻⁶))

Here n is the segment count, L is segment length in seconds, and floor means round down. The small 10⁻⁶ allowance and minimum of one are current generator rules. Segment s starts at lo + sL and ends at lo + (s + 1)L. There is no shortened final segment: an ordinary remainder is dropped.

Within each segment, coefficients retain their shape order. Append all d + 1 coefficients from segment 0, then those from segment 1, and continue. With toy degree-1 lists [20, 2] and [24, 2], the flat stored array is [20, 2, 24, 2]. The series stores t0 = lo, L and d + 1 so the runtime can find the right segment and coefficient group.

series end = lo + nL; coefficient count = n(d + 1)

Different series, or different retry attempts, can have different actual ends. Packaging must use those ends. During assembly, the generator also measures residuals at 60 points per segment and compares adjoining endpoints; the next lesson explains those measurements. The coefficient ordering and count rules here belong to HDCHEB01 rather than the SPK format.