Put time and distance on a scale

1. Convert Julian dates to ephemeris seconds

A calendar organizes time into days, months and years. That makes dates convenient to read, but finding the time between two dates means accounting for different month lengths and leap years. When calculating a planet’s motion, we can simplify that bookkeeping by counting days continuously. The elapsed time between two readings then becomes a subtraction: the later day count minus the earlier one.

Astronomers call this continuous day count a Julian date (JD)Julian dateClose A Julian date is a continuous count of days with a fractional part for the time within a day. Whole-number boundaries fall at noon. For example, 2451545.25 is a quarter day after 2451545.0.The count must carry a time-scale label. A Julian date on one scale cannot be treated as a reading on another scale simply because both are day counts. Explanatory Supplement: Julian date history, section 15.1.10Jet Propulsion Laboratory: time conventions for planetary calculations. The fractional part records the time within a day, so we can describe part of a day without introducing hours and minutes into the calculation. Unlike an ordinary calendar day, each whole-number boundary falls at noon. Adding 0.25 to the count therefore moves us six hours beyond that boundary.

There is one more detail to settle before we subtract. The count tells us how a date is written, but we also need to know which time scale supplies the reading. For this calculation, the dates must use Barycentric Dynamical Time (TDB)Barycentric Dynamical TimeClose This time scale is used to describe motion in the solar system. “Barycentric” refers to the solar system’s center of mass; “dynamical” refers to calculations of motion.For this step, you need to recognize the TDB label and supply dates already expressed on that scale. Converting from everyday clock time is a separate operation. Jet Propulsion Laboratory: time conventions for planetary calculations, the scale used to model motion in the solar system. We will keep that label beside the numbers so the assumption stays visible.

Suppose our reference reading is JD 2451545.0 TDB, which is noon on January 1, 2000 on that scale. A later reading of JD 2451545.25 TDB is a quarter day beyond it. Subtraction gives 0.25 days. Each TDB day contains 86,400 seconds, so multiplying by 86,400 gives an elapsed time of 21,600 seconds, or six hours.

(2451545.25 − 2451545.0) × 86400 = 21600 seconds

The same reasoning works for dates before the reference. JD 2451544.5 TDB is half a day earlier, so subtraction gives −0.5 days and multiplication gives −43,200 seconds. The minus sign tells us which side of the reference instant we are on.

This offset is called ephemeris time (ET)Ephemeris time in this calculationClose An ephemeris is a table or model describing celestial positions over time. Here, ephemeris time means elapsed TDB seconds from noon January 1, 2000 TDB, following the convention used by the Jet Propulsion Laboratory’s planetary-calculation software.ET has also named an older astronomical time standard. In this formula it is the TDB offset, not a request to convert to that historical standard. Jet Propulsion Laboratory: time conventions for planetary calculations. We can now write the calculation for any supplied TDB Julian date. In the equation, JD_TDB is that input day count, ET is the resulting offset, and s/day means seconds per day. The source-curve fitter uses this offset to describe its time intervals.

ET = (JD_TDB − 2451545.0) × 86400 s/day

const jdToEt = (jdTdb: number): number => (jdTdb - 2451545.0) * 86400;

2. Construct the engine time model

We have turned a supplied date into seconds for the source-curve fitter. The chart calculations need another view of time, because the sky we see depends on both the motion of celestial bodies and the rotation of Earth. These processes do not give us exactly the same kind of clock.

Imagine using Earth itself as a clock, reading how far it has turned. That is the idea behind Universal Time (UT1)Universal Time measured from Earth’s rotationClose UT1 measures time through Earth’s rotation, using astronomical observations. The numeral identifies this particular form of Universal Time.Because Earth does not rotate at a perfectly constant rate, UT1 and atomic-clock time gradually diverge. UTC is kept close to UT1, but their readings are not identical. United States Naval Observatory: Universal Time. Earth’s rotation varies slightly, so this clock does not advance at a perfectly steady rate. For calculations that need a uniform time scale, astronomers also use Terrestrial Time (TT)Terrestrial TimeClose “Terrestrial” means Earth-related. This astronomical time scale is used for calculations of celestial positions as seen from Earth.In practical timekeeping, TT is International Atomic Time (TAI) plus 32.184 seconds. TAI takes its abbreviation from the French Temps atomique international. TDB and TT are distinct time scales, even though their readings are close. United States Naval Observatory: Terrestrial TimeJet Propulsion Laboratory: time conventions for planetary calculations, which follows atomic-clock time.

The difference between the two readings is called Delta T, written ΔT. The Greek letter delta means a difference here: subtract Earth-rotation time from Terrestrial Time, and express the result in seconds. This is the quantity the engine estimates when building its approximate time model.

ΔT = TT − UT1, in seconds

To follow the engine’s arithmetic, start with an instant represented by hifitime, the time library used by the runtime. Its epochEpochClose An epoch is a specified instant. It can be the reference instant from which elapsed time is counted. In hifitime, Epoch is also the name of the value that represents an instant and supports expressing it on different time scales. value can express that instant on different time scales. The engine asks it for a UTC Julian date, subtracts 2451545.0, and stores the result in a field named ut. For a supplied UTC Julian date of 2451559.0, subtraction gives ut = 14 days.

At ut = 14 days, the model estimates ΔT as 63.86 seconds; the next section explains where that number comes from. We cannot add seconds directly to a value in days. Dividing by 86,400 first gives the correction in days, which we can then add to 14.

tt = 14 + 63.86 / 86400 = 14.00073912037 days (rounded)

Written generally, the two steps are below. JD_UTC means the Julian date returned on the UTC scale, and ΔT(ut) means the model’s estimate evaluated at that day offset. Later rotation and orientation formulas consume these ut and tt fields. This example begins after the library has supplied the UTC date; the arithmetic shown here does not implement the time-scale conversion itself.

ut = JD_UTC − 2451545.0

tt = ut + ΔT(ut) / 86400

const modeledTt = (ut: number, deltaSeconds: number): number => ut + deltaSeconds / 86400;

3. Evaluate the piecewise Delta T approximation

The preceding step needed an estimate of the difference between atomic and Earth-rotation time. That difference changes over history, so one short formula would struggle to describe the whole range. The engine instead chooses among several formulas, each intended for a particular interval of years. This is called a piecewise approximation: the date determines which piece we evaluate.

Before choosing a piece, the engine turns its ut day offset into an approximate year. Return to ut = 14 days from the previous example. Subtracting 14 leaves zero, so the year is 2000. For other offsets, the engine divides the shifted day count by 365.24217 and adds it to 2000. These constants are choices in the current model, rather than a general way to convert every calendar date.

y = 2000 + (ut − 14) / 365.24217

Now suppose this calculation has given us y = 2020. That year lies in the piece covering 2005 up to, but not including, 2050. Within this piece, we measure years from 2000, so the local variable u is 2020 − 2000 = 20. The formula combines a constant, a term proportional to u, and a term proportional to u². A sum of terms like these is a polynomial.

ΔT = 62.92 + 0.32217 × 20 + 0.005589 × 20²

Work through the terms separately: 20² is 400, the linear term is 6.4434, and the squared term is 2.2356. Adding them to 62.92 gives 71.599 seconds. The coefficients are ordered [62.92, 0.32217, 0.005589]: the constant first, then the coefficient of u, then the coefficient of u².

ΔT = 62.92 + 6.4434 + 2.2356 = 71.599 seconds

We can also arrange that arithmetic from the inside out. Multiply the last coefficient by 20, add the next coefficient, multiply by 20 again, and finally add the constant. This is Horner evaluation. It gives the same answer while avoiding a separate calculation of each power, which is useful for the longer polynomials in the model.

ΔT = (0.005589 × 20 + 0.32217) × 20 + 62.92 = 71.599 seconds

The general rule for this piece is ΔT = 62.92 + 0.32217u + 0.005589u², with u = y − 2000. Each other piece has its own expression for u and its own coefficients. At y = 2000, the engine uses the earlier 1986–2005 piece. Its u is zero, leaving the constant 63.86 seconds used in the preceding section.

At a boundary, the upper limit belongs to the next piece. For y = 2005, we therefore use the 2005–2050 formula with u = 5, giving 64.670575 seconds. At exactly 2050, we move on again and must recompute u using the next piece’s rule. The complete interval table below lets you inspect those choices.

const polynomial = (u: number, coefficients: readonly number[]): number =>
  coefficients.reduceRight((acc, c) => acc * u + c, 0);
// Coefficients are ordered [a0, a1, a2, ...].

4. Convert state units

We now have time expressed in the units required by the calculations. Position and velocity need the same attention. A vector measured in kilometers and a vector measured in another length unit can describe the same location, but the numerical components will differ. We need a conversion that changes their units while preserving their signs, axes and origin.

For solar-system distances, astronomers often use the Earth–Sun distance as a convenient scale. The corresponding length unit is called an astronomical unit (AU)Astronomical unitClose An astronomical unit is a length roughly equal to the average Earth–Sun distance. Since 2012 its definition has been fixed at exactly 149,597,870,700 meters, or 149,597,870.7 kilometers.AU/day combines this length unit with a day to express velocity. The unit applies separately to each component of a position or velocity vector. International Astronomical Union: the 2012 astronomical-unit definition. Its fixed value is 149,597,870.7 kilometers. A position component of that many kilometers is therefore 1 AU; dividing a kilometer component by this value gives the component in AU.

Take an invented position vector (149597870.7, 0, 0) km. Converting each component gives (1, 0, 0) AU. A component on the other side of the origin keeps its minus sign: −299195741.4 km becomes −2 AU. These numbers illustrate the unit conversion rather than an observed orbit.

position_AU = position_km / 149597870.7

Velocity involves two units, so it needs two changes. Suppose a component is 1 kilometer per second. Dividing by 149,597,870.7 expresses the same component in AU per second. To express it per day, multiply by the 86,400 seconds in a day. Both steps are needed: the first changes the length unit, while the second changes the time unit.

1 km/s × (1 AU / 149597870.7 km) × (86400 s/day) ≈ 0.000577548327 AU/day

Notice how the units cancel in that expression: kilometers cancel with kilometers, and seconds cancel with seconds, leaving AU/day. For an invented velocity vector (0, 1, 0) km/s, the result is approximately (0, 0.000577548327, 0) AU/day. Written for any component, the conversion is:

velocity_AU/day = velocity_km/s × 86400 / 149597870.7

The engine’s KmStateProvider applies these conversions to all three components and keeps the frame label. Multiplying them by the same positive factor preserves direction. The source uses a precomputed reciprocal for the length conversion; division here makes the relationship between the units easier to follow.

There is one other kind of unit to distinguish before moving on. Small angles are often measured in arcseconds. An arcsecondArcsecondClose An arcsecond is 1/3,600 of a degree. The “arc” in its name indicates an angle, so it is distinct from a second of elapsed time. A full turn contains 360 degrees. International Bureau of Weights and Measures: units of measurement is 1/3,600 of a degree, so 7,200 arcseconds is 2 degrees. Despite the name, this measures an angle rather than elapsed time. Turning a distance into an angle requires geometry and a distance to the observer; neither of the time or length conversions above performs that calculation.

const kmToAu = (km: number): number => km / 149597870.7;
const speedToAuPerDay = (kmPerSecond: number): number => kmToAu(kmPerSecond) * 86400;