forked from github.com/concept-collection/turing-surface
Accept a surface handed over as coefficients (reharm's export button)
A geometry can now arrive as spherical-harmonic coefficients instead of a
.m: Geometry.create with `coeffs` set re-indexes them onto the session's
band (relayoutCoeffs — truncation and zero-padding, both exact) and runs
the unchanged synthesis/metric tail. reharm's "Export to turing-surface"
opens the page with ?import and posts {lmax, mmax, Gx, Gy, Gz, provenance}
by postMessage once the page answers ready (only after boot has a live
session). The surface lands as an 'imported' dropdown entry and the run
adopts its lmax (non-stock values join the lmax dropdown). The .m editor,
bench command line and MATLAB export say "imported surface" instead of
applying. Checks: same-band import bit-identical, padding matches the
closed form on a finer grid, truncation exact on the sphere.
7 changed files+389−11
README.mdmodified+10−0View file
@@ -68,6 +68,16 @@ solver and does not disturb the run: the geometry is data whose shape in the
6868 bindings depends only on the grid, so a swap is sixteen buffer writes and the
6969 pattern carries straight on.
7070
71+A sixth kind of geometry arrives from outside: reharm (the sibling
72+spherical-harmonic fitting app) can hand a fitted surface — a scanned mesh, a
73+brain, a point cloud — to this page with its **Export to turing-surface**
74+button. The surface arrives as coefficients by `postMessage` (there is no .m;
75+`Geometry` is built from them directly, re-indexed onto whatever band is
76+selected), lands as an **imported** entry in the geometry dropdown at its own
77+lmax, and behaves like any other shape from there — except that the editor,
78+the benchmark command line and the MATLAB export, all of which need a .m, say
79+so instead.
80+
7181 A **morph** slider blends the drawn surface back to the unit sphere. The
7282 parametrization is the sphere's either way, so sweeping it shows which point
7383 went where.
src/geom/geometry.tsmodified+28−9View file
@@ -40,8 +40,9 @@ import {
4040 type RuntimeValue,
4141 } from 'numbl-src/numbl-core/runtime/types.ts';
4242 import { ShtPlan } from '../sht/sht.ts';
43-import type { ShtConfig } from '../sht/layout.ts';
43+import { relayoutCoeffs, type ShtConfig } from '../sht/layout.ts';
4444 import type { DerivPlan } from '../sht/deriv.ts';
45+import type { ImportedCoeffs } from './registry.ts';
4546 import { computeMetric, computeFluxMetric } from './metric.ts';
4647 import { toolFiles } from '../tools.ts';
4748 import { inFunction, inModel, ModelCompileError } from '../mgpu/errors.ts';
@@ -61,6 +62,13 @@ export interface GeometryOptions {
6162 params: ModelParams;
6263 /** Computes the theta/phi derivatives the inverse metric quantities need. */
6364 deriv: DerivPlan;
65+ /**
66+ * The surface as coefficients instead of a .m. When set, `source` is never
67+ * evaluated: the coefficients are re-indexed onto this plan's band
68+ * (truncation is exact projection, padding is exact) and everything from
69+ * synthesis onward runs unchanged.
70+ */
71+ coeffs?: ImportedCoeffs;
6472 }
6573
6674 export class Geometry {
@@ -189,16 +197,27 @@ export class Geometry {
189197 const { sht, cfg, source, paramNames, params, deriv } = opts;
190198 const npts = cfg.nlat * cfg.nphi;
191199
192- const { theta, phi } = gridAngles(sht, cfg);
193- const raw = evaluateShape(source, paramNames, params, theta, phi, npts);
194-
195200 // Coefficients first, then back to the grid: what the solver and the
196201 // renderer both see is the band-limited surface, not the raw .m output.
197- const [X, Y, Z] = [
198- await sht.analys(raw[0]),
199- await sht.analys(raw[1]),
200- await sht.analys(raw[2]),
201- ];
202+ // An imported surface arrives already as coefficients, so its .m stage is
203+ // a re-indexing onto this plan's band rather than an evaluation.
204+ let X: Float32Array;
205+ let Y: Float32Array;
206+ let Z: Float32Array;
207+ if (opts.coeffs) {
208+ const from = { lmax: opts.coeffs.lmax, mmax: opts.coeffs.mmax };
209+ X = relayoutCoeffs(opts.coeffs.X, from, cfg);
210+ Y = relayoutCoeffs(opts.coeffs.Y, from, cfg);
211+ Z = relayoutCoeffs(opts.coeffs.Z, from, cfg);
212+ } else {
213+ const { theta, phi } = gridAngles(sht, cfg);
214+ const raw = evaluateShape(source, paramNames, params, theta, phi, npts);
215+ [X, Y, Z] = [
216+ await sht.analys(raw[0]),
217+ await sht.analys(raw[1]),
218+ await sht.analys(raw[2]),
219+ ];
220+ }
202221 const [x, y, z] = [
203222 await sht.synth(X),
204223 await sht.synth(Y),
src/geom/registry.tsmodified+46−0View file
@@ -23,6 +23,23 @@ export interface MGeometry {
2323 params: ParamSpec[];
2424 /** MATLAB source — the shape itself. */
2525 source: string;
26+ /**
27+ * When set, the shape IS these spherical-harmonic coefficients and `source`
28+ * is a display-only stub: Geometry.create re-indexes them onto the session's
29+ * band instead of evaluating any .m. How a surface handed over by another
30+ * app (reharm's "Export to turing-surface") enters the pipeline.
31+ */
32+ coeffs?: ImportedCoeffs;
33+}
34+
35+/** A surface as coefficients of x, y, z — the same three arrays a reference
36+ * file's /geometry group holds, in the shared SHTNS layout (src/sht/layout.ts). */
37+export interface ImportedCoeffs {
38+ lmax: number;
39+ mmax: number;
40+ X: Float32Array;
41+ Y: Float32Array;
42+ Z: Float32Array;
2643 }
2744
2845 const sphere: MGeometry = {
@@ -92,6 +109,35 @@ export const defaultGeometryParams = (g: MGeometry): Params =>
92109 * none gets: the case where the solver is exact. */
93110 export const SPHERE_KEY = 'sphere';
94111
112+/** The dropdown key an imported surface lives under. Not in `mGeometries`:
113+ * the entry exists only after an import arrives, and the page owns it. */
114+export const IMPORTED_GEOMETRY_KEY = 'imported';
115+
116+/**
117+ * Wrap a handed-over surface as a geometry the rest of the app can hold like
118+ * any registry entry. No parameters — the shape is whatever the coefficients
119+ * say — and the .m slot carries a note rather than code, since the import
120+ * path in Geometry.create never evaluates it.
121+ */
122+export function makeImportedGeometry(
123+ name: string,
124+ blurb: string,
125+ coeffs: ImportedCoeffs,
126+): MGeometry {
127+ return {
128+ key: IMPORTED_GEOMETRY_KEY,
129+ label: name,
130+ blurb,
131+ params: [],
132+ source:
133+ `% ${name}\n` +
134+ `% This surface was handed over as spherical-harmonic coefficients\n` +
135+ `% (lmax ${coeffs.lmax}), so there is no .m to edit: the solver uses the\n` +
136+ `% coefficients directly, and edits here have no effect.\n`,
137+ coeffs,
138+ };
139+}
140+
95141 /** What the app and the benchmark start on. Not the sphere: this project
96142 * exists for the other shapes, and opening on the reference case would hide
97143 * the one thing it adds. */
src/main.tsmodified+155−2View file
@@ -1,5 +1,5 @@
11 import { requestShtDevice, describeAdapter } from './sht/sht.ts';
2-import { gridForLmax } from './sht/layout.ts';
2+import { gridForLmax, nlmCalc } from './sht/layout.ts';
33 import { ModelSession } from './mgpu/session.ts';
44 import { mModelByKey, presets, type MModel, type Params } from './mgpu/registry.ts';
55 import { ModelCompileError, formatFailure } from './mgpu/errors.ts';
@@ -18,6 +18,9 @@ import {
1818 mGeometryByKey,
1919 defaultGeometryParams,
2020 DEFAULT_GEOMETRY_KEY,
21+ IMPORTED_GEOMETRY_KEY,
22+ makeImportedGeometry,
23+ type ImportedCoeffs,
2124 type MGeometry,
2225 } from './geom/registry.ts';
2326 import {
@@ -283,6 +286,9 @@ let posBuf: Float32Array | null = null;
283286 * mode. While it is non-null there is no `session`: the study owns one per
284287 * variant, and the panels area is its grid. */
285288 let compareRun: CompareRun | null = null;
289+/** The surface handed over by another app (reharm), when one has arrived —
290+ * what the geometry dropdown's 'imported' entry resolves to. */
291+let importedGeometry: MGeometry | null = null;
286292 /** `session`'s spectral state as of the last (re-)seed — what "Restart"
287293 * rewinds to. Captured fresh each time a new field is actually established
288294 * (rebuild/reseed), not just once, so Restart reflects the run's current
@@ -396,7 +402,10 @@ function applyPreset(presetKey: string): void {
396402 }
397403
398404 function applyGeometryChoice(key: string): void {
399- const next = mGeometryByKey(key);
405+ // The imported entry is the page's own, not the registry's — it exists only
406+ // after a surface has been handed over (see the reharm import block below).
407+ const next =
408+ key === IMPORTED_GEOMETRY_KEY ? importedGeometry : mGeometryByKey(key);
400409 if (!next) {
401410 elErr.textContent = `No .m geometry for '${key}'`;
402411 return;
@@ -449,6 +458,10 @@ function updateCommand(): void {
449458 // equivalent is the ref checker, not the benchmark.
450459 if (compareRun?.refFile) {
451460 elCmd.textContent = `npm run ref -- --in ${compareRun.refFile.label}`;
461+ } else if (geometry.coeffs) {
462+ // An imported surface exists only as coefficients in this page, so no
463+ // command line can regenerate the run.
464+ elCmd.textContent = '(imported surface — no desktop equivalent)';
452465 } else {
453466 elCmd.textContent = formatCommand(currentSpec());
454467 }
@@ -459,6 +472,13 @@ function updateCommand(): void {
459472 * variant, the same choice `currentSpec` makes. Throws on a working copy
460473 * the export cannot parse (no init/step/shape function). */
461474 function matlabScriptText(): string {
475+ if (geometry.coeffs) {
476+ throw new Error(
477+ 'The imported surface has no .m, so the MATLAB export cannot bake it in. ' +
478+ 'Export the surface from reharm as .h5 and load its /geometry ' +
479+ 'coefficients with h5read instead.',
480+ );
481+ }
462482 const spec = currentSpec();
463483 return generateMatlabScript({
464484 model,
@@ -1738,6 +1758,118 @@ async function rebuildCompare(): Promise<void> {
17381758 await startCompare();
17391759 }
17401760
1761+// ------------------------------------------------------ import from reharm
1762+// reharm's "Export to turing-surface" button opens this page with ?import and
1763+// hands the fitted surface over by postMessage — the two apps live on
1764+// different origins, so a message between opener and opened tab is the only
1765+// channel, and structured clone carries the Float32Arrays without encoding.
1766+// Protocol: the opener pings {type:'reharm-import-ping'} until this page
1767+// answers {type:'reharm-import-ready'} (only once boot has a live session to
1768+// swap the surface into), then posts {type:'reharm-geometry', version, lmax,
1769+// mmax, Gx, Gy, Gz, name, provenance}; this page validates, installs it as
1770+// the dropdown's 'imported' entry, and acks. The payload is data, not code:
1771+// three coefficient arrays and label strings, the strings stripped of markup
1772+// before they reach any innerHTML (updateGeomNote).
1773+
1774+/** null until boot resolves; then '' when imports can proceed, or the reason
1775+ * they cannot (no WebGPU). */
1776+let importGate: string | null = null;
1777+
1778+const cleanLabel = (v: unknown, max: number): string =>
1779+ String(v ?? '')
1780+ .replace(/[<>&"'`]/g, ' ')
1781+ .replace(/\s+/g, ' ')
1782+ .trim()
1783+ .slice(0, max);
1784+
1785+/** Validate a handed-over payload into coefficients, or say what is wrong. */
1786+function importedCoeffsOf(d: Record<string, unknown>): ImportedCoeffs | string {
1787+ const version = Number(d.version);
1788+ if (version !== 1) return `unknown import payload version '${String(d.version)}'`;
1789+ const lmax = Number(d.lmax);
1790+ const mmax = Number(d.mmax);
1791+ if (!Number.isInteger(lmax) || lmax < 1 || lmax > 4096) return `bad lmax '${String(d.lmax)}'`;
1792+ if (!Number.isInteger(mmax) || mmax < 0 || mmax > lmax) return `bad mmax '${String(d.mmax)}'`;
1793+ const nlm = nlmCalc(lmax, mmax);
1794+ const arrays: Float32Array[] = [];
1795+ for (const name of ['Gx', 'Gy', 'Gz']) {
1796+ const q = d[name];
1797+ if (!(q instanceof Float32Array)) return `'${name}' is not a Float32Array`;
1798+ if (q.length !== 2 * nlm) return `'${name}' has ${q.length} values, expected 2*nlm = ${2 * nlm}`;
1799+ if (!q.every(Number.isFinite)) return `'${name}' has non-finite values`;
1800+ arrays.push(q);
1801+ }
1802+ return { lmax, mmax, X: arrays[0], Y: arrays[1], Z: arrays[2] };
1803+}
1804+
1805+/** Adopt the surface: it becomes (or replaces) the dropdown's imported entry
1806+ * and is selected. The run also adopts the surface's own band — truncating a
1807+ * fit to a coarser default would silently smooth the shape that was just
1808+ * exported — so the lmax dropdown gains the exact value when it is not a
1809+ * stock one, and a band change goes through the full rebuild (new grid)
1810+ * rather than the in-place swap. */
1811+function installImportedGeometry(g: MGeometry): void {
1812+ importedGeometry = g;
1813+ let opt = elGeometry.querySelector<HTMLOptionElement>(
1814+ `option[value="${IMPORTED_GEOMETRY_KEY}"]`,
1815+ );
1816+ if (!opt) {
1817+ opt = document.createElement('option');
1818+ opt.value = IMPORTED_GEOMETRY_KEY;
1819+ elGeometry.append(opt);
1820+ }
1821+ opt.textContent = g.label;
1822+ elGeometry.value = IMPORTED_GEOMETRY_KEY;
1823+ applyGeometryChoice(IMPORTED_GEOMETRY_KEY);
1824+ const lmax = g.coeffs!.lmax;
1825+ if (Number(elLmax.value) !== lmax) {
1826+ const options = [...elLmax.options];
1827+ if (!options.some((o) => Number(o.value) === lmax)) {
1828+ const o = document.createElement('option');
1829+ o.value = String(lmax);
1830+ o.textContent = String(lmax);
1831+ elLmax.add(o, options.find((x) => Number(x.value) > lmax) ?? null);
1832+ }
1833+ elLmax.value = String(lmax);
1834+ void rebuild();
1835+ } else {
1836+ viewChange = viewChange.then(() => applyGeometry());
1837+ }
1838+}
1839+
1840+window.addEventListener('message', (ev: MessageEvent) => {
1841+ const d = ev.data as Record<string, unknown> | null;
1842+ if (!d || typeof d !== 'object' || typeof d.type !== 'string') return;
1843+ const reply = (msg: Record<string, unknown>): void => {
1844+ (ev.source as Window | null)?.postMessage(msg, { targetOrigin: ev.origin });
1845+ };
1846+ if (d.type === 'reharm-import-ping') {
1847+ // Silence before boot resolves is fine — the opener keeps pinging.
1848+ if (importGate === '') reply({ type: 'reharm-import-ready' });
1849+ else if (importGate) reply({ type: 'reharm-import-error', message: importGate });
1850+ return;
1851+ }
1852+ if (d.type !== 'reharm-geometry') return;
1853+ if (importGate !== '') {
1854+ reply({ type: 'reharm-import-error', message: importGate ?? 'the page is still starting up' });
1855+ return;
1856+ }
1857+ const coeffs = importedCoeffsOf(d);
1858+ if (typeof coeffs === 'string') {
1859+ reply({ type: 'reharm-import-error', message: coeffs });
1860+ return;
1861+ }
1862+ const p = (d.provenance ?? {}) as Record<string, unknown>;
1863+ const name = `imported: ${cleanLabel(d.name, 40) || 'surface'}`;
1864+ const blurb =
1865+ `Handed over by reharm (lmax ${coeffs.lmax}` +
1866+ (p.map ? `, map ${cleanLabel(p.map, 60)}` : '') +
1867+ (p.filter ? `, filter ${cleanLabel(p.filter, 60)}` : '') +
1868+ `).`;
1869+ installImportedGeometry(makeImportedGeometry(name, blurb, coeffs));
1870+ reply({ type: 'reharm-import-ack' });
1871+});
1872+
17411873 // ---------------------------------------------------------------- boot
17421874 async function boot(): Promise<void> {
17431875 enterMode('simulate');
@@ -1758,6 +1890,8 @@ async function boot(): Promise<void> {
17581890 elErr.textContent =
17591891 `WebGPU is not available (${e instanceof Error ? e.message : e}). ` +
17601892 `Use a WebGPU-capable browser such as Chrome or Edge.`;
1893+ importGate = 'WebGPU is not available in the opened tab';
1894+ announceImportGate();
17611895 return;
17621896 }
17631897 device.lost.then((info) => {
@@ -1766,6 +1900,25 @@ async function boot(): Promise<void> {
17661900 }
17671901 });
17681902 await rebuild();
1903+ // Imports wait for a live run to swap the surface into; a tab opened by
1904+ // reharm's export button (?import) is told without waiting to be pinged.
1905+ importGate = '';
1906+ announceImportGate();
1907+}
1908+
1909+/** Tell the opener that navigated here with ?import how imports stand. Only
1910+ * a signal — no data — so the opener's origin need not be known. */
1911+function announceImportGate(): void {
1912+ if (!window.opener || !new URLSearchParams(location.search).has('import')) return;
1913+ const msg =
1914+ importGate === ''
1915+ ? { type: 'reharm-import-ready' }
1916+ : { type: 'reharm-import-error', message: importGate };
1917+ try {
1918+ (window.opener as Window).postMessage(msg, { targetOrigin: '*' });
1919+ } catch {
1920+ // The opener may already be gone; the ping path covers the rest.
1921+ }
17691922 }
17701923
17711924 void boot();
src/mgpu/session.tsmodified+2−0View file
@@ -155,6 +155,7 @@ export class ModelSession {
155155 paramNames: geometryModel.params.map((p) => p.key),
156156 params: geometryParams,
157157 deriv,
158+ coeffs: geometryModel.coeffs,
158159 });
159160 const gpu = await GpuModel.create({
160161 device,
@@ -211,6 +212,7 @@ export class ModelSession {
211212 paramNames: geometryModel.params.map((p) => p.key),
212213 params,
213214 deriv: this.#deriv,
215+ coeffs: geometryModel.coeffs,
214216 });
215217 this.#geometry = next;
216218 this.#geometryModel = geometryModel;
src/sht/layout.tsmodified+27−0View file
@@ -46,6 +46,33 @@ export function isPowerOfTwo(n: number): boolean {
4646 return n > 0 && (n & (n - 1)) === 0;
4747 }
4848
49+/**
50+ * Re-index coefficients from one (lmax, mmax) layout into another, dropping
51+ * the modes the target cannot hold and zero-filling the ones the source does
52+ * not have. Both truncation and padding are exact spectral operations —
53+ * truncation is orthogonal projection onto the smaller band, padding changes
54+ * nothing — so a surface carried as coefficients can be moved onto any plan.
55+ */
56+export function relayoutCoeffs(
57+ q: Float32Array,
58+ from: { lmax: number; mmax: number },
59+ to: { lmax: number; mmax: number },
60+): Float32Array {
61+ if (from.lmax === to.lmax && from.mmax === to.mmax) return q;
62+ const out = new Float32Array(2 * nlmCalc(to.lmax, to.mmax));
63+ const mTop = Math.min(from.mmax, to.mmax);
64+ const lTop = Math.min(from.lmax, to.lmax);
65+ for (let m = 0; m <= mTop; m++) {
66+ for (let l = m; l <= lTop; l++) {
67+ const src = 2 * lmIndex(from.lmax, l, m);
68+ const dst = 2 * lmIndex(to.lmax, l, m);
69+ out[dst] = q[src];
70+ out[dst + 1] = q[src + 1];
71+ }
72+ }
73+ return out;
74+}
75+
4976 /** Grid sizes for a given lmax, dealiased for a reaction of polynomial degree
5077 * `pdeg` (the rule from websph's reference implementation):
5178 * nlat >= ((pdeg+1)*lmax+1)/2, nphi >= (pdeg+1)*lmax+1. nphi is rounded up to
test/geometryChecks.tsmodified+121−0View file
@@ -282,6 +282,127 @@ export async function geometryChecks(
282282 sht.destroy();
283283 }
284284
285+ // ---- a surface handed over as coefficients is the same surface -----------
286+ // The import path (reharm's "Export to turing-surface"): Geometry.create
287+ // with `coeffs` set evaluates no .m — the coefficients are re-indexed onto
288+ // the session's band and everything from synthesis onward runs unchanged.
289+ // Checked in all three band relations: same (must reproduce the original
290+ // geometry, arrays and all), wider (zero-padding; the surface must still
291+ // match its closed form on the finer plan's own grid), and narrower
292+ // (truncation; exact when the content fits the smaller band, as the
293+ // degree-1 sphere does).
294+ {
295+ const importGeometry = async (
296+ coeffs: { lmax: number; mmax: number; X: Float32Array; Y: Float32Array; Z: Float32Array },
297+ lmax: number,
298+ ) => {
299+ const { nlat, nphi } = gridForLmax(lmax, 3);
300+ const cfg = { lmax, mmax: lmax, nlat, nphi };
301+ const sht = await ShtPlan.create(device, cfg);
302+ const deriv = await DerivPlan.create(device, sht);
303+ const geometry = await Geometry.create({
304+ sht, cfg, source: '', paramNames: [], params: {}, deriv, coeffs,
305+ });
306+ return { sht, deriv, cfg, geometry };
307+ };
308+
309+ const peanut = await buildGeometry(device, 'peanut');
310+ const coeffs = {
311+ lmax: peanut.cfg.lmax,
312+ mmax: peanut.cfg.mmax,
313+ X: peanut.geometry.X,
314+ Y: peanut.geometry.Y,
315+ Z: peanut.geometry.Z,
316+ };
317+
318+ // Same band: the identical pipeline on identical inputs, so the surface
319+ // and both metric formulations must come back numerically unchanged.
320+ const same = await importGeometry(coeffs, peanut.cfg.lmax);
321+ let worst = 0;
322+ const pairs: [Float32Array, Float32Array][] = [
323+ [same.geometry.x, peanut.geometry.x],
324+ [same.geometry.y, peanut.geometry.y],
325+ [same.geometry.z, peanut.geometry.z],
326+ [same.geometry.p1, peanut.geometry.p1],
327+ [same.geometry.p2, peanut.geometry.p2],
328+ [same.geometry.q2, peanut.geometry.q2],
329+ [same.geometry.r, peanut.geometry.r],
330+ [same.geometry.Vtx, peanut.geometry.Vtx],
331+ [same.geometry.Vpx, peanut.geometry.Vpx],
332+ ];
333+ for (const [a, b] of pairs) {
334+ for (let i = 0; i < a.length; i++) {
335+ worst = Math.max(worst, Math.abs(a[i] - b[i]));
336+ }
337+ }
338+ check(
339+ 'import: coefficients rebuild the identical geometry on the same band',
340+ worst <= 1e-7,
341+ `max |difference| = ${worst.toExponential(2)} over surface and metric`,
342+ );
343+ same.deriv.destroy();
344+ same.sht.destroy();
345+
346+ // Wider band: zero-padded onto a finer plan, then checked against the
347+ // .m's own closed form on that plan's grid — points the original build
348+ // never evaluated (Gauss latitudes do not nest).
349+ const p = defaultGeometryParams(peanut.g);
350+ const peanutRadius = (ct: number): number => {
351+ const st2 = Math.max(0, 1 - ct * ct);
352+ const r = 1 - p.waist * st2;
353+ return r * Math.hypot(Math.sqrt(st2), (1 + p.stretch) * ct);
354+ };
355+ const wide = await importGeometry(coeffs, 2 * peanut.cfg.lmax + 1);
356+ let dr = 0;
357+ for (let i = 0; i < wide.cfg.nlat; i++) {
358+ const want = peanutRadius(wide.sht.cosTheta[i]);
359+ for (let j = 0; j < wide.cfg.nphi; j++) {
360+ const k = i * wide.cfg.nphi + j;
361+ const got = Math.hypot(wide.geometry.x[k], wide.geometry.y[k], wide.geometry.z[k]);
362+ dr = Math.max(dr, Math.abs(got - want));
363+ }
364+ }
365+ check(
366+ 'import: zero-padding onto a wider band keeps the surface',
367+ dr < 1e-3,
368+ `max |dr| = ${dr.toExponential(2)} at lmax ${wide.cfg.lmax}`,
369+ );
370+ wide.deriv.destroy();
371+ wide.sht.destroy();
372+ peanut.deriv.destroy();
373+ peanut.sht.destroy();
374+
375+ // Narrower band: truncation is orthogonal projection, so a surface whose
376+ // content fits the smaller band comes through exactly — the sphere is
377+ // degree 1 and must still have radius 1 after import at any lmax.
378+ const sphere = await buildGeometry(device, SPHERE_KEY);
379+ const narrow = await importGeometry(
380+ {
381+ lmax: sphere.cfg.lmax,
382+ mmax: sphere.cfg.mmax,
383+ X: sphere.geometry.X,
384+ Y: sphere.geometry.Y,
385+ Z: sphere.geometry.Z,
386+ },
387+ Math.floor(sphere.cfg.lmax / 2),
388+ );
389+ let radiusErr = 0;
390+ for (let k = 0; k < narrow.geometry.x.length; k++) {
391+ const r = Math.hypot(narrow.geometry.x[k], narrow.geometry.y[k], narrow.geometry.z[k]);
392+ radiusErr = Math.max(radiusErr, Math.abs(r - 1));
393+ }
394+ check(
395+ 'import: truncation onto a narrower band is exact when the content fits',
396+ radiusErr < 2e-3,
397+ `sphere from lmax ${sphere.cfg.lmax} coefficients at lmax ${narrow.cfg.lmax}: ` +
398+ `max |r - 1| = ${radiusErr.toExponential(2)}`,
399+ );
400+ narrow.deriv.destroy();
401+ narrow.sht.destroy();
402+ sphere.deriv.destroy();
403+ sphere.sht.destroy();
404+ }
405+
285406 // ---- the unrolled loop: more ops, identical answer ----------------------
286407 {
287408 const model = mModelByKey('schnakenberg')!;