[Written by Claude] GW150914 + pink noise = Still.Pink

# gw150914.py
# Gets the GW150914 waveform (LIGO Hanford, numerical relativity)
# from GWOSC, and writes every 6th sample to gw150914-h1.json,
# next to this file.
import json
from pathlib import Path
from urllib.request import urlopen

URL = "https://gwosc.org/GW150914data/P150914/fig1-waveform-H.txt"

text = urlopen(URL).read().decode()
rows = [line.split() for line in text.splitlines()
        if line.strip() and not line.startswith("#")]
strain = [round(float(h), 4) for _, h in rows[::6]]  # 574 points

out = Path(__file__).with_name("gw150914-h1.json")
out.write_text(json.dumps({"dt": 6 / 16384, "strain": strain}) + "\n")

// signal.ts
// Pink noise for the Still.Pink waveform. Runs at build time.

// FNV-1a: a string to a 32-bit seed.
export function hash(text: string) {
	let h = 2166136261;
	for (let i = 0; i < text.length; i++) {
		h = Math.imul(h ^ text.charCodeAt(i), 16777619);
	}
	return h >>> 0;
}

// mulberry32: a tiny seeded random number generator, in [0, 1).
export function mulberry32(a: number) {
	return () => {
		a = (a + 0x6d2b79f5) | 0;
		let t = Math.imul(a ^ (a >>> 15), 1 | a);
		t = (t + Math.imul(t ^ (t >>> 7), 61 | t)) ^ t;
		return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
	};
}

// Box–Muller: uniform numbers to standard normal ones.
export function gaussian(random: () => number) {
	return () => {
		const u = 1 - random(); // in (0, 1], so log(u) is finite
		const v = random();
		return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v);
	};
}

type Wave = { f: number; a: number; b: number };

// Waves a·cos + b·sin at f = 1/span, 2/span, ..., fMax cycles per
// window, with a and b drawn from Normal(0, 1/√f), so power ∝ 1/f.
function pinkWaves(seed: string, span: number, fMax: number) {
	const normal = gaussian(mulberry32(hash(seed)));
	const waves: Wave[] = [];
	for (let k = 1; k <= span * fMax; k++) {
		const f = k / span;
		const std = 1 / Math.sqrt(f);
		waves.push({ f, a: normal() * std, b: normal() * std });
	}
	return waves;
}

// The noise at position x, in windows.
function noiseAt(waves: Wave[], x: number) {
	let y = 0;
	for (const { f, a, b } of waves) {
		const t = 2 * Math.PI * f * x;
		y += a * Math.cos(t) + b * Math.sin(t);
	}
	return y;
}

// Pink noise, `span` windows of n samples. Every wave fits a whole
// number of cycles in it, so it loops.
export function pinkNoiseLoop(seed: string, n: number, span = 8) {
	const waves = pinkWaves(seed, span, 64);
	const L = span * n;
	return Array.from({ length: L }, (_, i) => noiseAt(waves, i / n));
}

// n samples of pink noise, starting `shift` samples into the loop.
// Same as a slice of pinkNoiseLoop(), but computes only those n.
export function pinkNoise(seed: string, n: number, shift = 0, span = 8) {
	const waves = pinkWaves(seed, span, 64);
	const L = span * n;
	const start = ((shift % L) + L) % L;
	const ys = Array.from({ length: n }, (_, i) =>
		noiseAt(waves, ((start + i) % L) / n),
	);
	return standardize(ys);
}

// n samples of a precomputed loop, starting at `shift`, standardized.
export function windowOf(loop: number[], n: number, shift: number) {
	const L = loop.length;
	const start = ((shift % L) + L) % L;
	const ys = Array.from({ length: n }, (_, i) => loop[(start + i) % L]);
	return standardize(ys);
}

// A GPS time (s) as a shift in samples of spacing dt (s).
export function gpsToShift(gps: number, dt: number) {
	return Math.round(gps / dt);
}

// Zero mean, unit standard deviation.
export function standardize(ys: number[]) {
	const mean = ys.reduce((a, y) => a + y, 0) / ys.length;
	const v = ys.reduce((a, y) => a + (y - mean) ** 2, 0) / ys.length;
	return ys.map((y) => (y - mean) / Math.sqrt(v));
}

// Scale so the biggest swing is 1.
export function scaleToPeak(ys: number[]) {
	const peak = Math.max(...ys.map(Math.abs));
	return ys.map((y) => y / peak);
}

// still-pink.ts
// GW150914 + pink noise = Still.Pink.
import gw from './gw150914-h1.json' with { type: 'json' };
import { gpsToShift, pinkNoise, scaleToPeak } from './signal.ts';

export const SEED = '150914'; // GW150914's date
export const LEVEL = 0.2; // noise std; the chirp's peak is 1
// 2026-11-13 00:05:49 UTC, the day this blog launched,
// nudged 5.9 ms later so the sum ends on a climb.
export const GPS_TIME = 1478563567.0059;

export function stillPink(seed = SEED, level = LEVEL, gps = GPS_TIME) {
	const chirp = scaleToPeak(gw.strain); // peak = 1
	const shift = gpsToShift(gps, gw.dt);
	const noise = pinkNoise(seed, chirp.length, shift)
		.map((x) => x * level);
	const sum = chirp.map((h, i) => h + noise[i]);
	// One y-scale for all three, big enough for each of them.
	const lim = Math.max(...[...chirp, ...noise, ...sum].map(Math.abs));
	return { chirp, noise, sum, range: [-lim, lim] as [number, number] };
}

// The waveform itself, computed once and shared by every page.
export const waveform = stillPink();