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.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121// SPDX-FileCopyrightText: © 2026 Jeffrey C. Ollie <jeff@ocjtech.us>// SPDX-License-Identifier: MIT
//! Positions on the celestial sphere, and the conversions between the//! ecliptic, equatorial and horizontal frames (Meeus, *Astronomical//! Algorithms*, 2nd ed., ch. 13).//!//! Neither type says which equinox it is referred to; that is a property of//! where the coordinates came from, and every function returning one says.//! A conversion takes the obliquity belonging to that equinox: the true//! obliquity for apparent coordinates of date, from `astro.earth`.
const std = @import("std");const math = std.math;
const Angle = @import("Angle.zig");
/// A direction measured from the ecliptic and the vernal equinox.pub const Ecliptic = struct { /// Celestial longitude λ, measured eastwards along the ecliptic from the /// equinox, in [0°, 360°). longitude: Angle, /// Celestial latitude β, positive north of the ecliptic, in /// [−90°, +90°]. latitude: Angle,
/// The same direction in equatorial coordinates (Meeus, eqs. 13.3 and /// 13.4). pub fn toEquatorial(self: Ecliptic, obliquity: Angle) Equatorial { const l = self.longitude.radians(); const b = self.latitude.radians(); const e = obliquity.radians(); // atan2 rather than the book's tangent, which loses the quadrant. const ra = math.atan2(@sin(l) * @cos(e) - @tan(b) * @sin(e), @cos(l)); const dec = math.asin(clamp(@sin(b) * @cos(e) + @cos(b) * @sin(e) * @sin(l))); return .{ .right_ascension = Angle.fromRadians(ra).normalized(), .declination = Angle.fromRadians(dec), }; }};
/// A direction measured from the celestial equator and the vernal equinox.pub const Equatorial = struct { /// Right ascension α, measured eastwards along the equator from the /// equinox, in [0°, 360°). `Angle.hours` reads it in hours. right_ascension: Angle, /// Declination δ, positive north of the equator, in [−90°, +90°]. declination: Angle,
/// The same direction in ecliptic coordinates (Meeus, eqs. 13.1 and /// 13.2). pub fn toEcliptic(self: Equatorial, obliquity: Angle) Ecliptic { const a = self.right_ascension.radians(); const d = self.declination.radians(); const e = obliquity.radians(); const lon = math.atan2(@sin(a) * @cos(e) + @tan(d) * @sin(e), @cos(a)); const lat = math.asin(clamp(@sin(d) * @cos(e) - @cos(d) * @sin(e) * @sin(a))); return .{ .longitude = Angle.fromRadians(lon).normalized(), .latitude = Angle.fromRadians(lat), }; }
/// The same direction in an observer's horizontal coordinates, given /// its hour angle there (positive west of the meridian) and the /// observer's geodetic latitude (Meeus, eqs. 13.5 and 13.6, with the /// azimuth turned to count from the north; see `Horizontal`). Give it /// topocentric coordinates and the topocentric hour angle, or the /// result is out by the body's parallax. No refraction is applied: the /// altitude is the airless one. pub fn toHorizontal(self: Equatorial, hour_angle: Angle, latitude: Angle) Horizontal { const h = hour_angle.radians(); const d = self.declination.radians(); const phi = latitude.radians(); // The direction's components towards the north point, the east // point and the zenith. Meeus's tan A = sin H / (cos H sin φ − // tan δ cos φ), from the south, is the first two multiplied // through by cos δ, so that a body at a pole needs no tangent, and // with both signs reversed to count from the north. His // sin h (eq. 13.6) is the third; the altitude is taken by atan2 // rather than asin, which loses a millionth of a degree near the // zenith. const north = @sin(d) * @cos(phi) - @cos(d) * @cos(h) * @sin(phi); const east = -@cos(d) * @sin(h); const up = @sin(phi) * @sin(d) + @cos(phi) * @cos(d) * @cos(h); const azimuth = math.atan2(east, north); const altitude = math.atan2(up, @sqrt(north * north + east * east)); return .{ .azimuth = Angle.fromRadians(azimuth).normalized(), .altitude = Angle.fromRadians(altitude), }; }};
/// A direction as an observer on the Earth sees it, measured from the/// horizon and the north point.////// Azimuth counts from the north through the east: 0° north, 90° east,/// 180° south, 270° west, as compasses, almanacs and JPL Horizons count it./// Meeus counts from the *south*, westwards (p. 91); add 180° to his.pub const Horizontal = struct { /// Azimuth, from the north through the east, in [0°, 360°). At the /// zenith and the nadir it has no meaning, and comes out as whatever /// the rounding gives. azimuth: Angle, /// Altitude above the horizon, negative below it, in [−90°, +90°]. It /// is the airless altitude unless it says otherwise: the atmosphere /// raises what is seen by up to about 0.6°, at the horizon, and /// `astro.observer.observedAltitude` says by how much. altitude: Angle,};
/// Rounding can carry a sine a hair past ±1 at the poles, where `asin`/// would give NaN for a perfectly good direction. NaN stays NaN.fn clamp(x: f64) f64 { if (x > 1) return 1; if (x < -1) return -1; return x;}