feat(saturn): a limb law per display channel, no dashes along the limb, and Cassini's 2011 storm map as a dataset

docs/surface-preparation.md

@@ -785,7 +785,8 @@ same published law, so the limb in the app is the limb the instrument saw.
   | Venus | Minnaert, k 1.32 to 1.36 | [Pérez-Hoyos et al. 2018](https://doi.org/10.1002/2017JE005406), MESSENGER MASCS |
   | Mars | Hapke, surface only | [Vincendon 2013](https://doi.org/10.1016/j.pss.2012.12.005), OMEGA and CRISM |
   | Jupiter | Minnaert per channel | [Simon et al. 2015](https://doi.org/10.1088/0004-637X/812/1/55), OPAL |
-  | Saturn, Uranus, Neptune | Minnaert per channel | the OPAL README of each map |
+  | Saturn | Minnaert per display channel: the README's exponents weighted by the light each channel shows | the [OPAL README](https://archive.stsci.edu/missions/hlsp/opal/cycle32/saturn/hlsp_opal_hst_wfc3-uvis_saturn-2025_all_v1_readme.txt) of its map and the spectrum of [Karkoschka 1998](https://doi.org/10.1006/icar.1998.5913), weighted by [display-limb.mts](../packages/bake/authoring/saturn/display-limb.mts) |
+  | Uranus, Neptune | Minnaert per channel | the OPAL README of each map |
   | Earth | Minnaert per channel | fitted here to six [DSCOVR EPIC](https://epic.gsfc.nasa.gov/about) Level 1B frames ([fit-epic-limb.mts](../packages/bake/cli/fit-epic-limb.mts)) |
   | Moon | Hapke at 643 nm | [Sato et al. 2014](https://doi.org/10.1002/2013JE004580), the correction of the LROC WAC mosaic; w, b and h_S are medians of its PDS parameter map |
   | Ceres (dwarf planet) | Hapke at 749 nm | [Li et al. 2019](https://doi.org/10.1016/j.icarus.2018.12.038), Dawn Framing Camera |

packages/bake/authoring/saturn/display-limb.mts

@@ -0,0 +1,122 @@
+/**
+ * Saturn's limb law for each display channel.
+ *
+ * Hubble OPAL publishes one Minnaert exponent per filter, in the README of its Saturn maps. A display channel shows a
+ * range of wavelengths, so no single filter is its law. Each channel takes the README's exponents weighted by the light
+ * that channel shows: Saturn's full-disc albedo (Karkoschka 1998, PDS 1995LOW.TAB) times a 5,772 K Planck spectrum,
+ * through the CIE 1931 observer into linear sRGB, the computation behind the body's whole-disc color. Between two
+ * filters the exponent is interpolated in wavelength, and outside them the nearest one holds; that join is ours. The
+ * methane-band filters are left out: a narrow absorption band says nothing of the continuum beside it.
+ *
+ * A weighted sum of Minnaert laws is not a Minnaert law. Each channel records the exponent that fits its weighted
+ * profile best over the flood-lit disc, weighted by projected area, out to the emission angle the map's data reached,
+ * and the command prints how far that power law departs from the profile.
+ *
+ *   node packages/bake/authoring/saturn/display-limb.mts --albedo=<1995low.tab> [--write]
+ *
+ * Reads the per-filter records in src/objects/saturn/source/photometry/ and checks the spectrum against the recorded
+ * whole-disc color; --write stores the three display records beside them. Preparation only.
+ */
+import { projectRoot as checkoutProjectRoot } from '@cssearth/core/node';
+import { readFile, readdir, writeFile } from 'node:fs/promises';
+import { resolve } from 'node:path';
+import { pathToFileURL } from 'node:url';
+import { parseCieTable } from '@cssearth/bake/objects/color';
+import { readCie1931ColorMatching } from '@cssearth/bake/objects/sources';
+import { spectrumLinearSrgb } from '@cssearth/bake/objects/stellar';
+import { PHOTOMETRIC_MODEL_SCHEMA, loadPhotometricModelRecord } from '@cssearth/bake/photometry';
+
+/** The Sun's effective temperature in NASA's Sun fact sheet: the illuminant of the whole-disc color record. */
+const SOLAR_KELVIN = 5772;
+const PLANCK_H = 6.62607015e-34, LIGHT_C = 299792458, BOLTZMANN_K = 1.380649e-23;
+const VISIBLE_NM = Array.from({ length: 401 }, (_, i) => 380 + i);
+const CHANNELS = ['red', 'green', 'blue'] as const;
+
+export interface FilterExponent { readonly wavelengthNm: number; readonly coefficient: number }
+
+/** The exponent at a wavelength: linear between two filters, the nearest filter's outside them. */
+export function exponentAt(filters: readonly FilterExponent[], wavelengthNm: number): number {
+  const sorted = [...filters].sort((a, b) => a.wavelengthNm - b.wavelengthNm);
+  if (!sorted.length) throw new TypeError('A display limb law needs at least one filter.');
+  if (wavelengthNm <= sorted[0]!.wavelengthNm) return sorted[0]!.coefficient;
+  for (let i = 1; i < sorted.length; i++) {
+    const low = sorted[i - 1]!, high = sorted[i]!;
+    if (wavelengthNm <= high.wavelengthNm) return low.coefficient + (high.coefficient - low.coefficient) * (wavelengthNm - low.wavelengthNm) / (high.wavelengthNm - low.wavelengthNm);
+  }
+  return sorted.at(-1)!.coefficient;
+}
+
+/**
+ * Each display channel's flood-lit limb profile, mu^(2k - 1) weighted by the light the channel shows, and the Minnaert
+ * exponent that fits it best over the disc out to `maximumEmissionDegrees`.
+ */
+export function displayLimbExponents({ filters, light, colorMatching, maximumEmissionDegrees, steps = 1024 }: {
+  filters: readonly FilterExponent[]; light: (wavelengthNm: number) => number; colorMatching: Map<number, readonly number[]>; maximumEmissionDegrees: number; steps?: number;
+}) {
+  const whole = spectrumLinearSrgb(VISIBLE_NM, light, colorMatching), limit = Math.sin(maximumEmissionDegrees * Math.PI / 180);
+  const radii = Array.from({ length: steps }, (_, step) => (step + 0.5) / steps * limit);
+  const profiles = radii.map(radius => {
+    const mu = Math.sqrt(1 - radius * radius), lit = spectrumLinearSrgb(VISIBLE_NM, wavelength => light(wavelength) * mu ** (2 * exponentAt(filters, wavelength) - 1), colorMatching);
+    return { radius, mu, factors: [lit[0] / whole[0], lit[1] / whole[1], lit[2] / whole[2]] as const };
+  });
+  return CHANNELS.map((channel, index) => {
+    let best = { coefficient: Number.NaN, error: Infinity };
+    for (let thousandth = 300; thousandth <= 1300; thousandth++) {
+      const coefficient = thousandth / 1000;
+      let error = 0;
+      for (const { radius, mu, factors } of profiles) error += (mu ** (2 * coefficient - 1) - factors[index]!) ** 2 * radius;
+      if (error < best.error) best = { coefficient, error };
+    }
+    const departure = Math.max(...profiles.map(({ mu, factors }) => Math.abs(mu ** (2 * best.coefficient - 1) - factors[index]!)));
+    return { channel, coefficient: best.coefficient, departure, discMean: 2 / (2 * best.coefficient + 1) };
+  });
+}
+
+/** Saturn's albedo from the PDS table's eight columns: air wavelength is the second, Saturn's albedo the fifth. */
+export function readSaturnAlbedo(table: string): (wavelengthNm: number) => number {
+  const rows = table.split(/\r?\n/u).filter(line => line.trim()).map(line => { const cells = line.trim().split(/\s+/u).map(Number); return { cells: cells.length, wavelength: cells[1]!, albedo: cells[4]! }; });
+  if (rows.length < 2 || rows.some(row => row.cells !== 8 || !Number.isFinite(row.wavelength) || !Number.isFinite(row.albedo))) throw new TypeError('The albedo table is not the 1995LOW.TAB layout.');
+  if (rows[0]!.wavelength > VISIBLE_NM[0]! || rows.at(-1)!.wavelength < VISIBLE_NM.at(-1)!) throw new RangeError('The albedo table does not cover 380-780 nm.');
+  return wavelength => {
+    const next = Math.max(1, rows.findIndex(row => row.wavelength >= wavelength)), low = rows[next - 1]!, high = rows[next]!;
+    return low.albedo + (high.albedo - low.albedo) * (wavelength - low.wavelength) / (high.wavelength - low.wavelength);
+  };
+}
+
+const planck = (wavelengthNm: number) => { const metres = wavelengthNm * 1e-9; return 1 / (metres ** 5 * Math.expm1(PLANCK_H * LIGHT_C / (metres * BOLTZMANN_K * SOLAR_KELVIN))); };
+
+if (import.meta.url === pathToFileURL(process.argv[1] ?? '').href) {
+  const albedoPath = process.argv.find(argument => argument.startsWith('--albedo='))?.slice('--albedo='.length);
+  if (!albedoPath) throw new TypeError('Usage: display-limb.mts --albedo=<1995low.tab> [--write]');
+  const source = resolve(checkoutProjectRoot(import.meta.url), 'src/objects/saturn/source'), photometry = resolve(source, 'photometry');
+  // The README's continuum filters: every per-filter record beside the display records, by the wavelength its filter states.
+  const paths = (await readdir(photometry)).filter(name => /^opal-2025-minnaert-f\d+[nmw]\.json$/u.test(name)).sort();
+  const records = await Promise.all(paths.map(name => loadPhotometricModelRecord(source, `photometry/${name}`)));
+  const filters = records.map(record => {
+    const wavelengthNm = Number(/(\d+) nm/u.exec(record.filter)?.[1]), model = record.model;
+    if (!Number.isFinite(wavelengthNm) || model.family !== 'separable' || model.disk.family !== 'minnaert') throw new TypeError(`${record.id} is not a Minnaert record with a stated wavelength.`);
+    return { id: record.id, wavelengthNm, coefficient: model.disk.coefficient };
+  });
+  const ranges = records.map(record => JSON.stringify(record.fit));
+  if (new Set(ranges).size !== 1) throw new Error('The per-filter records state different fitted ranges.');
+  const fit = records[0]!.fit;
+  if (!fit.emissionDegrees) throw new TypeError('The per-filter records state no emission range.');
+  const albedo = readSaturnAlbedo(await readFile(resolve(albedoPath), 'latin1')), light = (wavelength: number) => planck(wavelength) * albedo(wavelength);
+  const colorMatching = parseCieTable((await readCie1931ColorMatching()).toString('utf8'), 3);
+  // The same spectrum gives the recorded whole-disc color; a different table or observer would not.
+  const recorded = JSON.parse(await readFile(resolve(photometry, 'karkoschka-1998-whole-disc-color.json'), 'utf8')).linearSrgb as number[];
+  const whole = spectrumLinearSrgb(VISIBLE_NM, light, colorMatching);
+  for (const channel of [1, 2]) if (Math.abs(whole[channel]! / whole[0] - recorded[channel]! / recorded[0]!) > 5e-4)
+    throw new Error(`The spectrum's ${CHANNELS[channel]}/red ratio ${(whole[channel]! / whole[0]).toFixed(4)} is not the recorded whole-disc color's ${(recorded[channel]! / recorded[0]!).toFixed(4)}.`);
+  const result = displayLimbExponents({ filters, light, colorMatching, maximumEmissionDegrees: fit.emissionDegrees[1] });
+  console.log('filters', filters.map(filter => `${filter.wavelengthNm} nm ${filter.coefficient}`).join(', '));
+  for (const { channel, coefficient, departure, discMean } of result) console.log(`${channel}: k ${coefficient.toFixed(3)}, flood-disc mean ${discMean.toFixed(4)}, largest departure from the weighted profile ${departure.toFixed(4)}`);
+  if (process.argv.includes('--write')) {
+    for (const { channel, coefficient } of result) await writeFile(resolve(photometry, `opal-2025-minnaert-display-${channel}.json`), JSON.stringify({
+      schema: PHOTOMETRIC_MODEL_SCHEMA, id: `opal-2025-minnaert-display-${channel}`, instrument: 'Hubble WFC3/UVIS (OPAL)',
+      filter: `sRGB ${channel}: the README's continuum filters weighted by Saturn's sunlit spectrum through the CIE 1931 observer`, quantity: 'radiance-factor',
+      model: { family: 'separable', disk: { family: 'minnaert', coefficient, coefficientPerDegree: 0 } }, fit,
+    }, null, 2) + '\n');
+    console.log('wrote the three display records');
+  }
+}

packages/bake/authoring/saturn/display-limb.test.mts

@@ -0,0 +1,40 @@
+import assert from 'node:assert/strict';
+import test from 'node:test';
+import { parseCieTable } from '@cssearth/bake/objects/color';
+import { readCie1931ColorMatching } from '@cssearth/bake/objects/sources';
+import { displayLimbExponents, exponentAt, readSaturnAlbedo } from './display-limb.mts';
+
+const colorMatching = parseCieTable((await readCie1931ColorMatching()).toString('utf8'), 3);
+
+test('an exponent is interpolated between two filters and held outside them', () => {
+  const filters = [{ wavelengthNm: 600, coefficient: 0.8 }, { wavelengthNm: 400, coefficient: 0.4 }];
+  assert.equal(exponentAt(filters, 380), 0.4);
+  assert.ok(Math.abs(exponentAt(filters, 450) - 0.5) < 1e-12);
+  assert.equal(exponentAt(filters, 700), 0.8);
+  assert.throws(() => exponentAt([], 500), /at least one filter/);
+});
+
+test('one filter gives every display channel its exponent', () => {
+  const result = displayLimbExponents({ filters: [{ wavelengthNm: 550, coefficient: 0.7 }], light: () => 1, colorMatching, maximumEmissionDegrees: 86.3 });
+  assert.deepEqual(result.map(channel => channel.coefficient), [0.7, 0.7, 0.7]);
+  for (const channel of result) { assert.ok(channel.departure < 1e-9); assert.ok(Math.abs(channel.discMean - 2 / 2.4) < 1e-12); }
+});
+
+test('a channel takes the exponents of the wavelengths it shows', () => {
+  const filters = [{ wavelengthNm: 450, coefficient: 0.6 }, { wavelengthNm: 650, coefficient: 0.9 }];
+  const [red, green, blue] = displayLimbExponents({ filters, light: () => 1, colorMatching, maximumEmissionDegrees: 86.3 });
+  // An sRGB primary weighs some wavelengths negatively, so a channel's exponent may lie outside the filters' range.
+  assert.ok(blue!.coefficient < green!.coefficient && green!.coefficient < red!.coefficient);
+  // The light's spectrum is part of the weight: another spectrum gives other exponents.
+  const reddened = displayLimbExponents({ filters, light: wavelength => wavelength ** 4, colorMatching, maximumEmissionDegrees: 86.3 });
+  assert.notDeepEqual(reddened.map(channel => channel.coefficient), [red!.coefficient, green!.coefficient, blue!.coefficient]);
+});
+
+test('the albedo table is read by column, past 1,000 nm too', () => {
+  const rows = [' 300.4  300.31   .0000 .2128 .2532 .5306 .6933 .0431', ' 600.0  599.83   .0010 .5000 .5000 .5000 .5000 .5000', '1049.2 1048.89   .2688 .3974 .5348 .0704 .0444 .1821'];
+  const albedo = readSaturnAlbedo(rows.join('\r\n') + '\r\n');
+  assert.ok(Math.abs(albedo(599.83) - 0.5) < 1e-12);
+  assert.ok(albedo(450) > 0.2532 && albedo(450) < 0.5);
+  assert.throws(() => readSaturnAlbedo(rows.slice(0, 2).join('\n')), /380-780 nm/);
+  assert.throws(() => readSaturnAlbedo('1 2 3\n4 5 6'), /1995LOW/);
+});

packages/bake/src/objects/geometry/ellipsoid.ts

@@ -235,3 +235,21 @@ export function planetographicRowsToMeshLatitude(data: Uint8Array, width: number
   }
   return output;
 }
+
+/** Resample an equirectangular map whose rows are samples of planetocentric latitude from +90 to -90 degrees inclusive (the first
+ * and last rows are the poles, as the Cassini ISS global maps of the PDS Atmospheres Node are gridded) onto rows evenly spaced in
+ * the mesh's own latitude: tan(planetocentric) = (b / a) tan(beta). The result has one row fewer than the samples. Rows are
+ * linearly interpolated; columns are unchanged. */
+export function planetocentricSampleRowsToMeshLatitude(data: Uint8Array, width: number, height: number, channels: number, axisRatio: number): Buffer {
+  if (!(axisRatio >= 1) || !Number.isFinite(axisRatio)) throw new RangeError('An axis ratio a / b must be finite and at least 1.');
+  if (!(height >= 2) || data.length !== width * height * channels) throw new RangeError('Map bytes do not match its dimensions.');
+  const rows = height - 1, rowBytes = width * channels, output = Buffer.alloc(rows * rowBytes);
+  for (let y = 0; y < rows; y += 1) {
+    const beta = Math.PI / 2 - (y + 0.5) / rows * Math.PI;
+    const planetocentric = Math.atan2(Math.sin(beta), axisRatio * Math.cos(beta));
+    const source = Math.max(0, Math.min(rows, (Math.PI / 2 - planetocentric) / Math.PI * rows));
+    const y0 = Math.min(rows - 1, Math.floor(source)), y1 = y0 + 1, fraction = source - y0;
+    for (let i = 0; i < rowBytes; i += 1) output[y * rowBytes + i] = Math.round(data[y0 * rowBytes + i]! * (1 - fraction) + data[y1 * rowBytes + i]! * fraction);
+  }
+  return output;
+}

packages/bake/src/objects/layers/material-composition/layered-oblate.test.mts

@@ -50,8 +50,8 @@ test("Saturn's tinted OPAL map is tied to Karkoschka's whole-disc color and keep
   const { tie: report } = await prepareSurfaceColor({ sourcePath: join(source, recipe.sources.surface), unobservedRows: recipe.surfaceUnobservedRows,
     width: parameters.planetSourceTextureWidth, height: parameters.planetSourceTextureHeight, equatorialToPolar: parameters.objectEquatorialRadiusKm / parameters.objectPolarRadiusKm,
     channelFactors: [tint.r / peak, tint.g / peak, tint.b / peak], tie });
-  assert.deepEqual(report, { reference: 'red', source: 'karkoschka-1998-whole-disc-color', measured: { green: 0.9317, blue: 0.6975 }, published: { green: 0.6759, blue: 0.5057 }, gains: [1, 0.7254, 0.725],
-    luminance: { factor: 1.2687, knee: 0.8, shoulderedTexels: 283798, shoulderedShare: 0.0684 } });
+  assert.deepEqual(report, { reference: 'red', source: 'karkoschka-1998-whole-disc-color', measured: { green: 0.9317, blue: 0.6975 }, published: { green: 0.7052, blue: 0.4573 }, gains: [1, 0.7569, 0.6556],
+    luminance: { factor: 1.2391, knee: 0.8, shoulderedTexels: 216448, shoulderedShare: 0.0522 } });
 });
 
 test("Saturn's caps sit on the oblate polar carrier, which faces out of the body at the lane's own tile size", async () => {