Account for light travel time

1. Iterate the delayed relative vector

When you look at an object, its light takes time to reach you. Across a room that delay is tiny; across the solar system it matters to a direction calculation. Recorded positions at a single instant therefore need another calculation before we can use them in the preparation tool’s angular fit. We begin with geometric positions from the preceding lesson, sharing a barycentric origin, EQJ axes and astronomical-unit measurements.

First find the distance between the target and Earth. Subtract Earth’s position component by component from the target’s. For a relative vector (3, 4, 0) AU, the distance is √(3² + 4² + 0²) = √25 = 5 AU. Squaring means multiplying a component by itself; adding the squares and taking the positive square root gives the vector’s length. Negative component signs do not make a distance negative.

distance = √(x² + y² + z²), in AU

Now divide distance by speed. The implementation uses c = 173.1446326846693 AU/day, so a 1 AU separation implies a delay of about 0.005775518331 day, or 499.004784 seconds. The units explain the division: AU divided by AU/day leaves days. This delay is a duration, rather than a new date.

delay in days = distance in AU / 173.1446326846693 AU/day

That first estimate uses the separation at the requested instant. But the objects move during the delay, so sampling an earlier instant can change the separation and hence the delay. The engine repeats this update, using each estimate to choose the next sampling time. This is a fixed-point iterationRepeating an estimate until it settlesClose An iteration repeats a calculation using the previous result as its next input. A fixed point is an input for which the update returns that same value.Here the delay determines the sampling time; those sampled positions determine the next delay. The routine stops when successive delays differ by less than its stated threshold, rather than requiring exact equality. Jet Propulsion Laboratory: light-time equations and reception corrections: look for a delay that scarcely changes when passed through the calculation again.

We will write the delay as τ (tau)The symbol used for the delayClose Tau is a Greek letter. In this lesson τ represents a duration in TDB days, measured backward from the requested instant. It is not a calendar date.The subscripts identify successive estimates: τ₀ is the starting estimate, τ₁ the next, and so on. . Its subscript counts the estimates: τ₁ is the first, τ₂ the second, and so on.

Follow the repeated estimate with invented straight-line motion. Write s for TDB days relative to the request, and set T(s) = (2 + 0.02s, 0, 0) AU and E(s) = (1 + 0.01s, 0, 0) AU. The coefficients 0.02 and 0.01 describe motion in AU/day. At s = 0, subtraction gives a separation of 1 AU and τ₁ = 0.005775518331 day.

For pass 2, set s = −τ₁. The target x coordinate becomes about 1.999884489633 AU, and Earth’s becomes about 0.999942244817 AU. Subtracting them gives R.x ≈ 0.999942244817 AU. Its delay is slightly smaller: τ₂ ≈ 0.005775184765 day. Repeating once more gives τ₃ ≈ 0.005775184784 day. These trajectories are teaching data, not real orbits; the next step decides when their estimate has settled enough.

Now write t for the requested TDB instant and τ for the current delay in days. Start with τ₀ = 0. On pass n, sample target T and Earth E at t − τₙ, subtract their positions to obtain Rₙ, then divide its length by c to produce τₙ₊₁. The letters T and E denote position functions of time; n counts the passes.

Rₙ = T(t − τₙ) − E(t − τₙ)

τₙ₊₁ = √(Rₙ.x² + Rₙ.y² + Rₙ.z²) / c

c = 173.1446326846693 AU/day

In apparent_lon_dec, base_offset = 0 means tdb_at(−τ) constructs that earlier instant relative to the requested one. The helper state_pos reads position only, leaving the provider’s velocity unused. Later rotations use the requested epoch’s time model, and the public runtime reads fitted angles rather than repeating this native preparation loop.

Matching origins, axes and AU units remain required before every subtraction. A provider failure propagates as an error. Earth as the target is a special case: geo_vector_eqj immediately returns (0, 0, 0), without provider calls or iteration. That zero vector has no direction to interpret.

type Position = readonly [number, number, number];
const relativePosition = (target: Position, earth: Position): Position =>
  [target[0] - earth[0], target[1] - earth[1], target[2] - earth[2]];
const delayDays = ([x, y, z]: Position): number =>
  Math.sqrt(x * x + y * y + z * z) / 173.1446326846693;

2. Stop or reject the iteration

Repeated estimation needs a stopping rule. Otherwise a preparation job could keep asking for states indefinitely, or continue with a delay outside the source margin reserved in the preceding lesson. The engine checks each new delay before deciding whether to return the relative vector from that pass.

Return to the straight-line example. Its first change is about 0.005775518331 day. The second change is about 3.33566e−7 day, and the third is about 1.92652e−11 day. The code requires a change smaller than 1e−9 day, so the second change is still too large while the third is small enough. Scientific notation 1e−9 means 0.000000001.

To measure a change regardless of its direction, subtract the previous estimate from the new one and take the absolute valueThe size of a change without its signClose Absolute value keeps a nonnegative number as it is and removes a negative number’s minus sign. Thus |−0.002| and |0.002| both equal 0.002.For successive delays, this lets the stopping rule measure how much the estimate changed regardless of whether it grew or shrank. . For example, a decrease has a negative difference, but removing the sign gives its size. The acceptance rule uses that size with a strict comparison:

accept if |τ_next − τ_previous| < 1e−9 day

We can express the threshold in more familiar time units by multiplying by the seconds in a day: 1e−9 day × 86400 s/day = 0.0000864 s, or 86.4 microseconds. This measures change between delay estimates. It is not a bound on angular error, the accuracy of the ephemeris, or agreement with observations.

reject first if τ_next > 1 day

If the change is too large, keep the new delay for the next pass. The loop allows at most ten passes, including the initial pass with zero delay. Success on pass 10 is allowed; ten unsuccessful passes return LightTimeDiverged. An unavailable state returns the provider’s error before the delay checks. The code does not separately validate finite numbers here, so these teaching traces assume finite states with the required labels.

otherwise update τ; fail after ten unsuccessful passes

Success returns the relative vector already evaluated on that pass. In our moving example, pass 3 returns the vector sampled at s = −τ₂, because τ₂ was the input to that pass. The code does not make an additional provider call at −τ₃ after accepting the new estimate. Keep that distinction when reproducing the final vector.

Two shorter examples show the boundaries. With a stationary 1 AU separation, pass 1 sets τ = 1/c and pass 2 repeats it exactly, so the change is zero and the vector returns. With an invented separation of 200 AU, the first estimate is 200/c ≈ 1.155103666 day. It fails the one-day guard on pass 1. A failed iteration provides no successful vector for the later angular fit.

type DelayDecision = 'reject' | 'accept' | 'continue';
const decideDelay = (previous: number, next: number): DelayDecision =>
  next > 1 ? 'reject' : Math.abs(next - previous) < 1e-9 ? 'accept' : 'continue';
// The caller still enforces the ten-pass budget and propagates provider errors.