The moon's position¶

Its right ascension and declination.

Requires

the time, as a Julian day (UT).

Provides

RightAscension and Declination (rad), the ecliptic longitude and latitude (rad).

Law

the mean longitude plus the six largest periodic terms in longitude and four in latitude; to a few tenths of a degree.

Source

the Astronomical Almanac's low-precision formulae (section D).

Proof

reference, invariant.

No Trait form: not written yet: its longitude, latitude and phase are lusine.sky.moon_longitude, moon_latitude and moon_phase; the obliquity and the turn from the ecliptic to the equator are not written, and a Julian day does not fit binary32 (sidereal's note); until then a reference.

sky.moon(jd)

The kernel follows the Almanac's low-precision series. Every number it runs on is a row's: the Moon's elements lusine.bodies.moon and the Earth's obliquity lusine.bodies.earth, bound as a game binds them: the kernel declares no constant.

In [1]:

import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
    root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import dmath as D, library
from lusine.kernels.sky._astro import julian_day, to_deg
from lusine.kernels.sky.moon import moon
from lusine.kernels.sky.sun import sun

import os, sys root = os.getcwd() while not os.path.isfile(os.path.join(root, 'docsite.json')): root = os.path.dirname(root) sys.path.insert(0, os.path.join(root, 'tools')) import numpy as np import matplotlib.pyplot as plt from lusine import dmath as D, library from lusine.kernels.sky._astro import julian_day, to_deg from lusine.kernels.sky.moon import moon from lusine.kernels.sky.sun import sun

The reference¶

The kernel's reference, in numpy: the body under the declaration above, and what every section below runs. It stays as the kernel's test beside its Trait form.

In [2]:

library.reference(moon)

library.reference(moon)

Out[2]:

the reference, moon.py

def moon(jd):
    """THE MOON'S POSITION: its right ascension and declination.

    Requires: the time, as a Julian day (UT).
    Provides: `RightAscension` and `Declination` (rad), the ecliptic longitude
    and latitude (rad).
    Law: the mean longitude plus the six largest periodic terms in longitude and
    four in latitude; to a few tenths of a degree.
    Source: the Astronomical Almanac's low-precision formulae (section D).
    Proof: reference, invariant."""
    m = _MOON
    t = (np.asarray(jd, dtype=np.float64) - m["epoch"]) / 36525.0

    def terms(kind):
        # the body's periodic terms in one coordinate, lusine.body_terms' rows (library_sources.bodies.TERMS)
        return sum(amplitude * _s(phase + rate * t) for coordinate, amplitude, phase, rate in bodies.TERMS
                   if coordinate == kind)
    lam = _wrap360(m["longitude_epoch"] + m["longitude_rate"] * t + terms("longitude"))
    beta = terms("latitude")
    e = _obliquity(np.asarray(jd, dtype=np.float64) - _EARTH["epoch"], _EARTH)
    ra = D.datan2(_s(lam) * _c(e) - (_s(beta) / _c(beta)) * _s(e), _c(lam))
    dec = D.dasin(_s(beta) * _c(e) + _c(beta) * _s(e) * _s(lam))
    return {"ra": ra, "dec": dec, "longitude": deg(lam), "latitude": deg(beta)}

Reference: the eclipse of 8 April 2024¶

At greatest eclipse (18:17 UT) the moon stood in front of the sun: their positions agree to within the discs' size.

In [3]:

jd = julian_day(2024, 4, 8, 18 + 17 / 60)
m, s = moon(jd), sun(jd)
sep = float(to_deg(np.arccos(np.clip(np.sin(m['dec']) * np.sin(s['dec']) + np.cos(m['dec']) * np.cos(s['dec']) * np.cos(m['ra'] - s['ra']), -1, 1))))
print(f'moon to sun: {sep:.3f}°')
assert sep < 0.5

jd = julian_day(2024, 4, 8, 18 + 17 / 60) m, s = moon(jd), sun(jd) sep = float(to_deg(np.arccos(np.clip(np.sin(m['dec']) * np.sin(s['dec']) + np.cos(m['dec']) * np.cos(s['dec']) * np.cos(m['ra'] - s['ra']), -1, 1)))) print(f'moon to sun: {sep:.3f}°') assert sep < 0.5

moon to sun: 0.425°

Invariant: within the moon's band, round the sky in a month¶

Its ecliptic latitude stays within the orbit's 5.1° inclination, and its longitude advances about 13.2° a day.

In [4]:

days = np.arange(0, 365, 0.5)
m = moon(julian_day(2024, 1, 1) + days)
assert np.degrees(np.max(np.abs(m['latitude']))) < 5.35
step = np.diff(np.unwrap(m['longitude'])) / 0.5
rate = float(np.degrees(np.mean(step)))
print(f'{rate:.2f}° a day')
assert abs(rate - 13.18) < 0.1
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(days, np.degrees(m['latitude']))
ax.set_xlabel('day of 2024')
ax.set_ylabel('ecliptic latitude (°)')
plt.show()

days = np.arange(0, 365, 0.5) m = moon(julian_day(2024, 1, 1) + days) assert np.degrees(np.max(np.abs(m['latitude']))) < 5.35 step = np.diff(np.unwrap(m['longitude'])) / 0.5 rate = float(np.degrees(np.mean(step))) print(f'{rate:.2f}° a day') assert abs(rate - 13.18) < 0.1 fig, ax = plt.subplots(figsize=(7, 3)) ax.plot(days, np.degrees(m['latitude'])) ax.set_xlabel('day of 2024') ax.set_ylabel('ecliptic latitude (°)') plt.show()

13.16° a day