The moon's phase¶
The moon's illuminated fraction, 0 new to 1 full, from its elongation from the sun, cos psi = cos(beta) cos(lambda - lambda_sun), with the phase angle taken as the elongation's supplement: k = (1 - cos psi) / 2, good to a few thousandths.
Its inputs are the moon's ecliptic latitude and longitude and the sun's longitude. Source: Meeus, Astronomical Algorithms, ch. 48. Checked against a reference in numpy, which has its own tests.
- Source
- Meeus, Astronomical Algorithms, ch. 48.
- Proof
- curve, reference, invariant.
moon_phase.trait
namespace lusine.sky;
float moon_phase(float<rad> moon_latitude, float<rad> moon_longitude, float<rad> sun_longitude)
{
return (1.0 - cos(moon_latitude) * cos(moon_longitude - sun_longitude)) / 2.0;
}
The moon's illuminated fraction from its elongation from the sun, k = (1 - cos psi) / 2 with cos psi = cos(beta) cos(lambda - lambda_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 trait_run
EPS = float(np.finfo(np.float32).eps)
The kernel¶
The law as the library states it, and as this page runs it.
print(trait_run.source('sky/moon_phase.trait'))
// lusine.sky.moon_phase: the moon's illuminated fraction, 0 new to 1 full, from its elongation from the sun, cos psi =
// cos(beta) cos(lambda - lambda_sun), with the phase angle taken as the elongation's supplement: k = (1 - cos psi) / 2,
// good to a few thousandths. Its inputs are the moon's ecliptic latitude and longitude and the sun's longitude.
// Source: Meeus, Astronomical Algorithms, ch. 48.
// Checked against a reference in numpy, which has its own tests.
namespace lusine.sky;
float moon_phase(float<rad> moon_latitude, float<rad> moon_longitude, float<rad> sun_longitude)
{
return (1.0 - cos(moon_latitude) * cos(moon_longitude - sun_longitude)) / 2.0;
}
The program¶
Each section runs the kernel over a lattice of items, one a case: the host sets every parameter, and main() binds the kernel over the items.
print(trait_run.law_program('sky/moon_phase.trait', 'moon_phase'))
namespace probe;
lusine.version version = { { number = 1 } };
struct item {
host float<rad> moon_latitude;
host float<rad> moon_longitude;
host float<rad> sun_longitude;
float result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { items.result_0 = lusine.sky.moon_phase(items.moon_latitude, items.moon_longitude, items.sun_longitude); }
Reference: Meeus's example 48.a¶
On 1992 April 12 at 0h TD the moon stands at ecliptic longitude 133.162655° and latitude -3.229126° (Meeus 47.a), the sun at right ascension 20.6579° and declination 8.6964°, and Meeus's full formula, with the distances, gives k = 0.6786. The sun's ecliptic longitude is its equatorial place turned by the obliquity, 23.440636°, in numpy; taking the phase angle as the elongation's supplement, as this law does, leaves it within a few thousandths.
d = np.pi / 180
alpha, delta, e = 20.6579 * d, 8.6964 * d, 23.440636 * d
sun = np.arctan2(np.sin(alpha) * np.cos(e) + np.tan(delta) * np.sin(e), np.cos(alpha))
k = trait_run.law('sky/moon_phase.trait', 'moon_phase', moon_latitude=-3.229126 * d, moon_longitude=133.162655 * d,
sun_longitude=sun)[0]
print(f'k = {k:.4f} against Meeus 0.6786')
assert abs(k - 0.6786) < 0.002
k = 0.6775 against Meeus 0.6786
Curve: a month of elongations¶
The fraction over a whole turn of the elongation, for a moon on the ecliptic and one 5° from it, against (1 - cos(beta) cos(psi)) / 2 in double precision from the same binary32 inputs, within binary32's rounding of the library's cosines (8 EPS).
psi = np.linspace(-np.pi, np.pi, 73)
lat = np.repeat([0.0, 5 * np.pi / 180], len(psi))
lon = np.tile(psi + 1.0, 2)
k = trait_run.law('sky/moon_phase.trait', 'moon_phase', moon_latitude=lat, moon_longitude=lon, sun_longitude=1.0)
f = lambda x: np.float32(x).astype(float)
want = (1 - np.cos(f(lat)) * np.cos(f(lon) - f(1.0))) / 2
assert np.max(np.abs(k - want)) <= 8 * EPS, np.max(np.abs(k - want)) / EPS
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(np.degrees(psi), want[:len(psi)], label='(1 - cos psi) / 2')
ax.plot(np.degrees(psi), k[:len(psi)], 'o', ms=3, label='moon_phase, on the ecliptic')
ax.plot(np.degrees(psi), k[len(psi):], 'x', ms=3, label='moon_phase, 5° off it')
ax.set_xlabel('elongation, moon less sun (deg)')
ax.set_ylabel('illuminated fraction')
ax.legend()
plt.show()
Limit: new, quarter, full, and the eclipse season¶
A moon with the sun on the ecliptic is dark, exactly; a quarter of the sky away it is half lit, and opposite it full, to binary32's rounding of the cosine at a right angle and a half turn. A new moon 5° off the ecliptic, as most are, keeps a sliver, (1 - cos 5°) / 2.
d = np.pi / 180
k = trait_run.law('sky/moon_phase.trait', 'moon_phase', moon_latitude=[0.0, 0.0, 0.0, 0.0, 5 * d],
moon_longitude=[2.0, 2.0 + np.pi / 2, 2.0 - np.pi / 2, 2.0 + np.pi, 2.0], sun_longitude=2.0)
assert k[0] == 0.0, k
assert abs(k[1] - 0.5) <= 2 * EPS and abs(k[2] - 0.5) <= 2 * EPS and abs(k[3] - 1.0) <= 2 * EPS, k
assert abs(k[4] - (1 - np.cos(float(np.float32(5 * d)))) / 2) <= 8 * EPS, k[4]
print('new', k[0], 'quarters', k[1], k[2], 'full', k[3], '5° off at new', k[4])
new 0.0 quarters 0.4999999701976776 0.5 full 1.0 5° off at new 0.0019026398658752441
Invariant: waxing as waning, and between none and all¶
The phase is the elongation's cosine, so a moon as far east of the sun as another is west of it shows the same fraction, bit for bit; the moon's latitude's sign does not matter either; and no fraction leaves [0, 1].
rng = np.random.default_rng(5)
m, s, b = rng.uniform(0, 2 * np.pi, 200), rng.uniform(0, 2 * np.pi, 200), rng.uniform(-0.09, 0.09, 200)
k = trait_run.law('sky/moon_phase.trait', 'moon_phase', moon_latitude=np.concatenate([b, -b]),
moon_longitude=np.concatenate([m, s]), sun_longitude=np.concatenate([s, m]))
assert np.array_equal(k[:200], k[200:])
assert k.min() >= 0.0 and k.max() <= 1.0
print(f'200 places mirrored, the same; fractions {k.min():.4f} to {k.max():.4f}')
200 places mirrored, the same; fractions 0.0006 to 0.9996