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>

3 files changed, 155 insertions, 0 deletions

diff --git a/tools/sky/astro_test.js b/tools/sky/astro_test.js
new file mode 100644
index 0000000..b299d16
--- /dev/null
+++ b/tools/sky/astro_test.js
@@ -0,0 +1,76 @@
+// Checks public/js/astro.js against JPL Horizons and Meeus.
+//
+//   deno test tools/sky/astro_test.js
+import '../../public/js/astro.js';
+
+const { Astro } = globalThis;
+const RAD = Math.PI / 180;
+const gap = (a, b) => Math.abs((((a - b + 180) % 360) + 360) % 360 - 180);
+const ok = (cond, msg) => { if (!cond) throw new Error(msg); };
+
+function eclipticOf(dir, ms) {
+  // Inverse of Astro.equatorial, to compare with Horizons' ecliptic longitudes.
+  const ob = (23.4393 - 3.563e-7 * (Astro.julianDay(ms) - 2451543.5)) * RAD;
+  const [x, y, z] = dir;
+  const ye = y * Math.cos(ob) + z * Math.sin(ob);
+  const ze = -y * Math.sin(ob) + z * Math.cos(ob);
+  return { lon: ((Math.atan2(ye, x) / RAD) + 360) % 360, lat: Math.asin(ze) / RAD };
+}
+
+Deno.test('planets are where Horizons puts them (geocentric, of date)', () => {
+  // From android-toolkit's AstroTest: apparent ecliptic longitude and latitude of the date.
+  const refs = {
+    '2026-10-31T20:00:00Z': { sun: [218.4102216, 0.0001252], moon: [115.1119618, 2.6704709], mercury: [226.8441830, -1.5976931], venus: [206.3301706, -4.9067015], mars: [138.5376903, 1.5969459], jupiter: [144.3028948, 0.6895096], saturn: [9.2969189, -2.6851943] },
+    '2028-08-17T06:00:00Z': { sun: [144.8636375, 0.0000135], moon: [101.0244024, 1.0776734], mercury: [165.1873864, 0.4638201], venus: [99.2952316, -2.9103916], mars: [108.2009535, 0.7458269], jupiter: [178.6149976, 1.1269943], saturn: [41.2847987, -2.5175988] },
+  };
+  for (const [iso, bodies] of Object.entries(refs)) {
+    const ms = Date.parse(iso);
+    for (const [name, [lon, lat]] of Object.entries(bodies)) {
+      const e = Astro.ecliptic(name, ms);
+      ok(gap(e.lon, lon) < 0.05, `${name} ${iso} lon ${e.lon} vs ${lon}`);
+      ok(Math.abs(e.lat - lat) < 0.05, `${name} ${iso} lat ${e.lat} vs ${lat}`);
+      const back = eclipticOf(Astro.equatorial(e.lon, e.lat, ms), ms);
+      ok(gap(back.lon, e.lon) < 1e-6 && Math.abs(back.lat - e.lat) < 1e-6, `${name} round trip`);
+    }
+  }
+});
+
+Deno.test("Seattle's sky matches Horizons (airless azimuth, elevation, magnitude)", () => {
+  // Horizons, observer at 47.6062 N, 122.3321 W, 2026-11-01 00:00 UT (JD 2461345.5).
+  const ms = Date.parse('2026-11-01T00:00:00Z');
+  const refs = {
+    sun: [239.356839801, 7.223273431],
+    moon: [338.838357584, -17.481568478, -10.609],
+    mercury: [231.106696194, 8.172959876, 3.268],
+    venus: [249.495509600, -1.327997551, -4.170],
+    mars: [316.936325646, -14.742817319, 0.789],
+    jupiter: [310.975220809, -14.507738850, -2.003],
+    saturn: [87.986528867, -0.198109819, 0.462],
+    uranus: [33.812741373, -14.521612899, 5.611],
+  };
+  const horizon = Astro.horizonMatrix(ms, 47.6062, -122.3321);
+  const sky = Astro.bodies(ms, 47.6062, -122.3321);
+  const found = { sun: { dir: sky.sun }, moon: sky.moon };
+  for (const p of sky.planets) found[p.name] = p;
+  for (const [name, [az, alt, mag]] of Object.entries(refs)) {
+    const h = Astro.altAz(found[name].dir, horizon);
+    ok(Math.abs(h.alt - alt) < 0.05, `${name} alt ${h.alt} vs ${alt}`);
+    ok(gap(h.az, az) < 0.1, `${name} az ${h.az} vs ${az}`);
+    // Schlyter's magnitudes are rough for Mercury near inferior conjunction, as it is here.
+    if (mag !== undefined) ok(Math.abs(found[name].magnitude - mag) < (name === 'mercury' ? 0.5 : 0.25), `${name} mag ${found[name].magnitude} vs ${mag}`);
+  }
+});
+
+Deno.test('precession matches Meeus example 21.b', () => {
+  // Theta Persei, J2000 to 2028 Nov 13.19 TD, with its proper motion (+0.99 s, -2.58") taken back out.
+  const ms = (2462088.69 - 2440587.5) * 86400000;
+  const ra = (2 + 44 / 60 + 11.986 / 3600) * 15 * RAD;
+  const dec = (49 + 13 / 60 + 42.48 / 3600) * RAD;
+  const v = Astro.apply(Astro.precessionMatrix(ms), [Math.cos(dec) * Math.cos(ra), Math.cos(dec) * Math.sin(ra), Math.sin(dec)]);
+  const ra2 = Math.atan2(v[1], v[0]) / RAD / 15;
+  const dec2 = Math.asin(v[2]) / RAD;
+  const wantRa = 2 + 46 / 60 + (11.331 - 0.99) / 3600;
+  const wantDec = 49 + 20 / 60 + (54.54 + 2.58) / 3600;
+  ok(Math.abs(ra2 - wantRa) * 3600 * 15 < 1, `ra ${ra2} vs ${wantRa}`);
+  ok(Math.abs(dec2 - wantDec) * 3600 < 1, `dec ${dec2} vs ${wantDec}`);
+});
diff --git a/tools/sky/build-milkyway.py b/tools/sky/build-milkyway.py
new file mode 100644
index 0000000..c9a63bd
--- /dev/null
+++ b/tools/sky/build-milkyway.py
@@ -0,0 +1,41 @@
+#!/usr/bin/env python3
+# Turns NASA SVS's "Deep Star Maps 2020" Milky Way background (the star map
+# without its Hipparcos and Tycho stars, which the sky draws itself) into the
+# sRGB JPEG the sky component maps onto the celestial sphere.
+#
+#   curl -O https://svs.gsfc.nasa.gov/vis/a000000/a004800/a004851/milkyway_2020_4k.exr
+#   python3 tools/sky/build-milkyway.py milkyway_2020_4k.exr public/assets/sky/milkyway-4k.jpg
+#
+# Needs vips (to read the OpenEXR), numpy and Pillow. The map is plate carrée
+# in J2000 right ascension and declination, centred on 0h with r.a. increasing
+# to the left; this keeps that layout. Values stay linear in brightness up to
+# WHITE (only the heart of the Sagittarius star cloud goes past it), and are
+# stored with the sRGB curve so the faint parts keep their detail in 8 bits.
+import os
+import subprocess
+import sys
+import tempfile
+
+import numpy as np
+from PIL import Image
+
+WHITE = 0.6
+
+
+def main(src, dst):
+    width, height = (int(subprocess.check_output(['vipsheader', '-f', f, src], stderr=subprocess.DEVNULL)) for f in ('width', 'height'))
+    with tempfile.TemporaryDirectory() as tmp:
+        raw = os.path.join(tmp, 'map.raw')
+        subprocess.run(['vips', 'rawsave', src, raw], check=True, stderr=subprocess.DEVNULL)
+        rgb = np.fromfile(raw, dtype=np.float32).reshape(height, width, -1)[..., :3]
+    x = np.clip(rgb / WHITE, 0.0, 1.0)
+    srgb = np.where(x <= 0.0031308, 12.92 * x, 1.055 * np.power(x, 1 / 2.4) - 0.055)
+    # A little noise before rounding keeps the dark gradients from banding.
+    noise = np.random.default_rng(1).uniform(-0.5, 0.5, srgb.shape)
+    out = np.clip(np.round(srgb * 255 + noise), 0, 255).astype(np.uint8)
+    Image.fromarray(out, 'RGB').save(dst, quality=88, optimize=True, progressive=True)
+    print(f'{width}x{height} -> {dst} ({os.path.getsize(dst) // 1024} KB)')
+
+
+if __name__ == '__main__':
+    main(sys.argv[1], sys.argv[2])
diff --git a/tools/sky/build-stars.py b/tools/sky/build-stars.py
new file mode 100644
index 0000000..061cea8
--- /dev/null
+++ b/tools/sky/build-stars.py
@@ -0,0 +1,38 @@
+#!/usr/bin/env python3
+# Packs the Yale Bright Star Catalog (5th ed.) into public/assets/sky/stars.bin
+# for the sky component: four little-endian float32s per star, right ascension
+# and declination (J2000, radians), V magnitude and B-V colour index.
+#
+#   curl -O http://tdc-www.harvard.edu/catalogs/bsc5.dat.gz
+#   python3 tools/sky/build-stars.py bsc5.dat.gz public/assets/sky/stars.bin
+#
+# Byte columns follow the catalog's ReadMe (CDS V/50).
+import gzip
+import math
+import struct
+import sys
+
+
+def main(src, dst):
+    opener = gzip.open if src.endswith('.gz') else open
+    stars = []
+    with opener(src, 'rt', encoding='ascii', errors='replace') as f:
+        for line in f:
+            line = line.rstrip('\n').ljust(197)
+            if not line[75:77].strip() or not line[102:107].strip():
+                continue  # the 14 entries that turned out not to be stars
+            ra = (int(line[75:77]) + int(line[77:79]) / 60 + float(line[79:83]) / 3600) * 15
+            dec = int(line[84:86]) + int(line[86:88]) / 60 + int(line[88:90]) / 3600
+            if line[83] == '-':
+                dec = -dec
+            vmag = float(line[102:107])
+            bv = float(line[109:114]) if line[109:114].strip() else 0.6
+            stars.append((math.radians(ra), math.radians(dec), vmag, bv))
+    with open(dst, 'wb') as out:
+        for s in stars:
+            out.write(struct.pack('<4f', *s))
+    print(f'{len(stars)} stars, {min(s[2] for s in stars):.2f} to {max(s[2] for s in stars):.2f} mag -> {dst}')
+
+
+if __name__ == '__main__':
+    main(sys.argv[1], sys.argv[2])