Sidereal time¶
The local mean sidereal time, Greenwich's mean sidereal time (IAU 1982) as a quadratic in the time since the Earth's epoch, theta0 + theta1 t + theta2 t^2, taken within a turn, and turned by the longitude east.
The three terms are the Earth's (its row of lusine.bodies), passed in; the time is given since the epoch they are taken at, since a binary32 Julian day cannot hold the hours. PRECISION: binary32 seconds since J2000 cannot hold the minute either. In 2026 (t near 8.4e8 s) the time's quantum is 64 s, 0.27 deg of sidereal time, and the product sidereal_rate * since_epoch (near 61,600 rad) rounds to 0.22 deg; the reference's notebook bounds the whole at 4 EPS |sidereal_rate t| + 4 EPS 2 pi, about 1.7 deg, seven minutes of sidereal time, growing with the time from the epoch. Good for the look of a sky, not for an ephemeris; a game that needs the minute hands the kernel a day count and a time of day, each held exactly. Source: Meeus, Astronomical Algorithms, eq. 12.4. Checked against a reference in numpy, which has its own tests.
- Source
- Meeus, Astronomical Algorithms, eq. 12.4 (IAU 1982).
- Proof
- reference, invariant.
sidereal.trait
namespace lusine.sky;
float<rad> sidereal(float<s> since_epoch, float<rad> longitude, float<rad> sidereal_epoch, float<rad/s> sidereal_rate,
float<rad/s2> sidereal_quadratic)
{
let gmst = sidereal_epoch + sidereal_rate * since_epoch + sidereal_quadratic * since_epoch * since_epoch;
return gmst - 6.283185307179586 * floor(gmst / 6.283185307179586) + longitude;
}
The local mean sidereal time, Greenwich's (IAU 1982) as theta0 + theta1 t + theta2 t^2 in the time since the Earth's epoch, within a turn, turned by the longitude east. Every number it runs on is the Earth's, passed in as a game binds its row: the kernel declares none but the turn, 2 pi. It takes the time since the epoch, not a Julian day: a binary32 Julian day is a quarter of a day apart, and a sidereal time from it would be off by a quarter of the sky.
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/sidereal.trait'))
// lusine.sky.sidereal: the local mean sidereal time, Greenwich's mean sidereal time (IAU 1982) as a quadratic in the
// time since the Earth's epoch, theta0 + theta1 t + theta2 t^2, taken within a turn, and turned by the longitude
// east. The three terms are the Earth's (its row of lusine.bodies), passed in; the time is given since the epoch they
// are taken at, since a binary32 Julian day cannot hold the hours.
// PRECISION: binary32 seconds since J2000 cannot hold the minute either. In 2026 (t near 8.4e8 s)
// the time's quantum is 64 s, 0.27 deg of sidereal time, and the product sidereal_rate * since_epoch (near 61,600 rad)
// rounds to 0.22 deg; the reference's notebook bounds the whole at 4 EPS |sidereal_rate t| + 4 EPS 2 pi, about 1.7
// deg, seven minutes of sidereal time, growing with the time from the epoch. Good for the look of a sky, not for an
// ephemeris; a game that needs the minute hands the kernel a day count and a time of day, each held exactly.
// Source: Meeus, Astronomical Algorithms, eq. 12.4.
// Checked against a reference in numpy, which has its own tests.
namespace lusine.sky;
float<rad> sidereal(float<s> since_epoch, float<rad> longitude, float<rad> sidereal_epoch, float<rad/s> sidereal_rate,
float<rad/s2> sidereal_quadratic)
{
let gmst = sidereal_epoch + sidereal_rate * since_epoch + sidereal_quadratic * since_epoch * since_epoch;
return gmst - 6.283185307179586 * floor + longitude;
}
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. The Earth's terms are its row of lusine.bodies, earth, in its own units: the epoch J2000.0 (JD 2451545.0), theta0 = 280.46061837°, theta1 = 360.98564736629° a day and theta2 = 0.000387933° a Julian century squared; the page turns them into radians and seconds, as a game binding the row does.
print(trait_run.law_program('sky/sidereal.trait', 'sidereal'))
d, DAY = np.pi / 180, 86400.0
EARTH = dict(sidereal_epoch=280.46061837 * d, sidereal_rate=360.98564736629 * d / DAY,
sidereal_quadratic=0.000387933 * d / (36525 * DAY) ** 2)
J2000 = 2451545.0
namespace probe;
lusine.version version = { { number = 1 } };
struct item {
host float<s> since_epoch;
host float<rad> longitude;
host float<rad> sidereal_epoch;
host float<rad/s> sidereal_rate;
host float<rad/s2> sidereal_quadratic;
float<rad> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { items.result_0 = lusine.sky.sidereal(items.since_epoch, items.longitude, items.sidereal_epoch, items.sidereal_rate, items.sidereal_quadratic); }
Reference: J2000.0¶
At 2000 January 1, 12:00 the Greenwich mean sidereal time is 280.46061837°, 18h 41m 50.5s (Meeus 12.4): theta0 itself, exactly as binary32 holds it.
g = trait_run.law('sky/sidereal.trait', 'sidereal', since_epoch=0.0, longitude=0.0, **EARTH)[0]
print(f'{np.degrees(g):.6f}°')
assert g == float(np.float32(280.46061837 * d))
280.460627°
Reference: Meeus's examples 12.a and 12.b¶
On 1987 April 10 at 0h UT Greenwich's mean sidereal time is 13h 10m 46.3668s, and at 19h 21m 0s it is 8h 34m 57.0896s. Thirteen years from the epoch theta1 t is near 30,000 rad, and binary32 holds it to a few of its own roundings: 4 EPS |theta1 t|, about 0.8°.
jd = np.array([2446895.5, 2446896.30625])
want = np.array([13 + 10 / 60 + 46.3668 / 3600, 8 + 34 / 60 + 57.0896 / 3600]) * 15 * d
t = (jd - J2000) * DAY
g = trait_run.law('sky/sidereal.trait', 'sidereal', since_epoch=t, longitude=0.0, **EARTH)
err = (g - want + np.pi) % (2 * np.pi) - np.pi
bound = 4 * EPS * EARTH['sidereal_rate'] * np.abs(t)
print('off by', np.degrees(err), 'deg, within', np.degrees(bound))
assert np.all(np.abs(err) <= bound)
off by [-0.06751801 -0.04618107] deg, within [0.80032481 0.80018603]
Curve: half a century of binary32 time¶
From 1970 to 2050, against the same quadratic in double precision at the binary32 time the kernel was given, taken within a turn. The difference is binary32's rounding of theta1 t, growing with the time from the epoch, and held within 4 EPS |theta1 t| + 4 EPS 2 pi.
t = np.linspace(-30, 50, 161) * 365.25 * DAY
g = trait_run.law('sky/sidereal.trait', 'sidereal', since_epoch=t, longitude=0.0, **EARTH)
t32 = np.float32(t).astype(float)
want = (EARTH['sidereal_epoch'] + EARTH['sidereal_rate'] * t32 + EARTH['sidereal_quadratic'] * t32 ** 2) % (2 * np.pi)
err = (g - want + np.pi) % (2 * np.pi) - np.pi
bound = 4 * EPS * EARTH['sidereal_rate'] * np.abs(t32) + 4 * EPS * 2 * np.pi
assert np.all(np.abs(err) <= bound) and g.min() >= 0.0 and g.max() <= 2 * np.pi
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(2000 + t / (365.25 * DAY), np.degrees(err), 'o', ms=3, label='sidereal - quadratic')
ax.plot(2000 + t / (365.25 * DAY), np.degrees(bound), lw=0.8, label='4 EPS |theta1 t|')
ax.plot(2000 + t / (365.25 * DAY), -np.degrees(bound), lw=0.8, color='C1')
ax.set_xlabel('year')
ax.set_ylabel('deg')
ax.legend()
plt.show()
Invariant: a sidereal day, and the longitude added¶
A solar day after the epoch the sky has turned 360.98564736629°, 0.98564736629° beyond a whole turn, and a place east of Greenwich is ahead by its longitude, west behind, to binary32's rounding of an angle within a turn.
lon = np.radians([0.0, 15.0, -74.0, 151.2])
g = trait_run.law('sky/sidereal.trait', 'sidereal', since_epoch=np.repeat([0.0, DAY], 4), longitude=np.tile(lon, 2),
**EARTH)
turn = g[4] - g[0]
print(f'{np.degrees(turn):.7f}° a day beyond a whole turn')
assert abs(turn - 0.98564736629 * d) <= 4 * EPS * 2 * np.pi
for day in (0, 4):
assert np.all(np.abs(g[day:day + 4] - g[day] - np.float32(lon)) <= 2 * EPS * (2 * np.pi + np.abs(lon))), g
0.9856509° a day beyond a whole turn