The Moon's longitude¶
The Moon's geocentric ecliptic longitude, its mean longitude (lusine.sky.mean_longitude over its row of lusine.bodies) and the sum of its periodic terms in longitude (its rows of lusine.body_terms whose coordinate is the longitude, each one lusine.sky.periodic_term), lambda = L + sum A sin(phase + rate t), taken within a turn.
The time is given since the Moon's epoch (its row's, J2000.0), in seconds, since a binary32 Julian day cannot hold the hours. Good to a few tenths of a degree. Source: The Astronomical Almanac, section D, low-precision formulae for the Moon; Meeus, Astronomical Algorithms, ch. 47. Checked against a reference in numpy, which has its own tests.
- Source
- The Astronomical Almanac, section D, low-precision formulae for the Moon; Meeus, Astronomical Algorithms, ch. 47.
- Proof
- reference, invariant.
moon_longitude.trait
namespace lusine.sky;
float<rad> moon_longitude(float<s> since_epoch)
{
let moon = lusine.bodies["moon"];
let terms = sum(lusine.sky.periodic_term(since_epoch, lusine.body_terms.amplitude, lusine.body_terms.phase,
lusine.body_terms.rate))
where lusine.body_terms.body == moon && lusine.body_terms.coordinate == lusine.coordinates["longitude"];
let l = lusine.sky.mean_longitude(since_epoch, moon.longitude_epoch, moon.longitude_rate) + terms;
return l - 360.0 deg * floor(l / 360.0 deg);
}
The Moon's geocentric ecliptic longitude, its mean longitude (lusine.sky.mean_longitude over its row of lusine.bodies) and its six largest periodic terms, the Almanac's low-precision series. One kernel sums the terms: a fold over the Moon's rows of lusine.body_terms whose coordinate is the longitude (a row of lusine.coordinates), applying lusine.sky.periodic_term to each. Every number is a row's; the kernel declares none but the turn.
In [1]:
Copied!
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 trait_run
from lusine import rows, binary32
F32 = np.float32
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 trait_run
from lusine import rows, binary32
F32 = np.float32
The kernel¶
The law as the library states it, and as this page runs it.
In [2]:
Copied!
print(trait_run.source('sky/moon_longitude.trait'))
print(trait_run.source('sky/moon_longitude.trait'))
// lusine.sky.moon_longitude: the Moon's geocentric ecliptic longitude, its mean longitude (lusine.sky.mean_longitude
// over its row of lusine.bodies) and the sum of its periodic terms in longitude (its rows of lusine.body_terms whose
// coordinate is the longitude, each one lusine.sky.periodic_term), lambda = L + sum A sin(phase + rate t), taken within
// a turn. The time is given since the Moon's epoch (its row's, J2000.0), in seconds, since a binary32 Julian day
// cannot hold the hours. Good to a few tenths of a degree.
// Source: The Astronomical Almanac, section D, low-precision formulae for the Moon; Meeus, Astronomical Algorithms,
// ch. 47.
// Checked against a reference in numpy, which has its own tests.
namespace lusine.sky;
float<rad> moon_longitude(float<s> since_epoch)
{
let moon = lusine.bodies["moon"];
let terms = sum(lusine.sky.periodic_term(since_epoch, lusine.body_terms.amplitude, lusine.body_terms.phase,
lusine.body_terms.rate))
where lusine.body_terms.body == moon && lusine.body_terms.coordinate == lusine.coordinates["longitude"];
let l = lusine.sky.mean_longitude(since_epoch, moon.longitude_epoch, moon.longitude_rate) + terms;
return l - 360.0 deg * floor(l / 360.0 deg);
}
The program¶
A lattice of items, each the host's time since the Moon's epoch (its row of lusine.bodies, J2000.0) in seconds, and main() binding both kernels over them. A binary32 Julian day cannot hold the hours (two days near 2,451,545 are a quarter of a day apart), so the time is given since the epoch. moon(jds) runs it once for a set of dates. The helpers below model the kernels in binary32 for the bit-for-bit section: term is lusine.sky.periodic_term rounded after each operation, its sin the library's (lusine.binary32's sin, the language's sequence), and tree Trait's fixed summing order.
In [3]:
Copied!
MOON = """namespace sky;
lusine.version version = { { number = 1 } };
struct item { host float<s> since_epoch; float<rad> moon_longitude; float<rad> moon_latitude; };
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main()
{
items.moon_longitude = lusine.sky.moon_longitude(items.since_epoch);
items.moon_latitude = lusine.sky.moon_latitude(items.since_epoch);
}
"""
print(MOON)
J2000 = 2451545.0 # the Moon's epoch, its row of lusine.bodies, as a Julian day
def moon(jds):
"""(longitude, latitude), rad, the kernels' at each Julian day, and the run; the time is given as binary32 seconds since the epoch."""
since = F32((np.atleast_1d(np.asarray(jds, dtype=float)) - J2000) * 86400.0)
run = trait_run.Program(['sky/moon_longitude.trait', 'sky/moon_latitude.trait'], MOON,
{'sky.items': [len(since)]})
run.store('sky.items', 'since_epoch', since)
lon, lat = run.steps('sky.main', 1, [('sky.items', 'moon_longitude'), ('sky.items', 'moon_latitude')])
return lon[0], lat[0], since, run
def tree(values):
"""A binary32 sum in Trait's fixed tree: the first ceil(n / 2) members, then the rest, each summed so."""
if len(values) == 1:
return values[0]
mid = (len(values) + 1) // 2
return F32(tree(values[:mid]) + tree(values[mid:]))
def term(r, since):
"""lusine.sky.periodic_term of a row, A sin(phase + rate t), in binary32 with the library's sin."""
return F32(F32(r['amplitude']) * binary32.evaluate('sin', F32(F32(r['phase']) + F32(F32(r['rate']) * since))))
def terms(run, coordinate):
"""The Moon's rows of lusine.body_terms in one coordinate, as the program holds them, in slot order."""
return [r for r in rows.rows(run.unit, 'lusine.body_terms') if r['key'].startswith(f'moon, {coordinate}, ')]
def ulp(x):
return np.spacing(np.abs(np.asarray(x, dtype=F32))).astype(float)
def published(jds):
"""astronomy-engine's geocentric ecliptic longitude and latitude of the Moon, deg."""
import astronomy
e = [astronomy.EclipticGeoMoon(astronomy.Time(float(jd) - J2000)) for jd in jds]
return np.array([x.lon for x in e]), np.array([x.lat for x in e])
CENTURY = J2000 + np.linspace(-36525.0, 36525.0, 4001) # 1900 to 2100, every 18 days
DATES = np.arange(2447892.5, 2469807.5, 3.7) # 1990 to 2050, every 3.7 days
ECLIPSE = 2460409.2618 # 8 April 2024, 18:17 UT, greatest eclipse
MOON = """namespace sky;
lusine.version version = { { number = 1 } };
struct item { host float<s> since_epoch; float<rad> moon_longitude; float<rad> moon_latitude; };
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main()
{
items.moon_longitude = lusine.sky.moon_longitude(items.since_epoch);
items.moon_latitude = lusine.sky.moon_latitude(items.since_epoch);
}
"""
print(MOON)
J2000 = 2451545.0 # the Moon's epoch, its row of lusine.bodies, as a Julian day
def moon(jds):
"""(longitude, latitude), rad, the kernels' at each Julian day, and the run; the time is given as binary32 seconds since the epoch."""
since = F32((np.atleast_1d(np.asarray(jds, dtype=float)) - J2000) * 86400.0)
run = trait_run.Program(['sky/moon_longitude.trait', 'sky/moon_latitude.trait'], MOON,
{'sky.items': [len(since)]})
run.store('sky.items', 'since_epoch', since)
lon, lat = run.steps('sky.main', 1, [('sky.items', 'moon_longitude'), ('sky.items', 'moon_latitude')])
return lon[0], lat[0], since, run
def tree(values):
"""A binary32 sum in Trait's fixed tree: the first ceil(n / 2) members, then the rest, each summed so."""
if len(values) == 1:
return values[0]
mid = (len(values) + 1) // 2
return F32(tree(values[:mid]) + tree(values[mid:]))
def term(r, since):
"""lusine.sky.periodic_term of a row, A sin(phase + rate t), in binary32 with the library's sin."""
return F32(F32(r['amplitude']) * binary32.evaluate('sin', F32(F32(r['phase']) + F32(F32(r['rate']) * since))))
def terms(run, coordinate):
"""The Moon's rows of lusine.body_terms in one coordinate, as the program holds them, in slot order."""
return [r for r in rows.rows(run.unit, 'lusine.body_terms') if r['key'].startswith(f'moon, {coordinate}, ')]
def ulp(x):
return np.spacing(np.abs(np.asarray(x, dtype=F32))).astype(float)
def published(jds):
"""astronomy-engine's geocentric ecliptic longitude and latitude of the Moon, deg."""
import astronomy
e = [astronomy.EclipticGeoMoon(astronomy.Time(float(jd) - J2000)) for jd in jds]
return np.array([x.lon for x in e]), np.array([x.lat for x in e])
CENTURY = J2000 + np.linspace(-36525.0, 36525.0, 4001) # 1900 to 2100, every 18 days
DATES = np.arange(2447892.5, 2469807.5, 3.7) # 1990 to 2050, every 3.7 days
ECLIPSE = 2460409.2618 # 8 April 2024, 18:17 UT, greatest eclipse
namespace sky;
lusine.version version = { { number = 1 } };
struct item { host float<s> since_epoch; float<rad> moon_longitude; float<rad> moon_latitude; };
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main()
{
items.moon_longitude = lusine.sky.moon_longitude(items.since_epoch);
items.moon_latitude = lusine.sky.moon_latitude(items.since_epoch);
}
Reference: astronomy-engine, 1990 to 2050¶
Against astronomy-engine's Moon (its full lunar theory, pinned), every 3.7 days for sixty years, the kernel is within 0.4°: the series' few tenths of a degree, the six terms it keeps. At the eclipse of 8 April 2024 it stood in front of the Sun: astronomy-engine's longitude there within 0.4°.
In [4]:
Copied!
lon, _, _, _ = moon(DATES)
want, _ = published(DATES)
miss = (np.degrees(lon) - want + 180.0) % 360.0 - 180.0
print(f'worst {np.max(np.abs(miss)):.3f} degrees over {len(DATES)} dates')
assert np.max(np.abs(miss)) < 0.4
at, _, _, _ = moon([ECLIPSE])
sun, _ = published([ECLIPSE])
print(f'8 April 2024, 18:17 UT: {np.degrees(at[0]):.3f} degrees; astronomy-engine {sun[0]:.3f}')
assert abs((np.degrees(at[0]) - sun[0] + 180.0) % 360.0 - 180.0) < 0.4
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(DATES[:300] - DATES[0], miss[:300])
ax.set_xlabel('days from 1 January 1990')
ax.set_ylabel('kernel - astronomy-engine (°)')
plt.show()
lon, _, _, _ = moon(DATES)
want, _ = published(DATES)
miss = (np.degrees(lon) - want + 180.0) % 360.0 - 180.0
print(f'worst {np.max(np.abs(miss)):.3f} degrees over {len(DATES)} dates')
assert np.max(np.abs(miss)) < 0.4
at, _, _, _ = moon([ECLIPSE])
sun, _ = published([ECLIPSE])
print(f'8 April 2024, 18:17 UT: {np.degrees(at[0]):.3f} degrees; astronomy-engine {sun[0]:.3f}')
assert abs((np.degrees(at[0]) - sun[0] + 180.0) % 360.0 - 180.0) < 0.4
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(DATES[:300] - DATES[0], miss[:300])
ax.set_xlabel('days from 1 January 1990')
ax.set_ylabel('kernel - astronomy-engine (°)')
plt.show()
worst 0.361 degrees over 5923 dates 8 April 2024, 18:17 UT: 19.541 degrees; astronomy-engine 19.360
Reference: bit for bit, and the float64 series within binary32's rounding¶
From 1900 to 2100 the kernel's longitude is a binary32 model's bit for bit: the mean longitude L0 + n t taken within a turn, the six terms summed in Trait's tree, the sum taken within a turn, each operation rounded once, from the rows as the program holds them. The numpy reference, the same series in float64 at the same binary32 time, agrees within 4 ulp of the mean longitude's argument |L0| + |n t|: n t rounds once and n was rounded once (an ulp), L0 and the wrap's turn add half an ulp each, and the terms, each under 0.11 rad, add less. A century from the epoch n t is 8,400 rad, whose ulp is a thousandth of a radian: that is 0.06°, the price of binary32 time a century out.
In [5]:
Copied!
sys.path.insert(0, os.path.join(root, 'kernels'))
from lusine.kernels.sky.moon_longitude import moon_longitude
lon, _, since, run = moon(CENTURY)
body = rows.row(run.unit, 'lusine.bodies', 'moon')
TURN = F32(2.0 * np.pi)
def within_turn(x):
return F32(x - F32(TURN * F32(np.floor(F32(x / TURN)))))
mean = within_turn(F32(F32(body['longitude_epoch']) + F32(F32(body['longitude_rate']) * since)))
model = within_turn(F32(mean + tree([term(r, since) for r in terms(run, 'longitude')])))
assert np.array_equal(F32(lon), model)
ref = moon_longitude(since.astype(float))
miss = np.abs((lon - ref + np.pi) % (2 * np.pi) - np.pi)
bound = 4 * ulp(abs(body['longitude_epoch']) + np.abs(body['longitude_rate'] * since.astype(float)))
print(f'{len(since)} dates bit for bit; the float64 series within {np.max(miss / bound) * 4:.2f} ulp, '
f'worst {np.degrees(np.max(miss)):.3f} degrees')
assert np.all(miss <= bound)
sys.path.insert(0, os.path.join(root, 'kernels'))
from lusine.kernels.sky.moon_longitude import moon_longitude
lon, _, since, run = moon(CENTURY)
body = rows.row(run.unit, 'lusine.bodies', 'moon')
TURN = F32(2.0 * np.pi)
def within_turn(x):
return F32(x - F32(TURN * F32(np.floor(F32(x / TURN)))))
mean = within_turn(F32(F32(body['longitude_epoch']) + F32(F32(body['longitude_rate']) * since)))
model = within_turn(F32(mean + tree([term(r, since) for r in terms(run, 'longitude')])))
assert np.array_equal(F32(lon), model)
ref = moon_longitude(since.astype(float))
miss = np.abs((lon - ref + np.pi) % (2 * np.pi) - np.pi)
bound = 4 * ulp(abs(body['longitude_epoch']) + np.abs(body['longitude_rate'] * since.astype(float)))
print(f'{len(since)} dates bit for bit; the float64 series within {np.max(miss / bound) * 4:.2f} ulp, '
f'worst {np.degrees(np.max(miss)):.3f} degrees')
assert np.all(miss <= bound)
4001 dates bit for bit; the float64 series within 1.75 ulp, worst 0.059 degrees
Invariant: round the sky in a month, within a turn¶
Over 2024 the longitude advances about 13.2° a day on average, a turn in the Moon's month, every value within a turn, [0, 2π).
In [6]:
Copied!
days = 2460310.5 + 0.5 * np.arange(731)
lon, _, _, _ = moon(days)
rate = float(np.degrees(np.mean(np.diff(np.unwrap(lon))))) / 0.5
print(f'{rate:.3f} degrees a day')
assert abs(rate - 13.18) < 0.05 and np.all((lon >= 0.0) & (lon < 2 * np.pi))
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(days[:120] - days[0], np.degrees(lon[:120]), '.', ms=3)
ax.set_xlabel('days from 1 January 2024')
ax.set_ylabel('ecliptic longitude (°)')
plt.show()
days = 2460310.5 + 0.5 * np.arange(731)
lon, _, _, _ = moon(days)
rate = float(np.degrees(np.mean(np.diff(np.unwrap(lon))))) / 0.5
print(f'{rate:.3f} degrees a day')
assert abs(rate - 13.18) < 0.05 and np.all((lon >= 0.0) & (lon < 2 * np.pi))
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(days[:120] - days[0], np.degrees(lon[:120]), '.', ms=3)
ax.set_xlabel('days from 1 January 2024')
ax.set_ylabel('ecliptic longitude (°)')
plt.show()
13.163 degrees a day
