The Moon's latitude¶

The Moon's geocentric ecliptic latitude, the sum of its periodic terms in latitude, beta = sum A sin(phase + rate t): the Moon's rows of lusine.body_terms whose coordinate is the latitude, each one lusine.sky.periodic_term, folded in one kernel.

The time is given since the Moon's epoch (its row of lusine.bodies, 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_latitude.trait

namespace lusine.sky;

float<rad> moon_latitude(float<s> since_epoch)
{
  return sum(lusine.sky.periodic_term(since_epoch, lusine.body_terms.amplitude, lusine.body_terms.phase,
                                      lusine.body_terms.rate))
         where lusine.body_terms.body == lusine.bodies["moon"]
               && lusine.body_terms.coordinate == lusine.coordinates["latitude"];
}

The Moon's geocentric ecliptic latitude, the sum of its four largest periodic terms in latitude, the Almanac's low-precision series. One kernel sums them: a fold over the Moon's rows of lusine.body_terms whose coordinate is the latitude (a row of lusine.coordinates), applying lusine.sky.periodic_term to each. Every number is a row's; the kernel declares none.

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 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]:

print(trait_run.source('sky/moon_latitude.trait'))

print(trait_run.source('sky/moon_latitude.trait'))

// lusine.sky.moon_latitude: the Moon's geocentric ecliptic latitude, the sum of its periodic terms in latitude,
// beta = sum A sin(phase + rate t): the Moon's rows of lusine.body_terms whose coordinate is the latitude, each one
// lusine.sky.periodic_term, folded in one kernel. The time is given since the Moon's epoch (its row of lusine.bodies,
// 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_latitude(float<s> since_epoch)
{
  return sum(lusine.sky.periodic_term(since_epoch, lusine.body_terms.amplitude, lusine.body_terms.phase,
                                      lusine.body_terms.rate))
         where lusine.body_terms.body == lusine.bodies["moon"]
               && lusine.body_terms.coordinate == lusine.coordinates["latitude"];
}

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]:

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.25°: the series' few tenths of a degree, the four terms it keeps. At the eclipse of 8 April 2024 the Moon was on the ecliptic: astronomy-engine's latitude there, 0.35°, within 0.1°.

In [4]:

_, lat, _, _ = moon(DATES)
_, want = published(DATES)
miss = np.degrees(lat) - want
print(f'worst {np.max(np.abs(miss)):.3f} degrees over {len(DATES)} dates')
assert np.max(np.abs(miss)) < 0.25
_, at, _, _ = moon([ECLIPSE])
_, there = published([ECLIPSE])
print(f'8 April 2024, 18:17 UT: {np.degrees(at[0]):+.3f} degrees; astronomy-engine {there[0]:+.3f}')
assert abs(np.degrees(at[0]) - there[0]) < 0.1

_, lat, _, _ = moon(DATES) _, want = published(DATES) miss = np.degrees(lat) - want print(f'worst {np.max(np.abs(miss)):.3f} degrees over {len(DATES)} dates') assert np.max(np.abs(miss)) < 0.25 _, at, _, _ = moon([ECLIPSE]) _, there = published([ECLIPSE]) print(f'8 April 2024, 18:17 UT: {np.degrees(at[0]):+.3f} degrees; astronomy-engine {there[0]:+.3f}') assert abs(np.degrees(at[0]) - there[0]) < 0.1

worst 0.180 degrees over 5923 dates
8 April 2024, 18:17 UT: +0.397 degrees; astronomy-engine +0.346

Reference: bit for bit, and the float64 series within binary32's rounding¶

From 1900 to 2100 the kernel's latitude is a binary32 model's bit for bit: the four terms, each rounded once an operation with the library's sin, summed in Trait's tree, from the rows as the program holds them. The numpy reference, the same series in float64 at the same binary32 time, agrees within Σ|A| 2 ulp(|phase| + |rate t|) + 1e-6 rad: each term's argument rounds twice and its rate once (two ulp of it), times its amplitude; the sin, the products and the sums, on values under 0.1 rad, add a few ulp of 1e-8.

In [5]:

sys.path.insert(0, os.path.join(root, 'kernels'))
from lusine.kernels.sky.moon_latitude import moon_latitude
_, lat, since, run = moon(CENTURY)
rows_ = terms(run, 'latitude')
model = tree([term(r, since) for r in rows_])
assert np.array_equal(F32(lat), model)
miss = np.abs(lat - moon_latitude(since.astype(float)))
bound = sum(abs(r['amplitude']) * 2 * ulp(abs(r['phase']) + np.abs(r['rate'] * since.astype(float)))
            for r in rows_) + 1e-6
print(f'{len(since)} dates bit for bit; the float64 series within {np.max(miss / bound):.2f} of the bound, '
      f'worst {np.degrees(np.max(miss)):.4f} degrees')
assert np.all(miss <= bound)

sys.path.insert(0, os.path.join(root, 'kernels')) from lusine.kernels.sky.moon_latitude import moon_latitude _, lat, since, run = moon(CENTURY) rows_ = terms(run, 'latitude') model = tree([term(r, since) for r in rows_]) assert np.array_equal(F32(lat), model) miss = np.abs(lat - moon_latitude(since.astype(float))) bound = sum(abs(r['amplitude']) * 2 * ulp(abs(r['phase']) + np.abs(r['rate'] * since.astype(float))) for r in rows_) + 1e-6 print(f'{len(since)} dates bit for bit; the float64 series within {np.max(miss / bound):.2f} of the bound, ' f'worst {np.degrees(np.max(miss)):.4f} degrees') assert np.all(miss <= bound)

4001 dates bit for bit; the float64 series within 0.41 of the bound, worst 0.0037 degrees

Invariant: within the orbit's band¶

Over sixty years the latitude stays within the orbit's 5.1° inclination and its largest terms' reach, 5.35°, and reaches both edges, crossing the ecliptic twice a month.

In [6]:

_, lat, _, _ = moon(DATES)
deg = np.degrees(lat)
print(f'from {deg.min():.2f} to {deg.max():.2f} degrees')
assert np.max(np.abs(deg)) < 5.35 and deg.max() > 5.0 and deg.min() < -5.0
days = 2460310.5 + 0.5 * np.arange(731)
_, lat, _, _ = moon(days)
crossings = int(np.sum(np.diff(np.sign(lat)) != 0))
print(f'{crossings} crossings of the ecliptic in 2024')
assert 24 <= crossings <= 28
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(days - days[0], np.degrees(lat))
ax.set_xlabel('days from 1 January 2024')
ax.set_ylabel('ecliptic latitude (°)')
plt.show()

_, lat, _, _ = moon(DATES) deg = np.degrees(lat) print(f'from {deg.min():.2f} to {deg.max():.2f} degrees') assert np.max(np.abs(deg)) < 5.35 and deg.max() > 5.0 and deg.min() < -5.0 days = 2460310.5 + 0.5 * np.arange(731) _, lat, _, _ = moon(days) crossings = int(np.sum(np.diff(np.sign(lat)) != 0)) print(f'{crossings} crossings of the ecliptic in 2024') assert 24 <= crossings <= 28 fig, ax = plt.subplots(figsize=(7, 3)) ax.plot(days - days[0], np.degrees(lat)) ax.set_xlabel('days from 1 January 2024') ax.set_ylabel('ecliptic latitude (°)') plt.show()

from -5.33 to 5.32 degrees
27 crossings of the ecliptic in 2024