skywatch algorithms

This document describes every computation skywatch performs, the worked examples it is validated against, and the accuracy you can expect. Chapter numbers refer to Jean Meeus, Astronomical Algorithms, 2nd edition (Willmann-Bell, 1998).

1. Time

Julian Day (ch. 7, eq. 7.1). For Gregorian dates with Y, M (Jan/Feb treated as months 13/14 of the previous year), A = ⌊Y/100⌋, B = 2 − A + ⌊A/4⌋:

JD = ⌊365.25 (Y + 4716)⌋ + ⌊30.6001 (M + 1)⌋ + D + B − 1524.5

Example 7.a: 1957 Oct 4.81 → JD 2436116.31 ✔.

ΔT = TT − UT uses the Espenak & Meeus (2006) piecewise polynomials, e.g. for 2005–2050 ΔT = 62.92 + 0.32217 t + 0.005589 t², t = y − 2000. All ephemerides take TT; observer-facing times are UT converted to the location's IANA zone with zoneinfo.

Sidereal time (ch. 12, eq. 12.4):

θ₀ = 280.46061837° + 360.98564736629° (JD − 2451545) + 0.000387933 T² − T³/38 710 000

Apparent sidereal time adds the equation of the equinoxes Δψ cos ε. Example 12.a: 1987 Apr 10 0h UT → 13h10m46.3668s ✔; 12.b: 19h21m UT → 128.7378734° ✔.

2. Coordinates

• Nutation (ch. 22, low-precision form, ±0.5″): Δψ = −17.20″ sin Ω − 1.32″ sin 2L − 0.23″ sin 2L′ + 0.21″ sin 2Ω, Δε = 9.20″ cos Ω + 0.57″ cos 2L + 0.10″ cos 2L′ − 0.09″ cos 2Ω. Example 22.a: Δψ = −3.788″, Δε = 9.443″ (within 0.5″) ✔.
• Precession (ch. 21, eqs. 21.2–21.4): rigorous ζ, z, θ between arbitrary epochs. Example 21.b (θ Persei to 2028 Nov 13.19): α = 41.547214°, δ = 49.348483° reproduced to < 0.1″ ✔.
• Nutation in RA/Dec (eq. 23.1) converts mean-of-date to apparent positions for catalog objects.
• Ecliptic ↔ equatorial ↔ horizontal (ch. 13). Azimuth is reported from north through east. Example 13.b (Venus from Washington): h = 15.1249°, A = 68.0337° (south-based) ✔ to 0.001°.
• Refraction — Bennett's formula (ch. 16, eq. 16.3) for an apparent altitude h₀: R = 1′ / tan(h₀ + 7.31/(h₀ + 4.4)). skywatch needs the inverse (true → apparent), which it obtains by fixed-point iteration h ← h_true + R(h) (six iterations converge to 1e-6°). R is scaled by exp(−elevation / 8434 m) for site pressure. Example 16.a: h₀ = 0.5° → R = 28.754′ ✔.
• Angular separation uses the Vincenty atan2 form, which is stable for both 0.1° conjunctions and near-180° oppositions.

3. The Sun (ch. 25, low precision)

L₀ = 280.46646 + 36000.76983 T + 0.0003032 T²
M  = 357.52911 + 35999.05029 T − 0.0001537 T²
C  = (1.914602 − 0.004817 T − 0.000014 T²) sin M + (0.019993 − 0.000101 T) sin 2M + 0.000289 sin 3M
λ  = L₀ + C − 0.00569 − 0.00478 sin Ω,  ε = ε₀ + 0.00256 cos Ω

Example 25.a (1992 Oct 13.0 TD): λ = 199.90895°, α = 198.38083°, δ = −7.78507°, R = 0.99766 au — all reproduced to 1e-5° ✔. Equation of time uses Smart's series (eq. 28.3): Example 28.b gives 13m42.7s ✔.

4. The Moon

Position (ch. 47). The fundamental arguments L′, D, M, M′, F are the full polynomials. skywatch sums the 60 terms of Table 47.A (Σl, Σr) and the 30 largest terms of Table 47.B (Σb), multiplying terms in M by E = 1 − 0.002516 T (and E² for 2M), plus the additive A₁, A₂, A₃ corrections. Then λ = L′ + Σl·10⁻⁶, β = Σb·10⁻⁶, Δ = 385000.56 + Σr/1000 km, nutation is added to λ.

Example 47.a (1992 Apr 12.0 TD):

Illumination (ch. 48): phase angle tan i = R sin ψ / (Δ − R cos ψ) with ψ the Sun–Moon elongation, k = (1 + cos i)/2. Example 48.a: k = 0.6786 ✔. Phase names use k and waxing/waning (new < 3 %, quarter 45–55 %, full > 97 %).

Phases (ch. 49): mean phase from lunation number k, then the 25-term periodic corrections for new/full or quarter phases, the W correction for quarters, and the 14 planetary corrections A₁…A₁₄. Example 49.a: new moon of 1977 Feb → JDE 2443192.65118 ✔ (< 1 s). Converting JDE (TT) to UT subtracts ΔT.

5. Planets

Heliocentric positions use JPL's Approximate Positions of the Planets, Table 1 (Standish), valid 1800–2050: six elements a, e, I, L, ϖ, Ω plus linear rates per Julian century. For each planet:

1. ω = ϖ − Ω, M = L − ϖ; solve Kepler's equation E − e sin E = M by Newton iteration (1e-12 rad).
2. Orbital-plane coordinates x′ = a(cos E − e), y′ = a√(1−e²) sin E, rotate by ω, I, Ω to J2000 ecliptic.
3. Geocentric vector = planet(t − τ) − Earth-Moon barycentre(t), iterating the light-time τ = Δ · 0.0057755 d/au.
4. Rotate to J2000 equatorial (ε = 23.4392911°), precess to date and apply nutation.

Example 33.a (Venus 1992 Dec 20.0 TD; Meeus uses full VSOP87): α = 21h04m41.454s, δ = −18°53′16.84″, Δ = 0.910845 au. skywatch gives 21h04m42.5s, −18°53′11″, 0.91100 au — 15″ and 6″ off, comfortably within the JPL table's stated error.

Magnitudes: simplified Mallama & Hilton (2018) phase curves, V = 5 log₁₀(rΔ) + f(i), e.g. Mars −1.601 + 0.02267 i − 0.0001302 i²; Saturn includes the ring term −1.825 sin B where B is the sub-Earth latitude on the ring plane computed from Saturn's pole (α 40.589°, δ 83.537°). Apparent diameter = 2·s₀/Δ with Meeus's semidiameters at 1 au.

6. Rise, transit, set and twilight

Instead of Meeus's three-point interpolation (ch. 15), skywatch evaluates the exact geometric altitude minus the standard altitude h₀ on a 10-minute grid and refines every sign change by bisection to 1 s. This works for any object (including the fast-moving Moon) and naturally reports multiple or zero events.

Transit is the zero crossing of the hour angle from negative to positive. If no crossing exists in the window the object is classified as circumpolar or never-rising from the altitude at the window's midpoint.

Example 15.a (Venus at Boston, 1988 Mar 20, using Meeus's tabulated positions with quadratic interpolation): rise 12h25m26s, transit 19h40m18s, set 2h54m40s UT — reproduced within 1 minute ✔.

Sunrise/sunset reference checks (tests, ±2 min): London 2024-06-21 and 2024-12-21 against timeanddate.com, Pune 2024-01-15 and 2024-06-21 against api.sunrise-sunset.org. Measured differences are 0–1.3 min.

7. Constellation identification

The IAU boundaries were drawn along B1875.0 meridians and parallels. Roman (1987, PASP 99, 695; CDS VI/42) tabulated them as 357 rows (RA_low, RA_high, Dec_low, constellation) sorted by descending Dec_low. A position is precessed to B1875.0 (JD 2405889.2586); the answer is the first row with Dec ≥ Dec_low and RA_low ≤ RA < RA_high. All eight test points from Roman's paper pass.

8. Planning: the dark window and target score

Night = local noon to next local noon. Dark window = astronomical night (Sun < −18°) intersected with "Moon below its rise/set altitude", sampled every 5 minutes. If it is empty, skywatch falls back to astronomical night with the Moon up, then nautical twilight, civil twilight and finally sunset→sunrise.

Every candidate (160 deep-sky objects, the 1st-magnitude named stars and the planets) is evaluated every 15 minutes across the window; its best sample gives the score

score = 100 · sin(alt) · (0.55·B + 0.45·S) · (1 − 0.8·P)     (0 if alt < 15°)

• B brightness: deep-sky clip((11 − m)/9, 0.05, 1); stars/planets clip((6 − m)/8, 0.1, 1).
• S size: deep-sky clip(log₁₀(size′ + 1)/log₁₀ 121, 0.1, 1) (saturates at 2°); planets 0.6, stars 0.4.
• P moonlight: if the Moon is up, k · max(0.2, 1 − sep/120°) where k is its illuminated fraction; multiplied by 0.25 for point sources, which suffer far less from sky brightness than nebulae.

sin(alt) approximates the reduction in airmass (extinction and seeing) toward the horizon.

Worked example — M31 from Pune at 23:41 IST on 2026-10-16 (Moon set): alt 67° → sin = 0.921; m = 3.4 → B = 0.844; size 177.8′ → S = 1.0; P = 0. Score = 100 · 0.921 · (0.464 + 0.45) = 84 ✔ (row 2 of the tonight example in the README).

9. Monthly events

• Phases: ch. 49 for every k in range.
• Conjunctions: all pairs of Moon + 7 planets are sampled every 2 h; local minima of separation are refined by golden-section search (1 min) and kept when ≤ 3°. The 2020 Dec 21 Jupiter–Saturn conjunction is found at 0.10° (true 0.102°).
• Elongations/oppositions: local maxima of Sun–planet elongation (12 h sampling, golden-section refinement). Mercury/Venus maxima are greatest elongations (east = evening if the planet's RA exceeds the Sun's); superior-planet maxima above 150° are oppositions. Mars 2020 opposition: 2020-10-13 (true 23:20 UT) ✔.
• Meteor showers: peak defined by the IMO solar longitude (J2000); the apparent solar longitude is converted to J2000 by subtracting general precession (1.397° per century) and the crossing is found by bisection.

10. Accuracy summary

Sources

• J. Meeus, Astronomical Algorithms, 2nd ed., 1998.
• E. M. Standish, Keplerian Elements for Approximate Positions of the Major Planets, JPL/SSD.
• F. Espenak & J. Meeus, Five Millennium Canon of Solar Eclipses, NASA TP-2006-214141 (ΔT polynomials).
• A. Mallama & J. Hilton, "Computing apparent planetary magnitudes for The Astronomical Almanac", Astronomy and Computing 25 (2018).
• N. Roman, "Identification of a constellation from a position", PASP 99, 695 (1987); CDS catalog VI/42.
• G. G. Bennett, "The calculation of astronomical refraction in marine navigation", J. Navigation 35 (1982).
• Catalog data: HYG database v3.8 (stars), OpenNGC (deep-sky), d3-celestial (constellation figures), IMO Meteor Shower Calendar (shower peaks).