fix(sky): keep the Sun and Moon, settle outside the ICRF tables

- Pass the custom sky box to the Viewer instead of `skyBox: false`, which
  also stopped CesiumWidget from creating the Sun and the Moon. Guarded by
  scripts/qa-sky-sun-moon.mjs.
- Outside the IAU 2006 tables (before 1975, after ~2050) the preload
  resolves without data; remember that span and stay on TEME instead of
  preloading and rendering every frame.
- The anime style no longer cel-shades the far plane, where four bands
  turned the faint Milky Way into blotches.
- Paint the planets out after the star-removing opening and with the
  ring median, so no brighter plateau is left; faces rebuilt.

Co-Authored-By: Claude Opus 5.5 <[email protected]>

scripts/build-milky-way.mjs

@@ -15,12 +15,12 @@
  *
  * Transform (deterministic):
  *   1. Resample to WORK_WIDTH×WORK_WIDTH/2 and convert sRGB to linear light.
- *   2. Paint out the planets the 2008–2010 exposures caught (PLANETS): they
- *      sit where today's sky has nothing and are far brighter than any star.
- *   3. Remove stars with a morphological opening (minimum, then maximum, over
+ *   2. Remove stars with a morphological opening (minimum, then maximum, over
  *      a (2·OPEN_RADIUS+1)² window): anything smaller than the window — every
  *      star — drops to the surrounding sky; the renderer draws catalogue stars
  *      itself, sharp and in the right place.
+ *   3. Paint out the planets the 2008–2010 exposures caught (PLANETS): they
+ *      sit where today's sky has nothing and are far brighter than any star.
  *   4. Blur with a three-pass box approximation of a Gaussian, σ = BLUR_DEG on
  *      the sky: rows wrap in longitude and widen by 1/cos(latitude) so the
  *      blur stays round towards the poles.
@@ -234,16 +234,18 @@ function galacticVector(l, b) {
 }
 
 /**
- * Replace a disc of the sky with the mean of the ring just outside it, per
- * channel. The blur that follows hides the seam.
+ * Replace a disc of the sky with the median of the ring just outside it, per
+ * channel. A median, not a mean: the planet's diffraction spikes reach far
+ * past the disc and would lift a mean well above the sky. The blur that
+ * follows hides the seam.
  */
 function paintOut(planes, width, height, { l, b, radiusDeg }) {
   const centre = galacticVector(l, b);
   const inner = Math.cos((radiusDeg * Math.PI) / 180);
-  const outer = Math.cos(((radiusDeg + 0.5) * Math.PI) / 180);
+  const outer = Math.cos(((radiusDeg + 1) * Math.PI) / 180);
   const disc = [];
   const ring = [];
-  const span = Math.ceil(((radiusDeg + 0.5) / 180) * height);
+  const span = Math.ceil(((radiusDeg + 1) / 180) * height);
   const row = Math.round((0.5 - b / 180) * height);
   for (let y = Math.max(0, row - span); y < Math.min(height, row + span); y++) {
     const latitude = 90 - ((y + 0.5) * 180) / height;
@@ -256,10 +258,9 @@ function paintOut(planes, width, height, { l, b, radiusDeg }) {
     }
   }
   for (const plane of planes) {
-    let sum = 0;
-    for (const index of ring) sum += plane[index];
-    const mean = sum / ring.length;
-    for (const index of disc) plane[index] = mean;
+    const values = Float32Array.from(ring, (index) => plane[index]).sort();
+    const median = values[values.length >> 1];
+    for (const index of disc) plane[index] = median;
   }
 }
 
@@ -323,12 +324,14 @@ async function main() {
   const { default: sharp } = await import('sharp');
   const { width, height, planes } = await decodePanorama(sharp, file);
   const sigmaPx = (BLUR_DEG / 180) * height;
-  for (const planet of PLANETS) paintOut(planes, width, height, planet);
   for (const plane of planes) {
     extremumPlane(plane, width, height, OPEN_RADIUS, Math.min);
     extremumPlane(plane, width, height, OPEN_RADIUS, Math.max);
-    blurPlane(plane, width, height, sigmaPx);
   }
+  // After the opening: it sinks the starry sky to the dark gaps between
+  // stars, and a fill taken before it would stay brighter than the sky.
+  for (const planet of PLANETS) paintOut(planes, width, height, planet);
+  for (const plane of planes) blurPlane(plane, width, height, sigmaPx);
   const black = luminancePercentile(planes, BLACK_PERCENTILE);
   const white = luminancePercentile(planes, WHITE_PERCENTILE);
   let total = 0;

scripts/format-scope.json

@@ -1067,5 +1067,6 @@
   "scripts/build-milky-way.mjs",
   "src/sky/stars.test.mjs",
   "src/sky/starrySky.test.mjs",
-  "src/tooling/skyAssets.test.mjs"
+  "src/tooling/skyAssets.test.mjs",
+  "src/tooling/skyAssets.paintOut.test.mjs"
 ]

scripts/qa-sky-sun-moon.mjs

@@ -0,0 +1,96 @@
+#!/usr/bin/env node
+/**
+ * qa-sky-sun-moon — the custom sky must keep Cesium's Sun and Moon.
+ *
+ * CesiumWidget creates `scene.sun` and `scene.moon` only when the `skyBox`
+ * option is not `false`, so replacing the sky box the wrong way silently
+ * removes the Sun disc and its bloom. Points the camera at the Sun from
+ * 22 000 km and fails unless both bodies exist; saves a screenshot.
+ *
+ * Usage: node scripts/qa-sky-sun-moon.mjs [--url http://localhost:4317] [--out qa-shots/sky-sun.png]
+ */
+import puppeteer from 'puppeteer';
+import { mkdirSync } from 'node:fs';
+import path from 'node:path';
+
+const argv = process.argv;
+const option = (name, fallback) =>
+  argv.includes(name) ? argv[argv.indexOf(name) + 1] : fallback;
+const url = option('--url', 'http://localhost:4317');
+const out = option('--out', 'qa-shots/sky-sun.png');
+mkdirSync(path.dirname(out), { recursive: true });
+
+const browser = await puppeteer.launch({
+  headless: 'new',
+  args: ['--no-sandbox', '--ignore-gpu-blocklist'],
+});
+try {
+  const page = await browser.newPage();
+  await page.setViewport({ width: 1280, height: 800 });
+  await page.evaluateOnNewDocument(() =>
+    localStorage.setItem('gev:first-run-mission:v1', 'suppressed'),
+  );
+  await page.goto(url, { waitUntil: 'domcontentloaded' });
+  await page.waitForFunction(() => window.__godsEyeView?.viewer?.scene, {
+    timeout: 90_000,
+  });
+  await new Promise((r) => setTimeout(r, 10_000));
+  const info = await page.evaluate(async () => {
+    // The page's own Cesium instance, as Vite served it.
+    const entry = performance
+      .getEntriesByType('resource')
+      .map((e) => e.name)
+      .find((n) => /\.vite\/deps\/cesium\.js/.test(n));
+    const C = await import(entry);
+    const { viewer } = window.__godsEyeView;
+    const scene = viewer.scene;
+    scene.camera.cancelFlight();
+    const time = viewer.clock.currentTime;
+    const sun = C.Cartesian3.normalize(
+      C.Matrix3.multiplyByVector(
+        C.Transforms.computeTemeToPseudoFixedMatrix(time),
+        C.Simon1994PlanetaryPositions.computeSunPositionInEarthInertialFrame(
+          time,
+        ),
+        new C.Cartesian3(),
+      ),
+      new C.Cartesian3(),
+    );
+    // Off to the side of the Earth, looking straight at the Sun.
+    const side = C.Cartesian3.normalize(
+      C.Cartesian3.cross(sun, C.Cartesian3.UNIT_Z, new C.Cartesian3()),
+      new C.Cartesian3(),
+    );
+    const destination = C.Cartesian3.add(
+      C.Cartesian3.multiplyByScalar(side, 2.2e7, new C.Cartesian3()),
+      C.Cartesian3.multiplyByScalar(sun, -0.6e7, new C.Cartesian3()),
+      new C.Cartesian3(),
+    );
+    const up = C.Cartesian3.normalize(
+      C.Cartesian3.cross(
+        C.Cartesian3.cross(sun, C.Cartesian3.UNIT_Z, new C.Cartesian3()),
+        sun,
+        new C.Cartesian3(),
+      ),
+      new C.Cartesian3(),
+    );
+    scene.camera.setView({ destination, orientation: { direction: sun, up } });
+    scene.requestRender();
+    return {
+      sun: Boolean(scene.sun),
+      moon: Boolean(scene.moon),
+      skyBox: scene.skyBox?.constructor?.name,
+    };
+  });
+  await new Promise((r) => setTimeout(r, 4_000));
+  await page.evaluate(() => window.__godsEyeView.viewer.scene.requestRender());
+  await new Promise((r) => setTimeout(r, 1_500));
+  await page.screenshot({ path: out });
+  console.log(JSON.stringify({ url, ...info, screenshot: out }));
+  if (!info.sun || !info.moon) {
+    console.error('FAIL: scene.sun and scene.moon must both exist');
+    process.exitCode = 1;
+  }
+} finally {
+  await browser.close();
+}

src/app/viewer.js

@@ -120,9 +120,10 @@ export function createApplicationViewer({ container, creditContainer }) {
     selectionIndicator: false,
     infoBox: false,
     baseLayer: false,
-    // Replaced below: the default Tycho-2 cube map is too coarse for the
-    // whole-globe view and its JPEG blocks show around the Earth.
-    skyBox: false,
+    // Instead of the default Tycho-2 cube map, whose JPEG blocks show around
+    // the whole globe. Passed here, not assigned later: `skyBox: false` would
+    // also stop CesiumWidget from creating the Sun and the Moon.
+    skyBox: new StarrySky(),
     creditContainer,
     msaaSamples: 4,
     contextOptions: { webgl: { preserveDrawingBuffer: true } },
@@ -134,7 +135,6 @@ export function createApplicationViewer({ container, creditContainer }) {
     // render loop. See app/atmosphereCompat.js.
     applyModelAtmosphereWorkaround(viewer.scene);
     viewer.scene.globe.show = false;
-    viewer.scene.skyBox = new StarrySky();
     viewer.scene.skyAtmosphere.show = true;
     viewer.scene.skyAtmosphere.atmosphereLightIntensity = 18;
     viewer.scene.skyAtmosphere.saturationShift = -0.12;

src/data/local_data/sky/README.md

@@ -42,10 +42,11 @@ Centre.
 (<https://www.eso.org/public/outreach/copyright/>). **Credit: ESO/S. Brunier.**
 
 **Changes made** (CC BY 4.0 asks for them to be stated): resampled to
-4096×2048, two planets (Jupiter, by position) painted over with the surrounding sky,
-stars removed with a 7×7 morphological opening, blurred (σ = 0.4°),
-reprojected from galactic coordinates to ICRF cube faces, black level set at
-the darkest 25 % of the sky, gamma-encoded and saved as WebP. Build:
+4096×2048, stars removed with a 7×7 morphological opening, two planets
+(Jupiter, by position) painted over with the median of the surrounding sky,
+blurred (σ = 0.4°), reprojected from galactic coordinates to ICRF cube faces,
+black level set at the darkest 25 % of the sky, gamma-encoded and saved as
+WebP. Build:
 `node scripts/build-milky-way.mjs <eso0932a.tif>`.
 
 The panorama is a photographic mosaic, not an astrometric map: bright stars

src/sky/starrySky.js

@@ -174,6 +174,7 @@ export class StarrySky {
     this._rotationPrecise = false;
     this._icrfRequest = undefined;
     this._icrfUnavailable = false;
+    this._icrfMissing = undefined;
     this._loading = false;
     this._starData = undefined;
     this._faceImages = undefined;
@@ -237,7 +238,7 @@ export class StarrySky {
   _updateRotation(time) {
     // An approximate rotation is final only while nothing better can come.
     const settled =
-      this._rotationPrecise || this._icrfRequest || this._icrfUnavailable;
+      this._rotationPrecise || this._icrfRequest || !this._canPreloadIcrf(time);
     if (settled && isSkyRotationCurrent(this._rotationTime, time)) return;
     const { precise } = computeSkyRotation(time, this._rotation);
     this._rotationPrecise = precise;
@@ -245,9 +246,20 @@ export class StarrySky {
     if (!precise) this._preloadIcrf(time);
   }
 
-  /** Fetch Cesium's IAU 2006 tables around `time` once; retry on new spans. */
+  /** Whether loading Cesium's IAU 2006 tables could still help at `time`. */
+  _canPreloadIcrf(time) {
+    return (
+      !this._icrfUnavailable &&
+      !(
+        this._icrfMissing &&
+        Cesium.TimeInterval.contains(this._icrfMissing, time)
+      )
+    );
+  }
+
+  /** Fetch Cesium's IAU 2006 tables around `time` once per span. */
   _preloadIcrf(time) {
-    if (this._icrfRequest || this._icrfUnavailable) return;
+    if (this._icrfRequest || !this._canPreloadIcrf(time)) return;
     const interval = new Cesium.TimeInterval({
       start: Cesium.JulianDate.addDays(time, -1, new Cesium.JulianDate()),
       stop: Cesium.JulianDate.addDays(time, 1, new Cesium.JulianDate()),
@@ -257,7 +269,14 @@ export class StarrySky {
     )
       .then(
         () => {
+          // The tables cover roughly 1975–2050; outside them the preload
+          // resolves with nothing, and TEME stays for this span.
+          if (!Cesium.Transforms.computeIcrfToFixedMatrix(time)) {
+            this._icrfMissing = interval;
+            return;
+          }
           this._rotationTime = undefined;
+          if (!this._destroyed) this._requestRender();
         },
         (error) => {
           // TEME is within a fraction of a degree; do not retry every frame.
@@ -267,7 +286,6 @@ export class StarrySky {
       )
       .finally(() => {
         this._icrfRequest = undefined;
-        if (!this._destroyed) this._requestRender();
       });
   }
 

src/sky/starrySky.test.mjs

@@ -87,3 +87,31 @@ test('the sky draws nothing when hidden, outside 3D, or in a pick pass', () => {
   assert.equal(sky.destroy(), undefined);
   assert.equal(sky.isDestroyed(), true);
 });
+
+const frame = () => new Promise((resolve) => setImmediate(resolve));
+
+for (const iso of ['1970-01-01T00:00:00Z', '2060-01-01T00:00:00Z']) {
+  test(`a clock outside the IAU 2006 tables (${iso}) settles on TEME instead of re-rendering every frame`, async () => {
+    const time = Cesium.JulianDate.fromIso8601(iso);
+    let renders = 0;
+    let preloads = 0;
+    const original = Cesium.Transforms.preloadIcrfFixed;
+    Cesium.Transforms.preloadIcrfFixed = (interval) => {
+      preloads += 1;
+      return original(interval);
+    };
+    try {
+      const sky = new StarrySky({ requestRender: () => (renders += 1) });
+      for (let i = 0; i < 20; i += 1) {
+        sky._updateRotation(time); // what update() does once per frame
+        await frame();
+      }
+      assert.ok(
+        renders <= 1 && preloads <= 1,
+        `20 frames at a frozen clock -> ${preloads} preloads, ${renders} render requests (expected <= 1 each)`,
+      );
+    } finally {
+      Cesium.Transforms.preloadIcrfFixed = original;
+    }
+  });
+}

src/styles/anime.js

@@ -12,6 +12,7 @@ export const animeShader = {
   },
   fragmentShader: /* glsl */ `
     uniform sampler2D colorTexture;
+    uniform sampler2D depthTexture;
     uniform vec2 colorTextureDimensions;
     uniform float intensity;
     uniform float saturation;
@@ -38,9 +39,12 @@ export const animeShader = {
       vec2 texel = 1.0 / colorTextureDimensions;
       vec4 color = texture(colorTexture, uv);
 
-      // Cel-shade: quantize luminance into bands
+      // Cel-shade: quantize luminance into bands. Not the far plane: four
+      // bands turn the faint Milky Way behind the globe into blotches. Its
+      // packed depth reads a few 1e-7 short of 1.0.
       float luma = dot(color.rgb, vec3(0.299, 0.587, 0.114));
-      float bands = mix(256.0, 4.0, intensity);
+      bool space = czm_readDepth(depthTexture, uv) > 0.99999;
+      float bands = space ? 256.0 : mix(256.0, 4.0, intensity);
       float celLuma = floor(luma * bands + 0.5) / bands;
       float lumaScale = (luma > 0.001) ? celLuma / luma : 1.0;
       vec3 celColor = color.rgb * lumaScale;

src/tooling/skyAssets.paintOut.test.mjs

@@ -0,0 +1,64 @@
+import assert from 'node:assert/strict';
+import test from 'node:test';
+import sharp from 'sharp';
+import {
+  CUBE_FACES,
+  FACE_FILE_SUFFIX,
+  PLANETS,
+  cubeFaceDirection,
+  icrfToGalactic,
+} from '../../scripts/build-milky-way.mjs';
+import { raDecToUnitVector } from '../sky/stars.js';
+
+const SKY = new URL('../data/local_data/sky/', import.meta.url);
+// Allowed excess of the painted disc over the surrounding sky, in 8-bit
+// levels of the encoded face: the blur and dither move values by a few.
+const TOLERANCE = 4;
+
+test('painted-out planets leave no disc brighter than the sky around them', async () => {
+  const centres = PLANETS.map((planet) => ({
+    planet,
+    vector: raDecToUnitVector(planet.l, planet.b),
+    disc: [],
+    ring: [],
+  }));
+  for (const face of CUBE_FACES) {
+    const { data, info } = await sharp(
+      new URL(`milkyway_${FACE_FILE_SUFFIX[face]}.webp`, SKY).pathname,
+    )
+      .raw()
+      .toBuffer({ resolveWithObject: true });
+    const size = info.width;
+    for (let j = 0; j < size; j++) {
+      for (let i = 0; i < size; i++) {
+        const d = cubeFaceDirection(face, (i + 0.5) / size, (j + 0.5) / size);
+        const n = Math.hypot(...d);
+        const g = icrfToGalactic(d.map((x) => x / n));
+        for (const centre of centres) {
+          const cos =
+            g[0] * centre.vector[0] +
+            g[1] * centre.vector[1] +
+            g[2] * centre.vector[2];
+          const degrees = (Math.acos(Math.min(1, cos)) * 180) / Math.PI;
+          const k = (j * size + i) * info.channels;
+          const luma = (data[k] + data[k + 1] + data[k + 2]) / 3;
+          if (degrees <= centre.planet.radiusDeg - 0.5) centre.disc.push(luma);
+          else if (
+            degrees >= centre.planet.radiusDeg + 1 &&
+            degrees <= centre.planet.radiusDeg + 3.5
+          )
+            centre.ring.push(luma);
+        }
+      }
+    }
+  }
+  const median = (values) => values.sort((a, b) => a - b)[values.length >> 1];
+  for (const { planet, disc, ring } of centres) {
+    const discMean = disc.reduce((s, v) => s + v, 0) / disc.length;
+    const ringMedian = median(ring);
+    assert.ok(
+      discMean <= ringMedian + TOLERANCE,
+      `l=${planet.l} b=${planet.b}: disc mean ${discMean.toFixed(1)} vs surrounding sky ${ringMedian.toFixed(1)}`,
+    );
+  }
+});