Rotate vectors into chart coordinates

1. Build the precession rotation

The goal

Input: a position vector in EQJ (equatorial J2000 axes) and the requested epoch’s AstroTime. Output: that vector in mean equatorial axes of date. Nutation will next move it to true equatorial axes of date. Rotation changes axes, not origin or AU units.

Where the idea comes from

The IAU adopted the P03 precession model in its 2006 Resolution B1. Its reference-system record names Capitaine and collaborators’ polynomial work and Hilton and collaborators’ matrix representations. That is relevant context for the polynomial-and-matrix approach below; matching this repository’s coefficients is a source check, not a claim of complete standards compliance.

• IAU Commission A3: reference-system resolutions, including P03 precession
• IAU SOFA: Earth-attitude and vector/matrix cookbooks

Calculate it

Set t = time.tt / 36525, in Julian centuries using the engine’s approximate TT-like offset from J2000. As the time lesson explained, from_tdb sets ut from UTC days and then tt = ut + modeled ΔT/86400; this is not a direct standards-grade TT conversion. A polynomial is a sum of constants times powers of t; t² means t × t. Horner form evaluates it by multiplying a running result by t and adding the next coefficient, which is how the Rust expressions are arranged. Evaluate the three polynomials in the equations below in arcseconds. ε₀ is the fixed 84381.406 arcseconds.

Convert ε₀, ψₐ, ωₐ and χₐ to radians by multiplying by ASEC2RAD = 4.848136811095359935899141e−6. Define sa = sin(ε₀), ca = cos(ε₀), sb = sin(−ψₐ), cb = cos(−ψₐ), sc = sin(−ωₐ), cc = cos(−ωₐ), sd = sin(χₐ), cd = cos(χₐ). The minus signs are part of the implementation; sin(−a) = −sin(a) while cos(−a) = cos(a).

Build all nine named coefficients using the equations below. A matrix is a rectangular table of numbers; this one has three rows and three columns. The names xx, xy and so on are source names, not instructions to use conventional row-major multiplication. From2000 stores rot = [[xx, xy, xz], [yx, yy, yz], [zx, zy, zz]]. Arrays count from zero: rot[1][0] is the first entry of the second stored row.

RotationMatrix::rotate reads down each stored column: x′ = rot[0][0]x + rot[1][0]y + rot[2][0]z; y′ = rot[0][1]x + rot[1][1]y + rot[2][1]z; z′ = rot[0][2]x + rot[1][2]y + rot[2][2]z. Thus From2000 gives x′ = xx·x + yx·y + zx·z, not xx·x + xy·y + xz·z. Prime (′) marks an output coordinate.

Into2000 instead stores [[xx, yx, zx], [xy, yy, zy], [xz, yz, zz]], exchanging rows and columns (the transpose). With the same rotate function this is the inverse rotation. Use From2000 for the forward date pipeline. In exact arithmetic rotations preserve √(x² + y² + z²); floating-point roundoff can produce tiny differences.

ψₐ = 5038.481507t − 1.0790069t² − 0.00114045t³ + 0.000132851t⁴ − 0.0000000951t⁵ (arcsec)

ωₐ = 84381.406 − 0.025754t + 0.0512623t² − 0.00772503t³ − 0.000000467t⁴ + 0.0000003337t⁵ (arcsec)

χₐ = 10.556403t − 2.3814292t² − 0.00121197t³ + 0.000170663t⁴ − 0.0000000560t⁵ (arcsec)

xx = cd·cb − sb·sd·cc; yx = cd·sb·ca + sd·cc·cb·ca − sa·sd·sc; zx = cd·sb·sa + sd·cc·cb·sa + ca·sd·sc

xy = −sd·cb − sb·cd·cc; yy = −sd·sb·ca + cd·cc·cb·ca − sa·cd·sc; zy = −sd·sb·sa + cd·cc·cb·sa + ca·cd·sc

xz = sb·sc; yz = −sc·cb·ca − sa·cc; zz = −sc·cb·sa + cc·ca

A worked example

At t = 0, ψₐ = χₐ = 0 and ωₐ = ε₀. Hence sb = sd = 0, cb = cd = 1, sc = −sa and cc = ca. Substitution gives xx = yy = zz = 1 and every off-diagonal coefficient zero, using sa² + ca² = 1. The synthetic vector (2, −1, 3) AU stays (2, −1, 3) AU. At t = 1 the polynomials give ψₐ = 5037.4014924059, ωₐ = 84381.4237831367 and χₐ = 8.173932437 arcseconds; these are polynomial checks, not measured sky positions.

type Vector = readonly [number, number, number];
type Matrix = readonly [Vector, Vector, Vector];
const rotate = (r: Matrix, [x, y, z]: Vector): Vector => [
  r[0][0] * x + r[1][0] * y + r[2][0] * z,
  r[0][1] * x + r[1][1] * y + r[2][1] * z,
  r[0][2] * x + r[1][2] * y + r[2][2] * z
];

2. Build the nutation rotation

The goal

Input: the precessed mean-equatorial vector, mean obliquity, true obliquity and nutation in longitude at the requested epoch. Output: the true-equatorial vector of date (eqd). Both angle extraction steps need this frame.

Where the idea comes from

The IAU’s 2000 Resolution B1.6 recommended a precession–nutation model and a shorter nutation version. This repository’s iau2000b routine has five periodic terms and fixed offsets, as the previous lesson showed. The name does not establish that the full recommended model is implemented. SOFA’s Earth-attitude cookbook, first released in 2007, documents the broader rotation framework.

• IAU Commission A3: reference-system resolutions, including P03 precession
• IAU SOFA: Earth-attitude and vector/matrix cookbooks

Calculate it

From e_tilt, take mobl and tobl in degrees and dpsi in arcseconds. Convert oblm = mobl × DEG2RAD, oblt = tobl × DEG2RAD, and ψ = dpsi × ASEC2RAD. Set cobm = cos(oblm), sobm = sin(oblm), cobt = cos(oblt), sobt = sin(oblt), cpsi = cos(ψ), spsi = sin(ψ). Each trigonometric input must be radians.

Evaluate the nine coefficients below. Forward nutation stores [[xx, xy, xz], [yx, yy, yz], [zx, zy, zz]] and uses the same column-reading rotate rule. Its enum label is From2000, but here the input has already been precessed: this step goes from mean equatorial of date to true equatorial of date, not directly from EQJ.

Into2000 stores [[xx, yx, zx], [xy, yy, zy], [xz, yz, zz]] and is the inverse nutation rotation. To reverse the entire pipeline, undo nutation first and precession second. Matrix order matters: apply precession to the vector, then nutation to that result. Nutation alone does not supply the long-term precession change.

xx = cpsi; yx = −spsi·cobm; zx = −spsi·sobm

xy = spsi·cobt; yy = cpsi·cobm·cobt + sobm·sobt; zy = cpsi·sobm·cobt − cobm·sobt

xz = spsi·sobt; yz = cpsi·cobm·sobt − sobm·cobt; zz = cpsi·sobm·sobt + cobm·cobt

eqd = nutation_From2000(precession_From2000(eqj))

A worked example

For easy arithmetic, invent mobl = tobl = 0° and dpsi = 108000 arcseconds = 30°. These are synthetic angles, not an Earth-tilt prediction. The stored matrix is [[cos30°, sin30°, 0], [−sin30°, cos30°, 0], [0, 0, 1]]. Applied to (2, 1, 0) AU, it gives x′ = 2×0.866025404 − 0.5 = 1.232050808 and y′ = 2×0.5 + 0.866025404 = 1.866025404 AU. The output length is still √5 AU. With dpsi = 0 and mobl = tobl, the matrix is identity.

3. Extract ecliptic longitude

The goal

Input: the original EQJ vector and requested time. Output from ecliptic_of_date: longitude and latitude in degrees and distance in AU. The builder keeps longitude for fitting; latitude and distance are intermediate outputs.

Where the idea comes from

SOFA publishes tools for both Cartesian vectors and spherical angles, and its miscellaneous cookbook covers ecliptic coordinates. This provides a contemporary reference for representing one direction in different coordinate systems. The repository’s exact sign conventions and guards come from frames.rs; no historical discoverer is assigned to these local implementation choices.

• IAU SOFA: Earth-attitude and vector/matrix cookbooks

Calculate it

First calculate eqd by precession and nutation. Let ε = true obliquity in radians. Tilt the equatorial axes into ecliptic axes with ex = eqd.x, ey = eqd.y·cosε + eqd.z·sinε, ez = −eqd.y·sinε + eqd.z·cosε. This rotation uses true, not mean, obliquity.

Compute xyproj = hypot(ex, ey) = √(ex² + ey²), the length of the projection onto the ecliptic plane. atan2(y, x) returns the signed angle of the pair (x, y), choosing the correct quadrant; its argument order is y first. Unlike atan(y/x), it can distinguish opposite directions and handle x = 0. JavaScript Math.atan2 returns radians in [−π, π]. Multiply by RAD2DEG = 57.295779513082321 to obtain degrees.

Start longitude at 0. If xyproj > 0, calculate longitude = RAD2DEG × atan2(ey, ex) and add 360° if it is negative. A vector straight above or below the ecliptic has no unique longitude; the source deliberately leaves it at 0. This guard is a convention, not a measured direction.

Calculate latitude = RAD2DEG × atan2(ez, xyproj), and distance = √(ex² + ey² + ez²). Latitude is the angle above or below the ecliptic plane; it lies between −90° and +90° for finite nonzero vectors. Distance retains the input AU units. A zero vector reaches longitude 0 and distance 0; the atan2 latitude is numerically evaluated, although a zero vector has no physical direction. The functions add no separate finite-value validation.

(ex, ey, ez) = (eqd.x, eqd.y·cosε + eqd.z·sinε, −eqd.y·sinε + eqd.z·cosε)

longitude = atan2(ey, ex) × RAD2DEG, wrapped if negative (when xyproj > 0)

latitude = atan2(ez, xyproj) × RAD2DEG; distance = √(ex² + ey² + ez²)

A worked example

Use a synthetic already-rotated eqd = (1, −√3, 1) AU and an invented ε = 0°, so the ecliptic vector is unchanged. xyproj = √(1 + 3) = 2 AU. atan2(−√3, 1) = −60°, so longitude is 300°. Latitude is atan2(1, 2) ≈ 26.565051177° and distance is √5 ≈ 2.236067977 AU. A different synthetic ecliptic vector (0, 0, 2) has xyproj = 0, longitude 0 by guard, latitude +90°, and distance 2 AU.

type EclipticVector = readonly [number, number, number];
const eclipticAngles = ([x, y, z]: EclipticVector): readonly [number, number, number] => {
  const xy = Math.hypot(x, y);
  const raw = xy > 0 ? Math.atan2(y, x) * 57.295779513082321 : 0;
  const lon = raw < 0 ? raw + 360 : raw;
  return [lon, Math.atan2(z, xy) * 57.295779513082321, Math.sqrt(x * x + y * y + z * z)];
};

4. Extract equatorial declination

The goal

Input: the original EQJ vector and requested time. Output: declination in degrees in the true equatorial frame of date. The builder stores this alongside ecliptic longitude; these two angles belong to different planes.

Where the idea comes from

SOFA’s vector/matrix cookbook covers conversions between Cartesian components and spherical angles. This is the mathematical setting for recovering an angle from a component ratio. The particular eqd_declination guard belongs to this implementation, and the source supplies no separate historical attribution for it.

• IAU SOFA: Earth-attitude and vector/matrix cookbooks

Calculate it

Repeat precession then nutation to obtain eqd. Do not apply the ecliptic tilt here: declination is measured from the true equator, not from the ecliptic. Compute r = √(eqd.x² + eqd.y² + eqd.z²), then the dimensionless ratio eqd.z / r.

asin means inverse sine: it finds the angle whose sine equals the supplied ratio. Its principal result runs from −π/2 to +π/2 radians (−90° to +90°). For example, asin(0.5) = π/6 radians = 30°. Use radians internally and multiply by RAD2DEG for the output.

If r > 0, return RAD2DEG × asin(eqd.z / r). Otherwise return 0. Thus a zero vector receives a numeric zero even though it has no meaningful declination. The source does not clamp the ratio into [−1, 1] or add a separate nonfinite-input check. The examples assume finite vectors with a valid ratio; do not silently introduce extra guards into an explanation of the source.

r = √(eqd.x² + eqd.y² + eqd.z²)

declination = RAD2DEG × asin(eqd.z / r) if r > 0; otherwise 0

A worked example

For the synthetic true-equatorial vector eqd = (√3, 0, 1) AU, r = √(3 + 1) = 2 AU and z/r = 0.5, so declination is +30°. Multiplying every component by 4 gives r = 8 AU and z = 4 AU, retaining the same declination. A synthetic vector (0, 0, −2) gives −90°; (0, 0, 0) returns 0 by the guard.

type EquatorialVector = readonly [number, number, number];
const declination = ([x, y, z]: EquatorialVector): number => {
  const r = Math.sqrt(x * x + y * y + z * z);
  return r > 0 ? 57.295779513082321 * Math.asin(z / r) : 0;
};

5. Produce angular samples

The goal

Input: a provider, target body, requested TDB epoch and a function constructing offset epochs. Output: (normalized ecliptic longitude, true-equatorial declination), or the propagated provider/light-time error. These samples feed native dataset fitting; the public runtime evaluates the resulting coefficients.

Where the idea comes from

JPL NAIF distinguishes reception light-time correction from stellar aberration and other effects. Its usual reception equation keeps the observer at the requested time. That distinction helps describe what this engine actually computes. The local apparent_lon_dec name does not establish a complete apparent-place model or compliance with NAIF’s correction recipes.

• JPL NAIF: reception light time and aberration corrections

Calculate it

Build time = AstroTime::from_tdb(requested_tdb). Call geo_vector_eqj with base_offset = 0. As in the light-time lesson, each pass subtracts target and Earth positions sampled at the same delayed epoch t − τ. Both states are in shared barycentric EQJ axes before subtraction. Return the vector from the successful pass; propagate provider errors or LightTimeDiverged.

Use that returned vector with the original requested time for both ecliptic_of_date and eqd_declination. The positions are delayed, but the date-frame precession, nutation and true obliquity are evaluated at the request. Do not replace the requested rotation epoch with the delayed sampling epoch.

ecliptic_of_date returns (lon, lat, dist). Discard lat and dist here. Normalize lon with the source’s Euclidean remainder by 360°: the intended range is [0°, 360°), with floating-point roundoff at the seam: a tiny negative remainder plus 360 can round to exactly 360. Treat that as the same direction as 0°. For example −10° becomes 350° and 370° becomes 10°. The ecliptic helper already wraps its usual atan2 result, so this final normalization also defines the wrapper’s boundary behavior.

Calculate declination separately from the same EQJ vector and requested time, repeating the equatorial rotations. Return (normalized lon, declination). This routine adds no stellar aberration, gravitational deflection, topocentric observer displacement or atmospheric refraction. Do not add those corrections to a source-faithful teaching trace.

Earth is special: geo_vector_eqj returns the zero vector immediately. The angle helpers then supply numeric zero values for the ordinary all-positive-zero components. This represents the source’s fallback, not an observable Earth direction from Earth. A failure in light time supplies no angular sample to fit.

time = AstroTime::from_tdb(requested_tdb)

geo = geo_vector_eqj(provider, body, base_offset = 0, tdb_at)

(lon, unused_lat, unused_dist) = ecliptic_of_date(geo, time)

sample = (normalize360(lon), eqd_declination(geo, time))

A worked example

For a deliberately simplified trace, suppose a successful delay pass returns geo = (1, −1, 0) AU and suppose the requested date rotations are identity with invented true obliquity 0°. These rotations are synthetic teaching assumptions, not the actual Earth model at any claimed date. Longitude is −45° before wrapping, hence 315°; declination is 0°. The sample is (315°, 0°). If the delay was 0.006 day, the positions were sampled at t − 0.006 day, while the assumed frame rotations still belong to t.

const normalize360 = (degrees: number): number => {
  const remainder = degrees % 360;
  return remainder < 0 ? remainder + 360 : remainder;
};
// Finite-input teaching equivalent of the Rust Euclidean wrapping step.