Something went wrong. Try again.
A Zig library for astronomical calculations: the Sun and Moon from anywhere on Earth, the Moon's phases and libration, and time scales.
Something went wrong. Try again.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259// SPDX-FileCopyrightText: © 2026 Jeffrey C. Ollie <jeff@ocjtech.us>// SPDX-License-Identifier: MIT
//! An observer on the Earth's surface, and what a body looks like from//! there: its topocentric position, shifted from the geocentric one by its//! parallax (Meeus, *Astronomical Algorithms*, 2nd ed., chs. 11 and 40),//! and the refraction that lifts it from its airless altitude to where it//! is seen (ch. 16).//!//! Longitude counts positive east of Greenwich, as the IAU, WGS 84 and JPL//! Horizons have it. Meeus counts it positive *west*; his longitudes change//! sign here (see ADR 0003).//!//! Left out, each well inside the accuracy of the Sun's and the Moon's//! positions: diurnal aberration, at most 0.32″; polar motion, a few tenths//! of a second of arc; the difference between UTC and UT1, under 0.9 s of//! the Earth's rotation; and the deflection of the vertical, by which the//! local plumb line departs from the ellipsoid's normal by up to about a//! minute of arc in mountains.
const std = @import("std");const math = std.math;
const Angle = @import("Angle.zig");const coord = @import("coord.zig");const earth = @import("earth.zig");const Instant = @import("time.zig").Instant;
/// A point on the Earth's surface.pub const Observer = struct { /// Geodetic latitude, the latitude on maps: the angle between the /// equator and the normal to the WGS 84 ellipsoid, positive north. It /// must be within [−90°, +90°]; beyond that, every result is NaN. latitude: Angle, /// Longitude, positive east of Greenwich. Any finite value is taken /// modulo 360°. longitude: Angle, /// Height above the WGS 84 ellipsoid, in meters. Height above sea level /// differs from it by the geoid, at most about 100 m, which moves the /// Moon by about 0.003″. Any finite height is accepted, including /// negative ones. height_m: f64 = 0,
/// The observer's place in the Earth's meridian plane: ρ sin φ′ and /// ρ cos φ′, where ρ is the distance from the Earth's center in units of /// the equatorial radius and φ′ the geocentric latitude (Meeus, p. 82, /// with the WGS 84 ellipsoid in place of the IAU 1976 one). NaN if the /// latitude is outside [−90°, +90°] or anything is not finite. pub fn geocentric(self: Observer) Geocentric { if (!self.isValid()) return .{ .rho_sin_phi = math.nan(f64), .rho_cos_phi = math.nan(f64) }; const phi = self.latitude.radians(); const u = math.atan2(polar_ratio * @sin(phi), @cos(phi)); const h = self.height_m / (equatorial_radius_km * 1000); return .{ .rho_sin_phi = polar_ratio * @sin(u) + h * @sin(phi), .rho_cos_phi = @cos(u) + h * @cos(phi), }; }
fn isValid(self: Observer) bool { return @abs(self.latitude.radians()) <= math.pi / 2.0 and math.isFinite(self.longitude.radians()) and math.isFinite(self.height_m); }};
/// An observer's place in the Earth's meridian plane; see/// `Observer.geocentric`.pub const Geocentric = struct { /// ρ sin φ′: positive in the northern hemisphere. rho_sin_phi: f64, /// ρ cos φ′: never negative. rho_cos_phi: f64,};
/// Where a body appears to an observer.pub const Topocentric = struct { /// Topocentric right ascension and declination, referred to the true /// equator and equinox of date. equatorial: coord.Equatorial, /// Topocentric ecliptic longitude and latitude, referred to the true /// equinox of date. ecliptic: coord.Ecliptic, /// The topocentric hour angle: how far west of the observer's meridian /// the body is, in (−180°, +180°]; negative while it is still rising /// towards the meridian. hour_angle: Angle, /// Azimuth, from the north through the east, and the airless altitude: /// no refraction is applied. `observedAltitude` gives where the body is /// seen. horizontal: coord.Horizontal, /// The distance from the observer to the body's center, in kilometers. distance_km: f64,};
/// The equatorial radius of the WGS 84 ellipsoid, in kilometers.pub const equatorial_radius_km = 6378.137;
/// b/a = 1 − f for the WGS 84 ellipsoid, f = 1/298.257223563.const polar_ratio = 1 - 1.0 / 298.257223563;
/// A body's topocentric position, from its apparent geocentric equatorial/// coordinates and distance at an instant, which may be in either time/// scale. This is what every body's own `topocentric` calls, and what a/// body without one can be passed through.////// The parallax is applied rigorously, as the difference of the body's and/// the observer's geocentric vectors. That is Meeus's eqs. 40.6 and 40.7,/// with sin π the ratio of the Earth's equatorial radius to the distance,/// and it gives the topocentric distance as well. The hour angle comes from/// the apparent sidereal time, the ecliptic coordinates from the/// topocentric equatorial ones with the true obliquity, and the horizontal/// ones from the topocentric hour angle and declination with the/// observer's geodetic latitude (Meeus, eqs. 13.5 and 13.6).////// NaN if the instant, the coordinates or the distance is not finite, or the/// observer is not valid (see `Observer`).pub fn topocentric( instant: Instant, observer: Observer, equatorial: coord.Equatorial, distance_km: f64,) Topocentric { if (!observer.isValid()) return nan_topocentric; const place = observer.geocentric(); // Reduced before anything is subtracted from it: a longitude of 10³⁰⁰° // is finite, and would swallow any right ascension. const local_sidereal_time = earth.greenwichApparentSiderealTime(instant).add(observer.longitude).normalized().radians(); const alpha = equatorial.right_ascension.radians(); const delta = equatorial.declination.radians(); const hour_angle = local_sidereal_time - alpha;
// The body as seen from the Earth's center, in kilometers, in a frame // whose x axis points to the observer's meridian on the equator, y to // the west point and z to the north celestial pole; then the observer's // own position taken away. const x = distance_km * @cos(delta) * @cos(hour_angle) - equatorial_radius_km * place.rho_cos_phi; const y = distance_km * @cos(delta) * @sin(hour_angle); const z = distance_km * @sin(delta) - equatorial_radius_km * place.rho_sin_phi;
const topocentric_hour_angle = math.atan2(y, x); const topocentric_equatorial: coord.Equatorial = .{ .right_ascension = Angle.fromRadians(local_sidereal_time - topocentric_hour_angle).normalized(), .declination = Angle.fromRadians(math.atan2(z, @sqrt(x * x + y * y))), }; const h = Angle.fromRadians(topocentric_hour_angle).signed(); return .{ .equatorial = topocentric_equatorial, .ecliptic = topocentric_equatorial.toEcliptic(earth.trueObliquity(instant)), .hour_angle = h, .horizontal = topocentric_equatorial.toHorizontal(h, observer.latitude), .distance_km = @sqrt(x * x + y * y + z * z), };}
/// The air an observer looks through, which sets how strongly it/// refracts. The defaults are the conditions the refraction formulas were/// written for: sea level, 1010 mb and 10 °C.pub const Atmosphere = struct { /// Pressure at the observer, in millibars (hectopascals). pressure_mb: f64 = 1010, /// Air temperature at the observer, in degrees Celsius. temperature_c: f64 = 10,
/// How much stronger or weaker refraction is than under the default /// conditions: P/1010 × 283/(273 + T) (Meeus, p. 107). Meeus calls it /// only approximately correct, and refraction near the horizon /// wanders by tenths of a degree with the weather whatever it is /// given. fn factor(self: Atmosphere) f64 { return self.pressure_mb / 1010 * 283 / (273 + self.temperature_c); }};
/// Where a body at a given airless altitude is seen: the altitude plus/// the refraction, by Sæmundsson's formula (Meeus, eq. 16.4), with/// Meeus's constant that makes it zero at the zenith, scaled for the/// atmosphere.////// Refraction lifts a body by about 29′ at the horizon, 5.4′ at 10° and/// 1′ at 45°. Below an airless altitude of −1° it is held at its value/// there, 38.80′ in the default atmosphere, as JPL Horizons holds it: the/// formula means nothing so far below the horizon, and holding it keeps/// the observed altitude increasing with the airless one, so that/// `airlessAltitude` has an answer for every observed altitude. Refraction/// a few degrees from the horizon is uncertain by much more than the/// formula's accuracy; see Meeus, p. 107.////// NaN if the altitude is not finite.pub fn observedAltitude(airless: Angle, atmosphere: Atmosphere) Angle { if (!math.isFinite(airless.radians())) return nan; const h = airless.degrees(); return Angle.fromDegrees(h + atmosphere.factor() * refraction(h) / 60);}
/// Where a body seen at a given observed altitude is without the/// atmosphere: the exact inverse of `observedAltitude`, in any atmosphere,/// found by Newton's method.////// Meeus gives Bennett's formula (eq. 16.3) for this direction, but it and/// Sæmundsson's are separate fits that agree only to about 4″ in the/// default atmosphere, and once the refraction is scaled for any other/// they are evaluated at different altitudes and part by minutes of arc/// near the horizon. Inverting Sæmundsson's keeps the two directions/// consistent, and consistent with JPL Horizons.////// NaN if the altitude is not finite.pub fn airlessAltitude(observed: Angle, atmosphere: Atmosphere) Angle { if (!math.isFinite(observed.radians())) return nan; const h0 = observed.degrees(); const f = atmosphere.factor(); // Refraction changes by at most about 0.17° a degree of altitude in // the default atmosphere, so h + f R(h) climbs steadily and Newton's // method, started from h₀ less the refraction there, converges in a // few steps; below −1°, where refraction is held, in one. var h = h0 - f * refraction(h0) / 60; for (0..8) |_| { const step = (h + f * refraction(h) / 60 - h0) / (1 + f * refractionSlope(h) / 60); h -= step; if (@abs(step) < 1e-12) break; } return Angle.fromDegrees(h);}
/// Sæmundsson's refraction at an airless altitude h in degrees, in the/// default atmosphere, in minutes of arc, held below −1°.fn refraction(h: f64) f64 { const a = @max(h, hold_below_degrees); const r = 1.02 / @tan(saemundssonArgument(a) * math.rad_per_deg) + 0.0019279; return @max(r, 0);}
/// The derivative of `refraction`, in minutes of arc a degree.fn refractionSlope(h: f64) f64 { if (h < hold_below_degrees) return 0; const x = saemundssonArgument(h) * math.rad_per_deg; const s = @sin(x); const dx_dh = 1 - 10.3 / ((h + 5.11) * (h + 5.11)); return -1.02 * math.rad_per_deg * dx_dh / (s * s);}
/// h + 10.3/(h + 5.11), in degrees.fn saemundssonArgument(h: f64) f64 { return h + 10.3 / (h + 5.11);}
/// The airless altitude in degrees below which refraction is held.const hold_below_degrees: f64 = -1.0;
const nan = Angle.fromRadians(math.nan(f64));
const nan_topocentric: Topocentric = .{ .equatorial = .{ .right_ascension = nan, .declination = nan }, .ecliptic = .{ .longitude = nan, .latitude = nan }, .hour_angle = nan, .horizontal = .{ .azimuth = nan, .altitude = nan }, .distance_km = math.nan(f64),};