Yesod: self-hosted WebXR rooms

Start a room in a world and hang out in VR from the headset's browser:
A-Frame + Networked-A-Frame over EasyRTC, with a coturn relay that can hide
everyone's IP address. Ships with the Campfire, Temple, Cabin, Balcony and
Test room worlds.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

1 files changed, 248 insertions, 0 deletions

diff --git a/public/js/astro.js b/public/js/astro.js
new file mode 100644
index 0000000..8019e06
--- /dev/null
+++ b/public/js/astro.js
@@ -0,0 +1,248 @@
+// 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);