path: root/public/js/astro.js

// Where the Sun, Moon and planets are, for the live sky. Ported from the
// android-toolkit :astro module (Astro.kt and Planets.kt): Paul Schlyter's "How
// to compute planetary positions" for the Moon and planets, good to a minute or
// two of arc, and the NOAA equations for the Sun's longitude. Uranus, which the
// Kotlin module leaves out, uses Schlyter's elements and perturbations too.
//
// Angles are in degrees unless a name says otherwise, and times are Unix
// milliseconds. Directions are unit vectors [x, y, z] in equatorial
// coordinates: x toward the equinox, z toward the north celestial pole. A
// matrix is three rows. Like the rest of public/js this is a plain script; it
// defines the global `Astro`, which tools/sky/astro_test.js checks.
(function (root) {
  const RAD = Math.PI / 180;
  const sin = (d) => Math.sin(d * RAD);
  const cos = (d) => Math.cos(d * RAD);
  const atan2d = (y, x) => Math.atan2(y, x) / RAD;
  const asind = (x) => Math.asin(Math.min(1, Math.max(-1, x))) / RAD;
  const mod = (a, n) => ((a % n) + n) % n;

  /** Earth radii in an astronomical unit. */
  const AU = 23454.8;

  const julianDay = (ms) => ms / 86400000 + 2440587.5;
  /** Days from Schlyter's epoch, 2000 January 0.0. */
  const dayNumber = (ms) => julianDay(ms) - 2451543.5;

  /** The Sun's apparent ecliptic longitude (NOAA solar calculator). */
  function sunLongitude(ms) {
    const jc = (julianDay(ms) - 2451545) / 36525;
    const l0 = mod(280.46646 + jc * (36000.76983 + jc * 0.0003032), 360);
    const m = 357.52911 + jc * (35999.05029 - 0.0001537 * jc);
    const c = sin(m) * (1.914602 - jc * (0.004817 + 0.000014 * jc)) +
      sin(2 * m) * (0.019993 - 0.000101 * jc) +
      sin(3 * m) * 0.000289;
    const omega = 125.04 - 1934.136 * jc;
    return mod(l0 + c - 0.00569 - 0.00478 * sin(omega), 360);
  }

  /**
   * Schlyter's elements for day d: node, inclination, argument of perihelion,
   * semi-major axis (AU; Earth radii for the Moon), eccentricity, mean anomaly.
   * The Sun's are really the Earth's orbit seen from the Earth.
   */
  const ELEMENTS = {
    sun: (d) => [0, 0, 282.9404 + 4.70935e-5 * d, 1, 0.016709 - 1.151e-9 * d, 356.0470 + 0.9856002585 * d],
    moon: (d) => [125.1228 - 0.0529538083 * d, 5.1454, 318.0634 + 0.1643573223 * d, 60.2666, 0.0549, 115.3654 + 13.0649929509 * d],
    mercury: (d) => [48.3313 + 3.24587e-5 * d, 7.0047 + 5.00e-8 * d, 29.1241 + 1.01444e-5 * d, 0.387098, 0.205635 + 5.59e-10 * d, 168.6562 + 4.0923344368 * d],
    venus: (d) => [76.6799 + 2.46590e-5 * d, 3.3946 + 2.75e-8 * d, 54.8910 + 1.38374e-5 * d, 0.723330, 0.006773 - 1.302e-9 * d, 48.0052 + 1.6021302244 * d],
    mars: (d) => [49.5574 + 2.11081e-5 * d, 1.8497 - 1.78e-8 * d, 286.5016 + 2.92961e-5 * d, 1.523688, 0.093405 + 2.516e-9 * d, 18.6021 + 0.5240207766 * d],
    jupiter: (d) => [100.4542 + 2.76854e-5 * d, 1.3030 - 1.557e-7 * d, 273.8777 + 1.64505e-5 * d, 5.20256, 0.048498 + 4.469e-9 * d, 19.8950 + 0.0830853001 * d],
    saturn: (d) => [113.6634 + 2.38980e-5 * d, 2.4886 - 1.081e-7 * d, 339.3939 + 2.97661e-5 * d, 9.55475, 0.055546 - 9.499e-9 * d, 316.9670 + 0.0334442282 * d],
    uranus: (d) => [74.0005 + 1.3978e-5 * d, 0.7733 + 1.9e-8 * d, 96.6612 + 3.0565e-5 * d, 19.18171 - 1.55e-8 * d, 0.047318 + 7.45e-9 * d, 142.5905 + 0.011725806 * d],
  };

  const anomaly = (body, d) => ELEMENTS[body](d)[5];

  /** Where a body is in its orbit, as ecliptic { lon, lat, r }: heliocentric for a planet, geocentric for the Sun and Moon. */
  function orbitPosition(body, d) {
    const [node, incl, peri, axis, e, mean] = ELEMENTS[body](d);
    const m = mod(mean, 360);
    let ea = m + (e * sin(m) * (1 + e * cos(m))) / RAD;
    for (let i = 0; i < 4; i++) ea -= (ea - (e * sin(ea)) / RAD - m) / (1 - e * cos(ea));
    const xv = axis * (cos(ea) - e);
    const yv = axis * Math.sqrt(1 - e * e) * sin(ea);
    const r = Math.hypot(xv, yv);
    const u = atan2d(yv, xv) + peri;
    const x = r * (cos(node) * cos(u) - sin(node) * sin(u) * cos(incl));
    const y = r * (sin(node) * cos(u) + cos(node) * sin(u) * cos(incl));
    const z = r * sin(u) * sin(incl);
    return { lon: mod(atan2d(y, x), 360), lat: atan2d(z, Math.hypot(x, y)), r };
  }

  function moonEcliptic(d) {
    const p = orbitPosition('moon', d);
    const [sNode, , sPeri, , , ms] = ELEMENTS.sun(d);
    const [mNode, , mPeri, , , mm] = ELEMENTS.moon(d);
    const elong = (mm + mPeri + mNode) - (ms + sPeri + sNode);
    const f = mm + mPeri;
    const lon = p.lon -
      1.274 * sin(mm - 2 * elong) + 0.658 * sin(2 * elong) - 0.186 * sin(ms) -
      0.059 * sin(2 * mm - 2 * elong) - 0.057 * sin(mm - 2 * elong + ms) + 0.053 * sin(mm + 2 * elong) +
      0.046 * sin(2 * elong - ms) + 0.041 * sin(mm - ms) - 0.035 * sin(elong) - 0.031 * sin(mm + ms) -
      0.015 * sin(2 * f - 2 * elong) + 0.011 * sin(mm - 4 * elong);
    const lat = p.lat -
      0.173 * sin(f - 2 * elong) - 0.055 * sin(mm - f - 2 * elong) - 0.046 * sin(mm + f - 2 * elong) +
      0.033 * sin(f + 2 * elong) + 0.017 * sin(2 * mm + f);
    const r = p.r - 0.58 * cos(mm - 2 * elong) - 0.46 * cos(2 * elong);
    return { lon: mod(lon, 360), lat, r };
  }

  /** Jupiter and Saturn pull on each other enough to move them by most of a degree, and both nudge Uranus. */
  function perturb(body, p, d) {
    const mj = anomaly('jupiter', d);
    const ms = anomaly('saturn', d);
    if (body === 'jupiter') {
      p.lon += -0.332 * sin(2 * mj - 5 * ms - 67.6) - 0.056 * sin(2 * mj - 2 * ms + 21) +
        0.042 * sin(3 * mj - 5 * ms + 21) - 0.036 * sin(mj - 2 * ms) + 0.022 * cos(mj - ms) +
        0.023 * sin(2 * mj - 3 * ms + 52) - 0.016 * sin(mj - 5 * ms - 69);
    } else if (body === 'saturn') {
      p.lon += 0.812 * sin(2 * mj - 5 * ms - 67.6) - 0.229 * cos(2 * mj - 4 * ms - 2) +
        0.119 * sin(mj - 2 * ms - 3) + 0.046 * sin(2 * mj - 6 * ms - 69) + 0.014 * sin(mj - 3 * ms + 32);
      p.lat += -0.020 * cos(2 * mj - 4 * ms - 2) + 0.018 * sin(2 * mj - 6 * ms - 49);
    } else if (body === 'uranus') {
      const mu = anomaly('uranus', d);
      p.lon += 0.040 * sin(ms - 2 * mu + 6) + 0.035 * sin(ms - 3 * mu + 33) - 0.015 * sin(mj - mu + 20);
    }
    return p;
  }

  /**
   * Where a body is in the zodiac of the date, seen from the Earth's centre:
   * { lon, lat, dist } with dist in AU (Earth radii for the Moon), plus the
   * planets' distance from the Sun, r.
   */
  function ecliptic(body, ms) {
    const d = dayNumber(ms);
    const sun = orbitPosition('sun', d);
    if (body === 'sun') return { lon: sunLongitude(ms), lat: 0, dist: sun.r };
    if (body === 'moon') {
      const m = moonEcliptic(d);
      return { lon: m.lon, lat: m.lat, dist: m.r };
    }
    const p = perturb(body, orbitPosition(body, d), d);
    // Move the origin from the Sun to the Earth.
    const x = p.r * cos(p.lon) * cos(p.lat) + sun.r * cos(sun.lon);
    const y = p.r * sin(p.lon) * cos(p.lat) + sun.r * sin(sun.lon);
    const z = p.r * sin(p.lat);
    return { lon: mod(atan2d(y, x), 360), lat: atan2d(z, Math.hypot(x, y)), dist: Math.hypot(x, y, z), r: p.r, sunDist: sun.r };
  }

  /** Ecliptic longitude and latitude of the date to an equatorial direction of the date. */
  function equatorial(lon, lat, ms) {
    const ob = 23.4393 - 3.563e-7 * dayNumber(ms);
    return [
      cos(lon) * cos(lat),
      sin(lon) * cos(lat) * cos(ob) - sin(lat) * sin(ob),
      sin(lon) * cos(lat) * sin(ob) + sin(lat) * cos(ob),
    ];
  }

  /** Greenwich mean sidereal time, in degrees. */
  function siderealTime(ms) {
    const d = julianDay(ms) - 2451545;
    const c = d / 36525;
    return mod(280.46061837 + 360.98564736629 * d + c * c * (0.000387933 - c / 38710000), 360);
  }

  const multiply = (a, b) => a.map((row) => [0, 1, 2].map((j) => row[0] * b[0][j] + row[1] * b[1][j] + row[2] * b[2][j]));
  const apply = (m, v) => m.map((row) => row[0] * v[0] + row[1] * v[1] + row[2] * v[2]);
  const rotZ = (a) => [[cos(a), -sin(a), 0], [sin(a), cos(a), 0], [0, 0, 1]];

  /**
   * Turns J2000 directions (the star catalogue and the Milky Way map) into
   * directions of the date: the IAU 1976 precession angles, as in Meeus,
   * Astronomical Algorithms, chapter 21.
   */
  function precessionMatrix(ms) {
    const t = (julianDay(ms) - 2451545) / 36525;
    const zeta = (2306.2181 * t + 0.30188 * t * t + 0.017998 * t ** 3) / 3600;
    const z = (2306.2181 * t + 1.09468 * t * t + 0.018203 * t ** 3) / 3600;
    const theta = (2004.3109 * t - 0.42665 * t * t - 0.041833 * t ** 3) / 3600;
    const tilt = [[cos(theta), 0, -sin(theta)], [0, 1, 0], [sin(theta), 0, cos(theta)]];
    return multiply(rotZ(z), multiply(tilt, rotZ(zeta)));
  }

  /**
   * Turns equatorial directions of the date into the local sky at a latitude
   * and (east-positive) longitude: rows give the components toward North,
   * East and the zenith.
   */
  function horizonMatrix(ms, latitude, longitude) {
    const lst = siderealTime(ms) + longitude;
    const sl = sin(lst);
    const cl = cos(lst);
    const sp = sin(latitude);
    const cp = cos(latitude);
    return [
      [-sp * cl, -sp * sl, cp],
      [-sl, cl, 0],
      [cp * cl, cp * sl, sp],
    ];
  }

  /** Altitude and azimuth (clockwise from North) of an equatorial direction, given a horizonMatrix. */
  function altAz(dir, horizon) {
    const [n, e, u] = apply(horizon, dir);
    return { alt: asind(u), az: mod(atan2d(e, n), 360) };
  }

  const PLANETS = ['mercury', 'venus', 'mars', 'jupiter', 'saturn', 'uranus'];

  /** Schlyter's visual magnitudes, from the distances to the Sun (r) and Earth (dist) in AU and the phase angle. */
  const MAGNITUDE = {
    mercury: (fv) => -0.36 + 0.027 * fv + 2.2e-13 * fv ** 6,
    venus: (fv) => -4.34 + 0.013 * fv + 4.2e-7 * fv ** 3,
    mars: (fv) => -1.51 + 0.016 * fv,
    jupiter: (fv) => -9.25 + 0.014 * fv,
    saturn: (fv) => -9.0 + 0.044 * fv,
    uranus: (fv) => -7.15 + 0.001 * fv,
  };

  /**
   * Everything that moves against the stars, at a time and place: equatorial
   * directions of the date and visual magnitudes. The Moon's direction is
   * topocentric (from the observer, not the Earth's centre: up to a degree
   * lower), and it comes with its distance and phase.
   */
  function bodies(ms, latitude, longitude) {
    const d = dayNumber(ms);
    const sunE = ecliptic('sun', ms);
    const sun = equatorial(sunE.lon, sunE.lat, ms);

    const moonE = ecliptic('moon', ms);
    const geo = equatorial(moonE.lon, moonE.lat, ms);
    const lst = siderealTime(ms) + longitude;
    const observer = [cos(latitude) * cos(lst), cos(latitude) * sin(lst), sin(latitude)];
    const topo = geo.map((g, i) => g * moonE.dist - observer[i]);
    const topoDist = Math.hypot(...topo);
    // The Sun is 400 times farther than the Moon, so the angle at the Moon is
    // the supplement of the Moon's elongation from the Sun, to a tenth of a degree.
    const elongation = Math.acos(Math.max(-1, Math.min(1, apply([sun], geo)[0]))) / RAD;
    const phaseAngle = 180 - elongation;
    const moon = {
      dir: topo.map((v) => v / topoDist),
      dist: topoDist,
      phaseAngle,
      illuminated: (1 + cos(phaseAngle)) / 2,
      magnitude: 0.23 + 5 * Math.log10(sunE.dist * (moonE.dist / AU)) + 0.026 * phaseAngle + 4.0e-9 * phaseAngle ** 4,
    };

    const planets = PLANETS.map((name) => {
      const e = ecliptic(name, ms);
      const cosFv = (e.r * e.r + e.dist * e.dist - e.sunDist * e.sunDist) / (2 * e.r * e.dist);
      const fv = Math.acos(Math.max(-1, Math.min(1, cosFv))) / RAD;
      let magnitude = MAGNITUDE[name](fv) + 5 * Math.log10(e.r * e.dist);
      if (name === 'saturn') {
        // The rings add up to a magnitude and a half when they're open; in 2025 they were edge-on.
        const b = asind(sin(e.lat) * cos(28.06) - cos(e.lat) * sin(28.06) * sin(e.lon - 169.51 - 3.82e-5 * d));
        magnitude += -2.6 * Math.abs(sin(b)) + 1.2 * sin(b) ** 2;
      }
      return { name, dir: equatorial(e.lon, e.lat, ms), magnitude };
    });

    return { sun, moon, planets };
  }

  root.Astro = { julianDay, siderealTime, ecliptic, equatorial, precessionMatrix, horizonMatrix, altAz, bodies, multiply, apply };
})(globalThis);