Describe the moving reference planes

1. Evaluate the implemented nutation approximation

The goal

Input: a date as TT days since J2000. Output: dpsi (Δψ, nutation in longitude) and deps (Δε, nutation in obliquity), both in arcseconds. These describe changes in reference-plane orientation, not a planet’s own longitude or declination. The next step adds deps to mean obliquity; the following lesson uses both angles in a rotation.

Where the idea comes from

Royal Museums Greenwich credits James Bradley with discovering nutation through observations between 1727 and 1748. Modern models express the wobble as sums of periodic terms. The SOFA manual attributes the abridged IAU 2000B model to McCarthy and Luzum (2003) and describes 77 terms plus fixed offsets. This repository’s function is named iau2000b but evaluates only five terms. Its name does not establish full implementation of that standard, and the standard’s published accuracy is not a claim for these five terms.

• Royal Museums Greenwich: James Bradley and nutation
• IAU SOFA manual: nut00b, obl06 and equation-of-equinoxes routines
• USNO Astronomical Almanac glossary: obliquity, nutation and equinoxes

Calculate it

Introduce trigonometry with a circle of radius 1 centered at (0, 0). Start on the positive horizontal axis and turn counterclockwise through angle a. The point reached has coordinates (cos a, sin a). Cosine is its horizontal coordinate; sine is its vertical coordinate. At 0° they are (1, 0); at 90° they are (0, 1); at 180° they are (−1, 0). Both functions repeat after a full turn, and each lies between −1 and 1. Thus A sin a + B cos a is a repeating signed correction, built from multiplication and addition.

A radian measures angle as arc length divided by circle radius. One full turn spans 2π radians = 360°; π ≈ 3.141592653589793. JavaScript Math.sin and Math.cos take radians. Convert degrees with a_rad = a_deg × π / 180. One degree has 3600 arcseconds, so a_rad = a_arcsec × π / 648000. An arcsecond is an angle, separate from a second of time.

Set t = tt / 36525, Julian centuries of TT from J2000. Form four arguments in arcseconds: l′ = 1287104.79305 + 129596581.0481t; F = 335779.526232 + 1739527262.8478t; D = 1072260.70369 + 1602961601.2090t; Ω = 450160.398036 − 6962890.5431t. Here l′ is the Sun’s mean anomaly, F is the Moon’s mean longitude minus its ascending-node longitude, D is the Moon–Sun mean elongation, and Ω is the mean longitude of the Moon’s ascending node. They are phases for periodic terms; you do not need to reconstruct orbits from them.

For each argument, take its remainder after division by 1296000 arcseconds (one turn), then multiply by ASEC2RAD = 4.848136811095359935899141e−6. Rust % and JavaScript % retain the input’s sign. A negative remainder is fine: sine and cosine give the same phase after adding a whole turn. Combine these radian arguments into a₁ = Ω, a₂ = 2(F − D + Ω), a₃ = 2(F + Ω), a₄ = 2Ω, a₅ = l′. Combined phases need no additional normalization.

For each row in the complete snippet, compute pᵢ = (Pᵢ + Ptᵢ × t) sin aᵢ + Pcᵢ cos aᵢ and eᵢ = (Eᵢ + Etᵢ × t) cos aᵢ + Esᵢ sin aᵢ. The order matters: the main longitude coefficient multiplies sine, while the main obliquity coefficient multiplies cosine. Add all five pᵢ and all five eᵢ. Coefficients are in 0.1 microarcseconds; multiplying each sum by 1e−7 converts to arcseconds. Finally add the fixed offsets: dpsi = −0.000135 + 1e−7Σpᵢ; deps = +0.000388 + 1e−7Σeᵢ. These offsets are part of this implemented approximation; do not omit them or round the per-term values before summing.

t = tt / 36525

a₁ = Ω; a₂ = 2(F − D + Ω); a₃ = 2(F + Ω); a₄ = 2Ω; a₅ = l′

pᵢ = (Pᵢ + Ptᵢt) sin aᵢ + Pcᵢ cos aᵢ

eᵢ = (Eᵢ + Etᵢt) cos aᵢ + Esᵢ sin aᵢ

Δψ = −0.000135 + 1e−7Σpᵢ arcsec; Δε = +0.000388 + 1e−7Σeᵢ arcsec

A worked example

At J2000, tt = 0 and t = 0. The radian arguments are l′ ≈ 6.240060127, F ≈ 1.627905082, D ≈ 5.198466589 and Ω ≈ 2.182439197. For row 1, sin Ω ≈ 0.8187057655 and cos Ω ≈ −0.5742132614. Its longitude contribution is (−172064161 × sin Ω + 33386 × cos Ω) × 1e−7 ≈ −14.088909133 arcsec. Rows 2–5 contribute +0.471841931, −0.221409943, −0.195030654 and −0.005182172 arcsec. The five obliquity contributions are −5.284507996, −0.535049200, +0.022753951, +0.030592192 and +0.007388527 arcsec. Sum with unrounded values, then add the offsets: dpsi ≈ −14.038824970 arcsec and deps ≈ −5.758434526 arcsec. These are outputs of the repository approximation, not measurements of the actual sky.

// Educational transcription of src/astro/frames.rs, not full IAU 2000B.
type Nutation = { dpsi: number; deps: number };
const nutationFiveTerms = (tt: number): Nutation => {
  const t = tt / 36525;
  const radians = (arcsec: number): number =>
    (arcsec % 1296000) * 4.848136811095359935899141e-6;
  const elp = radians(1287104.79305 + t * 129596581.0481);
  const f = radians(335779.526232 + t * 1739527262.8478);
  const d = radians(1072260.70369 + t * 1602961601.2090);
  const om = radians(450160.398036 - t * 6962890.5431);
  // Each row: [phase, P, Pt, Pc, E, Et, Es].
  const rows: readonly (readonly number[])[] = [
    [om, -172064161, -174666, 33386, 92052331, 9086, 15377],
    [2 * (f - d + om), -13170906, -1675, -13696, 5730336, -3015, -4587],
    [2 * (f + om), -2276413, -234, 2796, 978459, -485, 1374],
    [2 * om, 2074554, 207, -698, -897492, 470, -291],
    [elp, 1475877, -3633, 11817, 73871, -184, -1924]
  ];
  const sum = rows.reduce((total, [a, p, pt, pc, e, et, es]) => ({
    dp: total.dp + (p + pt * t) * Math.sin(a) + pc * Math.cos(a),
    de: total.de + (e + et * t) * Math.cos(a) + es * Math.sin(a)
  }), { dp: 0, de: 0 });
  return { dpsi: -0.000135 + sum.dp * 1e-7,
    deps: 0.000388 + sum.de * 1e-7 };
};

2. Compute mean and true obliquity

The goal

Input: TT centuries t plus dpsi and deps in arcseconds. Output: mobl (mean obliquity) and tobl (true obliquity) in degrees, and ee in seconds of sidereal time. Obliquity is the angle between equator and ecliptic; rotations need its radian conversion. Chart angle calculations also use true_obliquity, and sidereal_time uses ee.

Where the idea comes from

The IERS Conventions publish the IAU 2006 mean-obliquity polynomial used here. SOFA’s obl06 documentation cites Hilton and colleagues (2006) for that model. “Mean” omits the periodic nutation correction; “true” includes the implemented correction. The USNO glossary defines the equation of the equinoxes as apparent sidereal time minus mean sidereal time. This engine supplies its leading Δψ cos ε expression; it does not implement every complementary term in a full modern equation of the equinoxes.

• IERS Conventions (2010), chapter 5: mean obliquity polynomial
• IAU SOFA manual: nut00b, obl06 and equation-of-equinoxes routines
• USNO Astronomical Almanac glossary: obliquity, nutation and equinoxes

Calculate it

A polynomial is a sum of coefficients times powers: t² means t × t, and t³ means t × t × t. Compute mean obliquity in arcseconds as M = 84381.406 − 46.836769t − 0.0001831t² + 0.00200340t³ − 0.000000576t⁴ − 0.0000000434t⁵. Divide by 3600 to obtain mobl in degrees. The leading value is roughly 23.44°, not 84381 degrees. t remains TT centuries since J2000; do not use UTC years or a full Julian date.

Horner evaluation avoids constructing every power separately. Begin with the coefficient of t⁵, multiply by t, add the coefficient of t⁴, multiply by t again, and continue until adding the constant. The snippet shows the same nested expression as mean_obliq. At t = 0 all nonconstant terms vanish; at t = 1 every power is 1, making an easy arithmetic check.

True obliquity adds the signed wobble in the same unit: tobl = mobl + deps / 3600. A negative deps reduces the tilt. true_obliquity simply returns tobl from e_tilt; it introduces no additional correction. Longitude nutation dpsi is a separate angle and must not replace deps in this addition.

For ee, convert mobl from degrees to radians and evaluate dpsi × cos(mobl × π / 180) / 15. The numerator is an angle in arcseconds. A full turn is 360° = 24 hours of sidereal angle, so one hour corresponds to 15°, and one second of that time corresponds to 15 arcseconds. Dividing arcseconds by 15 therefore gives seconds of sidereal time, not hours. Divide ee by 3600 only if hours are required. In sidereal_time, 15 × ee converts it back to arcseconds before the code adds other angular terms and Earth rotation.

M = 84381.406 − 46.836769t − 0.0001831t² + 0.00200340t³ − 0.000000576t⁴ − 0.0000000434t⁵ arcsec

mobl = M / 3600 degrees

tobl = mobl + deps / 3600 degrees

ee = dpsi × cos(mobl × π / 180) / 15 seconds of sidereal time

A worked example

At J2000, M = 84381.406 arcsec and mobl = 23.439279444444°. With the previous step’s deps ≈ −5.758434526006 arcsec, add −0.001599565146° to obtain tobl ≈ 23.437679879298°. The same date’s dpsi ≈ −14.038824970274 arcsec gives ee ≈ −0.858691414656 seconds of sidereal time, or −0.000238525393 hours. One TT century later (tt = 36525), the polynomial gives M ≈ 84334.571050681 arcsec and mobl ≈ 23.426269736300°; you must recompute nutation for that date before finding its true tilt. Neither example is an independent astronomical accuracy check.

// Input t is TT centuries since J2000; corrections are arcseconds.
const meanObliquity = (t: number): number => {
  const arcsec = ((((-0.0000000434 * t - 0.000000576) * t
    + 0.00200340) * t - 0.0001831) * t - 46.836769) * t + 84381.406;
  return arcsec / 3600;
};
const earthTilt = (t: number, dpsi: number, deps: number) => {
  const mobl = meanObliquity(t);
  return { dpsi, deps, mobl, tobl: mobl + deps / 3600,
    ee: dpsi * Math.cos(mobl * 0.017453292519943296) / 15 };
};
// ee is seconds of sidereal time; mobl and tobl are degrees.