fix(nearby-universe): galaxies with a Hubble Cepheid distance sit at it, among their Cepheids

packages/bake/cli/prepare-catalogue-points.mts

@@ -33,7 +33,7 @@ import { readFitsHdus, binaryTable, tableColumn, numbers } from '@cssearth/bake/
 import { parseCieTable, linearToSrgb } from '@cssearth/bake/objects/color';
 import { spectrumLinearSrgb } from '@cssearth/bake/objects/stellar';
 import { readCie1931ColorMatching } from '@cssearth/bake/objects/sources';
-import { type CatalogueSizeBy, checkCatalogueSizeBy, checkCatalogueToneBy, placeGroupMembers, recipePublished, toneCataloguePalette, writeCatalogueBank } from '@cssearth/bake/volume/node';
+import { type CatalogueSizeBy, checkCatalogueSizeBy, checkCatalogueToneBy, placeGroupMembers, placeMeasuredRows, recipePublished, toneCataloguePalette, writeCatalogueBank } from '@cssearth/bake/volume/node';
 
 const [objectDirectoryArgument, id] = process.argv.slice(2);
 if (!objectDirectoryArgument || !id || !/^[a-z][a-z0-9-]*$/u.test(id)) throw new TypeError('Usage: prepare-catalogue-points.mts <object-directory> <id>');
@@ -65,7 +65,11 @@ const table = recipe.table as { path: string; bytes: number; format: 'whitespace
    * distance unit), spread along the line of sight as widely as the group spreads across the sky: `depth: 'isotropic'`
    * gives each member the sky-plane offset of the member half the group away in table order. A row's own distance is kept
    * where its group has one member here or no distance. `source` and `basis` say whose groups and why. */
-  groupDistance?: { group: Column; distance: Column; depth: 'isotropic'; source: string; basis: string } };
+  groupDistance?: { group: Column; distance: Column; depth: 'isotropic'; source: string; basis: string };
+  /** Rows a paper measures one by one: a CSV with a header, `name` (this table's name column) and `distance` in `unit`. Such a
+   * row sits at that distance on its own sight line, whatever its column, its redshift or its group says: a galaxy whose
+   * Cepheids Hubble measured is where they put it, not at its group's average. `source` and `basis` say whose measurements. */
+  measuredDistance?: { path: string; unit: 'distance-modulus'; source: string; basis: string } };
 const discPlacement = (recipe.frame as { placement?: unknown } | undefined)?.placement === 'image-layer-disc';
 // A galaxy with no disc (M87): each ICRS row at a depth drawn from a published spheroid's density along its sight line.
 const spheroidPlacement = (recipe.frame as { placement?: unknown } | undefined)?.placement === 'spheroid';
@@ -282,7 +286,7 @@ if r['disc']: xyz = []
 elif r['icrs']: xyz = SkyCoord(ra=[k[1] for k in kept] * u.deg, dec=[k[2] for k in kept] * u.deg, distance=[k[3] for k in kept] * scale, frame='icrs').cartesian.xyz.to(out).value.T
 else: xyz = SkyCoord(l=[k[1] for k in kept] * u.deg, b=[k[2] for k in kept] * u.deg, distance=[k[3] for k in kept] * scale, frame='galactic').icrs.cartesian.xyz.to(out).value.T
 json.dump({'rows': rows, 'selected': len(kept) + missing, 'missingDistance': missing, 'excluded': excluded, 'astropy': astropy.__version__,
-  'points': [[round(float(v), 4) for v in p] for p in xyz], 'groups': grouped if r['group'] else None, 'colors': [k[4] for k in kept], 'magnitudes': [k[6] for k in kept], 'bands': [k[7] for k in kept], 'redshifts': redshifts,
+  'points': [[round(float(v), 4) for v in p] for p in xyz], 'names': [k[0] for k in kept], 'groups': grouped if r['group'] else None, 'colors': [k[4] for k in kept], 'magnitudes': [k[6] for k in kept], 'bands': [k[7] for k in kept], 'redshifts': redshifts,
   'sigmas': [None if k[5] is None else round(min(k[5], 1e6), 4) for k in kept], 'sky': [[k[1], k[2]] for k in kept] if r['disc'] else None,
   'maxDistanceKpc': 0 if r['disc'] else max(k[3] for k in kept) * (.001 if r['unit'] in ('pc', 'distance-modulus', 'redshift-planck18') else 1)}, sys.stdout)`;
 const { astroqueryToolchainSync } = await import('@cssearth/telescope/node');
@@ -294,7 +298,7 @@ const run = spawnSync(toolchain.python, ['-c', python], { env: { ...process.env,
     disc: skyPlacement, filters, colorColumn: colorBy?.column ?? colorByClass?.column ?? colorByBv?.column ?? null, toneColumn: toneBy?.magnitudeColumn ?? null, nanomaggies: toneBy?.nanomaggies === true,
     bands: colorByBands ? [colorByBands.red, colorByBands.green, colorByBands.blue] : null, outUnit: frame.unit ?? 'kpc', ...(table.missingDistance === undefined ? {} : { missing: table.missingDistance }) }) });
 if (run.status !== 0) throw new Error(`Catalogue point conversion failed for ${table.path}: ${run.stderr.slice(-2000)}`);
-const converted = JSON.parse(run.stdout) as { rows: number; selected: number; missingDistance: number; excluded: number; astropy: string; points: number[][]; colors: (number | null)[];
+const converted = JSON.parse(run.stdout) as { rows: number; selected: number; missingDistance: number; excluded: number; astropy: string; points: number[][]; names: string[]; colors: (number | null)[];
   magnitudes: (number | null)[]; bands: (number[] | null)[]; redshifts: number[] | null;
   sigmas: (number | null)[]; sky: [number, number][] | null; maxDistanceKpc: number; groups: [string, string][] | null };
 /** A disc placement may spread the rows through the disc's published thickness: each row keeps its place in the disc (the
@@ -316,6 +320,10 @@ if (groupDistance) {
   }
   groupPlaced = placeGroupMembers(converted.points, converted.groups!.map(([group]) => group), group => groupDistances.get(group) ?? null);
 }
+const measuredDistance = table.measuredDistance;
+if (measuredDistance && (skyPlacement || measuredDistance.unit !== 'distance-modulus' || !measuredDistance.source || !measuredDistance.basis)) throw new TypeError(`${at('table.measuredDistance')} needs a path, unit distance-modulus, a source and a basis, and a table placed by distance.`);
+if (measuredDistance) await readFile(resolve(objectDirectory, '../../sources', `${measuredDistance.source}.json`)).catch(() => { throw new TypeError(`${at('table.measuredDistance.source')} ${measuredDistance.source} has no record in src/sources.`); });
+const measuredPlaced = measuredDistance ? placeMeasuredRows(converted.points, converted.names, await readFile(resolve(sourceDirectory, measuredDistance.path)), measuredDistance.path, frame.unit === 'Mpc' ? 1e-6 : 1e-3) : 0;
 const discThickness = (recipe.frame as { discThickness?: unknown } | undefined)?.discThickness as undefined | {
   profile: 'sech2' | 'exponential'; scaleHeightPc: number | { atCentrePc: number; perKpcPc: number }; source: string; basis: string };
 const flare = typeof discThickness?.scaleHeightPc === 'object' && discThickness.scaleHeightPc !== null ? discThickness.scaleHeightPc : null;
@@ -576,6 +584,7 @@ const bank = { schema: 'cssearth-catalogue-points@1', id, source, meaning: recip
   counts: { rows: converted.rows, selected: converted.selected, points: converted.points.length, missingDistance: converted.missingDistance,
     ...(table.exclude ? { excluded: converted.excluded } : {}), ...(bulgePlacement ? { bulge: bulgeMembers } : {}) },
   ...(bulgePlacement ? { bulge: { source: bulgePlacement.source, basis: bulgePlacement.basis } } : {}),
+  ...(measuredDistance ? { measuredDistance: { source: measuredDistance.source, basis: measuredDistance.basis, placed: measuredPlaced } } : {}),
   // Each point's group, which a merge that keeps groups first reads (merge-catalogue-points.mts).
   ...(groupDistance ? { groupDistance: { source: groupDistance.source, depth: groupDistance.depth, basis: groupDistance.basis, placed: groupPlaced },
     groups: converted.groups!.map(([group]) => group) } : {}),

packages/bake/src/volume/node/catalogue-groups.test.ts

@@ -1,7 +1,7 @@
 import { test } from 'node:test';
 import assert from 'node:assert/strict';
 import { isDeepStrictEqual } from 'node:util';
-import { placeGroupMembers, selectByShell } from './catalogue-groups.ts';
+import { placeGroupMembers, placeMeasuredRows, selectByShell } from './catalogue-groups.ts';
 
 test('a group sits at its distance, as deep as it is wide; a lone member or a group without a distance keeps its place', () => {
   // Four members around +x at 10 units, spread along y (east of +x is +y), plus a lone galaxy and a group with no distance.
@@ -37,3 +37,13 @@ test('sky bands split each shell\'s room evenly over the sky, so a densely cover
   assert.equal(kept.filter(point => point.reference[1]! < 0).length, 1);
   assert.throws(() => selectByShell(east, 10, () => 4, undefined, 0), /skyBands/);
 });
+
+test('a row a paper measures sits at that distance on its own sight line; a listed row the table lacks is refused', () => {
+  // Three galaxies in Mpc: the first a group member moved to 19.7, the second placed by its redshift at 32, the third untouched.
+  const points = [[19.7, 0, 0, 5], [0, 32, 0], [3, 4, 0]], names = ['13179', '12342712+0211163', '50063'];
+  const placed = placeMeasuredRows(points, names, 'name,distance\n13179,31.307\n12342712+0211163,30.906\n', 'cepheid-hosts.csv', 1e-6);
+  assert.equal(placed, 2);
+  assert.deepEqual(points, [[18.2558, 0, 0, 5], [0, 15.1775, 0], [3, 4, 0]]);
+  assert.throws(() => placeMeasuredRows(points, names, 'name,distance\n99999,31\n', 'cepheid-hosts.csv', 1e-6), /99999 is not among the table's rows/u);
+  assert.throws(() => placeMeasuredRows(points, names, 'PGC,DM\n13179,31.307\n', 'cepheid-hosts.csv', 1e-6), /the header is name,distance/u);
+});

packages/bake/src/volume/node/catalogue-groups.ts

@@ -3,6 +3,7 @@
  * `packages/bake/cli/prepare-catalogue-points.mts` (`table.groupDistance`) and `merge-catalogue-points.mts`
  * (`densityCap.groupsFirst`) call these.
  */
+import { gunzipSync } from 'node:zlib';
 
 /** Members of each group of two or more at the group's distance, `distance(group)` in the points' unit: each on its own
  * sight line, moved along it by the sky-plane offset (from the members' mean direction, along the axis east of it) of the
@@ -84,3 +85,30 @@ export function selectByShell<T extends ShellPoint>(order: readonly T[], width:
     return true;
   });
 }
+
+/** Rows a paper measures one by one, at that distance on their own sight line whatever their column, redshift or group says:
+ * a galaxy whose Cepheids Hubble measured is where they put it, not at its group's average (NGC 4536 was drawn 17 Mpc from
+ * its Cepheids, by its redshift, 2026-10-01). `file` is a CSV, gzipped or not, of `name,distance` with distance moduli; `perPc` turns parsecs into the
+ * points' unit. A listed row the table does not hold is a mistake in the list. Returns how many rows moved. */
+export function placeMeasuredRows(points: number[][], names: readonly string[], file: Uint8Array | string, path: string, perPc: number,
+  round = (value: number) => Math.round(value * 1e4) / 1e4): number {
+  // The tracked list is gzipped, as the tables are (a bare .csv under source/ is a download, and ignored).
+  const csv = typeof file === 'string' ? file : (path.endsWith('.gz') ? gunzipSync(file) : Buffer.from(file)).toString('utf8');
+  const [header, ...lines] = csv.trim().split('\n');
+  if (header !== 'name,distance') throw new TypeError(`${path}: the header is name,distance, not "${header}".`);
+  const wanted = new Map(lines.map(line => {
+    const [name, modulus] = line.split(','), value = Number(modulus);
+    if (!name || !modulus || !Number.isFinite(value)) throw new TypeError(`${path}: "${line}" is not a name and a distance modulus.`);
+    return [name, 10 ** (value / 5 + 1) * perPc] as const;
+  }));
+  let placed = 0;
+  for (const [index, name] of names.entries()) {
+    const distance = wanted.get(name);
+    if (distance === undefined) continue;
+    const point = points[index]!, length = Math.hypot(point[0]!, point[1]!, point[2]!);
+    for (const axis of [0, 1, 2]) point[axis] = round(point[axis]! / length * distance);
+    wanted.delete(name); placed++;
+  }
+  if (wanted.size) throw new TypeError(`${path}: ${[...wanted.keys()].join(', ')} ${wanted.size === 1 ? 'is' : 'are'} not among the table's rows.`);
+  return placed;
+}

packages/telescope-cli/authoring/nearby-universe/cepheid-hosts.mts

@@ -0,0 +1,43 @@
+#!/usr/bin/env node
+/** The galaxies of the Nearby Universe field whose distance Hubble measured from their Cepheids, as the rows the field's two
+ * galaxy tables hold them by:
+ *
+ *   node packages/telescope-cli/authoring/nearby-universe/cepheid-hosts.mts
+ *
+ * Riess et al. (2016), table 5, gives each of 19 type Ia supernova hosts a Cepheid distance modulus (`RIESS_2016_HOSTS`, the same
+ * values that place their Cepheids' own packages). SIMBAD resolves each host's position, and the row within `MATCH_ARCSEC` of it
+ * in the Cosmicflows-4 table (named by PGC number) or, when Cosmicflows-4 does not hold the galaxy, in the 2MRS table (named by
+ * its 2MASS identifier) is listed with that modulus. `packages/bake/cli/prepare-catalogue-points.mts` (`table.measuredDistance`)
+ * then places the row there, instead of at its group's average or its redshift's distance.
+ *
+ * Output: `source/galaxies/cepheid-hosts.csv.gz` and `source/galaxies-2mrs/cepheid-hosts.csv.gz`, `name,distance`. It prints each match. */
+import { readFile, writeFile } from 'node:fs/promises';
+import { resolve } from 'node:path';
+import { gunzipSync, gzipSync } from 'node:zlib';
+import { resolveSkyTarget, WORKSPACE } from '@cssearth/telescope/node';
+import { hostName, RIESS_2016_HOSTS } from '../../src/new-object/sh0es.mts';
+
+/** A galaxy's catalogue position and SIMBAD's agree to a few arcseconds; the nearest other galaxy in these tables is arcminutes away. */
+const MATCH_ARCSEC = 30;
+const source = resolve(WORKSPACE, 'src/objects/nearby-universe/source');
+const TABLES = [{ bank: 'galaxies', path: 'cf4-hyperleda.csv.gz' }, { bank: 'galaxies-2mrs', path: 'twomrs-sample.csv.gz' }] as const;
+
+const radians = (degrees: number) => degrees * Math.PI / 180;
+const separationArcsec = (ra1: number, dec1: number, ra2: number, dec2: number) => Math.acos(Math.min(1,
+  Math.sin(radians(dec1)) * Math.sin(radians(dec2)) + Math.cos(radians(dec1)) * Math.cos(radians(dec2)) * Math.cos(radians(ra1 - ra2)))) * 180 / Math.PI * 3600;
+// Both tables are name, RAJ2000, DEJ2000, … with a header (their points.json columns 1 to 3).
+const rows = await Promise.all(TABLES.map(async table => gunzipSync(await readFile(resolve(source, table.bank, table.path))).toString('utf8').trim().split('\n').slice(1)
+  .map(line => line.split(',')).filter(cells => cells[1] && cells[2]).map(cells => ({ name: cells[0]!, ra: Number(cells[1]), dec: Number(cells[2]) }))));
+
+const listed: string[][] = TABLES.map(() => []);
+for (const [host, [modulus]] of Object.entries(RIESS_2016_HOSTS)) {
+  const galaxy = hostName(host), resolved = await resolveSkyTarget(WORKSPACE, galaxy);
+  if (!resolved) throw new Error(`SIMBAD does not know ${galaxy}.`);
+  const { raDegrees, decDegrees } = resolved.target;
+  const found = rows.map(table => table.map(row => ({ row, arcsec: separationArcsec(row.ra, row.dec, raDegrees, decDegrees) })).reduce((nearest, next) => next.arcsec < nearest.arcsec ? next : nearest));
+  const bank = found.findIndex(match => match.arcsec <= MATCH_ARCSEC);
+  if (bank < 0) throw new Error(`${galaxy}: no row within ${MATCH_ARCSEC}" in ${TABLES.map(table => table.path).join(' or ')} (nearest ${found.map(match => `${match.arcsec.toFixed(0)}"`).join(', ')}).`);
+  listed[bank]!.push(`${found[bank]!.row.name},${modulus}`);
+  console.log(`${galaxy}: ${TABLES[bank]!.bank} ${found[bank]!.row.name}, ${found[bank]!.arcsec.toFixed(1)}" from SIMBAD's position, modulus ${modulus}`);
+}
+for (const [index, table] of TABLES.entries()) await writeFile(resolve(source, table.bank, 'cepheid-hosts.csv.gz'), gzipSync(`name,distance\n${listed[index]!.join('\n')}\n`, { level: 9 }));