Find the Moon orbit crossing direction

1. Form the geocentric Moon state

The goal

Input: Moon and Earth geometric states from the same provider at one requested TDB epoch, in shared-origin EQJ axes. Output: relative position in AU and relative velocity in AU/day. Both are needed to describe an instantaneous orbital plane.

Where the idea comes from

NASA’s explanation by eclipse specialist Fred Espenak connects eclipses to the Moon’s tilted orbit and the two places it crosses the ecliptic. This supplies the astronomical context for using Earth as the origin. The repository’s subtraction and provider interface have no documented individual historical inventor; they are coordinate bookkeeping, not a new eclipse model.

• NASA: Eclipses and the Moon’s orbit, by Fred Espenak

Calculate it

Geocentric means Earth-centered. Subtract Earth’s x, y and z from the corresponding Moon coordinates. Do the same for velocity: an Earth-relative rate must remove Earth’s motion as well as its position. The units of each subtraction must match.

Both state_at calls use the same requested epoch. Unlike the body-direction lesson, this path uses no light-travel-time iteration. It describes instantaneous geometry rather than the direction of received light. A provider failure is returned to the caller; no invented state replaces it.

The source then constructs AstroTime::from_tdb at that epoch. Its approximate tt drives the date-frame rotations from the previous lessons. Do not substitute a different time convention when comparing with this implementation.

r = Moon.pos(t) − Earth.pos(t) [AU]

v = Moon.vel(t) − Earth.vel(t) [AU/day]

A worked example

Synthetic shared-origin states: Moon.pos = (4, 6, 1), Earth.pos = (3, 4, 1) AU, Moon.vel = (0.1, 0.3, 0.2), Earth.vel = (0.1, 0.1, 0) AU/day. The relative state is r = (1, 2, 0) AU and v = (0, 0.2, 0.2) AU/day. These large distances simplify arithmetic; they are not a realistic lunar orbit.

type Vector = readonly [number, number, number];
const subtract = (a: Vector, b: Vector): Vector =>
  [a[0] - b[0], a[1] - b[1], a[2] - b[2]];

2. Rotate both state vectors

The goal

Input: relative r and v in EQJ. Output: both vectors expressed in true ecliptic-of-date axes (ECT). A cross product would be meaningless if its two inputs used different axes.

Where the idea comes from

NAIF’s rotation guide distinguishes rotating a vector from transforming a complete state in a moving frame. That distinction is useful here: the local helper rotates velocity components but does not add the derivative of its time-dependent rotation matrix. The reference explains the general mathematics, not the historical origin or accuracy of this implementation.

• JPL NAIF: rotations and derivatives of rotations

Calculate it

For each vector independently, apply precession From2000, then nutation From2000, then the true-obliquity rotation. Reuse the same requested epoch and tilt for position and velocity. Let Q(t) stand for that combined rotation; then the returned pair is (Q r, Q v).

For a rotated equatorial vector p and true obliquity ε in radians: x stays p.x, y becomes p.y cos ε + p.z sin ε, and z becomes −p.y sin ε + p.z cos ε. Velocity uses exactly the same arithmetic, retaining AU/day rather than AU.

The derivative of rotated position Q(t)r(t) would be Q v + Q̇ r, where Q̇ describes changing axes. eqj_to_ect_state supplies Q v only. Treat its result as velocity components expressed along date axes, not the full rate of changing date-frame coordinates. The node formula uses this implemented convention.

Our 30° example isolates the last rotation and assumes precession and nutation have already been applied. A real date uses the calculated tilt, not 30°. This is a component demonstration, not a replacement time model.

r_ECT = Q(t) r_EQJ; v_ECT = Q(t) v_EQJ

z_ECT = −y_EQD sin ε + z_EQD cos ε

d(Qr)/dt = Qv + Q̇r; this helper omits Q̇r

A worked example

Take an already rotated equatorial r = (1, 0, 0) AU and v = (0, 0, 2) AU/day, with synthetic ε = 30°. Since sin 30° = 0.5 and cos 30° = √3/2, r_ECT = (1, 0, 0) AU and v_ECT = (0, 1, √3) ≈ (0, 1, 1.732050808) AU/day. Both vectors now describe the same axes.

3. Extract the node longitude

The goal

Input: Earth-relative position and velocity in ECT. Output: ascending-node longitude wrapped toward [0, 360) degrees (with a floating-point boundary caveat below). The builder fits this derived angle separately from the source-body longitudes.

Where the idea comes from

The ascending crossing is northward through the reference plane; the opposite crossing is descending. NASA explains why these directions matter for eclipses. NAIF’s orbital-element routine also derives node information from a state, but its inertial-frame requirements and degeneracy checks differ from this repository’s helper. Neither reference establishes an individual inventor for this exact local formula; none is claimed.

• NASA: Eclipses and the Moon’s orbit, by Fred Espenak
• JPL NAIF: orbital elements from a position and velocity; degeneracy conventions

Calculate it

A cross product builds a vector perpendicular to two input vectors. For r and v, compute h.x = r.y v.z − r.z v.y, h.y = r.z v.x − r.x v.z, h.z = r.x v.y − r.y v.x. This h is specific angular momentum, with units AU²/day. Its direction describes the oriented instantaneous orbital plane.

Let k = (0, 0, 1), the positive ecliptic z direction. The intersection direction n = k × h = (−h.y, h.x, 0). For a nondegenerate state this points toward the ascending crossing: the direction of velocity there has a positive z component. Its opposite is the descending direction. This is geometry of the instantaneous state, not a prediction of a crossing time.

Use atan2(n.y, n.x) = atan2(h.x, −h.y), not atan2(h.y, h.x). atan2 keeps the quadrant; multiply its radian result by 180/π. Add or subtract whole turns to put the answer into [0, 360). The source uses Rust rem_euclid(360); the snippet uses a mathematically equivalent double remainder for ordinary finite angles. Floating-point arithmetic can round a tiny negative angle to exactly 360 in Rust rem_euclid; the JavaScript double remainder may instead yield 0. Treat those as the same direction at the seam, not bit-identical implementations.

If h.x and h.y are both zero, the projected node direction has no length and its longitude is geometrically undefined. This helper has no degeneracy rejection: it still evaluates atan2, whose signed-zero behavior can yield a numerical angle. Do not interpret that as a physical node or add a safety check to the description that the source does not have. Provider errors still propagate.

h = r × v

n = (−h.y, h.x, 0)

Ω = wrap360(atan2(h.x, −h.y) × 180/π)

A worked example

Use synthetic ECT r = (1, 0, 0) AU and v = (0, 1, 1) AU/day. Then h = (0, −1, 1) AU²/day and n = (1, 0, 0), so Ω = 0°. The object is at that crossing with positive z velocity. A second synthetic pair r = (0, −1, 0), v = (1, 0, 1) gives h = (−1, 0, 1), n = (0, −1, 0), and Ω = 270° after wrapping −90°. Neither example is an actual Moon ephemeris.

type Vector = readonly [number, number, number];
const cross = (a: Vector, b: Vector): Vector => [
  a[1] * b[2] - a[2] * b[1],
  a[2] * b[0] - a[0] * b[2],
  a[0] * b[1] - a[1] * b[0]
];
const nodeLongitude = (r: Vector, v: Vector): number => {
  const [hx, hy] = cross(r, v);
  const angle = Math.atan2(hx, -hy) * 180 / Math.PI;
  return ((angle % 360) + 360) % 360;
};