Climate Pair

Sixty-eight years of carbon, ocean heat, ENSO and storms, heard from the Pacific and from the Atlantic

1958–2025 · one year per bar · two pieces on one grammar

Two pieces built from the global climate record since the Keeling curve began in 1958. One bar is one year: seven seconds, four beats for the seasons, twelve pulses for the months. The things that only rise, like carbon dioxide, ocean heat and sea-surface temperature, live in the low drone and the strings, and the whole mix gets about 4 dB louder from the first decade to the last. Storms are the deep drums. ENSO sets the chords. The two pieces share their global layers and differ in their ocean.

Status: a first render. These pieces were built and balanced by measurement, and checked with the plots below, but they have not yet been refined by ear. Graham chose the data and gave direction; the musical choices are Claude Code’s first attempt. See How this was made.

Pacific

Key of D · 8:10

The Pacific Decadal Oscillation sets the mode. Typhoons in the West Pacific sound on the left and hurricanes in the Northeast Pacific on the right, as on a map. Each spring the Kettle River, which crosses the BC–Washington border near Grand Forks, blooms at the date of its annual peak.

Atlantic

Key of A · 8:10

The Atlantic Multidecadal Oscillation sets the mode, and winter wind follows the North Atlantic Oscillation. Atlantic hurricanes are the drums. From November 1978, when satellites began measuring it, Arctic sea ice enters as a high glass voice that breathes with the seasons and thins over the decades.

Watch and listen

How to read the spiral. Each ring is a year: 1958 at the centre, 2025 at the edge, January at the top, running clockwise. The ring is drawn as you hear it.

• Colour is the global ocean surface temperature anomaly, from cool blue to ember.
• Thickness is ENSO strength: the ring swells in strong El Niño and La Niña months.
• Sparks are storms at their peak, larger for higher categories. Each one flares as its drum sounds, and a ring marks a landfall. In the Pacific piece, West Pacific typhoons are pale blue and Northeast Pacific hurricanes pale amber.
• Pacific: green blooms mark the Kettle River’s annual peak, larger for bigger years. Atlantic: from November 1978, a white line along the inside of each ring is Arctic sea-ice extent. It swells each winter and thins each September.
• Grey marks are volcanic eruptions. A red halo appears while the high violin plays (global temperature above +0.8 °C). The glow at the centre follows the contrabass (ocean heat), and the twinkles at the playhead follow the shimmer (sunspots).

The strip under the spiral is a timeline: click or drag it to move through the years. The labels below it light up as each instrument plays. Tick both side by side to watch the two oceans at once while you listen to one. In the big El Niño years, the Pacific sparks crowd while the Atlantic ones thin.

The spiral and strip are drawn from the same data as the music (scripts/export_climate_viz.py). Sunspot numbers and ocean heat content are not published as values (see Sources and licences), so those two layers appear only as the loudness of their instruments.

How to listen

Each layer is tagged: R retrieved (a measured record), D derived (computed from a record, with a choice made along the way), S simulated. There is no simulated layer in these two pieces.

What you hear	Data
Drone (synth and pipe-organ pedal): louder and more open as it rises, and it breathes once a year	Mauna Loa CO₂, deseasonalized; the seasonal cycle is the Northern Hemisphere’s plants taking in and releasing carbon	R
Cellos on chord tones, one voice early and three by the late 1990s	Global ocean surface temperature anomaly	R
Contrabass below the cellos: the deep ocean filling, root then root and fifth	Ocean heat content, 0–2000 m (IAP v4.2)	R
Chords, two per year: El Niño lifts the harmony (VI, or III in a strong event), La Niña sinks it (iv); neutral years circle i–VII–i–v	Oceanic Niño Index	R
Monthly pulse (piano in the Pacific, harp in the Atlantic): louder in strong ENSO, higher in El Niño, lower in La Niña	Oceanic Niño Index	R
Mode and pad: Dorian in the warm phase, Aeolian in the cool; the pad brightens with the index	PDO or AMO, 3-year mean, with a small dead band so the mode doesn’t flicker	D
Storms: tropical storm = a rubbed bass drum; hurricane = tuned timpani; major hurricane = timpani, concert bass drum and a sub-bass boom that briefly ducks the drone; Category 5 = tam-tam; landfall = cymbal; a low rumble follows each season’s accumulated energy (ACE)	NOAA HURDAT2 (Atlantic, Northeast Pacific) and IBTrACS (West Pacific) best tracks, placed at each storm’s peak	R
Shimmer: a high pad and glockenspiel, fullest at each solar maximum	Monthly sunspot number	R
Heat: a held high violin once the 12-month global mean passes +0.8 °C	NASA GISTEMP v4	R
Volcanoes: a gong, and the drone darkens for about two years	Eruption dates; the cooling envelope is a chosen shape	D
Pacific: spring bloom of harp and piano, larger for bigger years	Kettle River near Laurier, WA: the date and size (rank) of each year’s peak daily flow	R
Atlantic: winter wind, stronger when the NAO is positive	North Atlantic Oscillation, December–March	R
Atlantic: glass voice: its loudness follows each month’s ice extent; its number of voices follows that year’s September minimum	Arctic sea-ice extent, NSIDC Sea Ice Index	R

The pair is a conversation. ENSO runs through both. In the big El Niño years (1982–83, 1997–98, 2015–16) the Pacific storm drums crowd and the Atlantic ones thin, a well-documented link between ENSO and Atlantic hurricane activity.

Check plots

Each year’s data (blue) against each instrument’s loudness (orange). A layer should speak when its data says it should.

Pacific

Atlantic

Listening guide

Every timed event in each score, generated from the render.

NotePacific: 73 events

Time	What you hear
0:00	Strong El Niño (1958-01, ONI +1.7)
0:05	Category 5 (WPAC): Ida (Sep 1958, 175 kt) - tam-tam
0:11	Category 5 (NEPAC): Patsy (Sep 1959, 150 kt) - tam-tam
0:21	PDO phase turns cool (1961) - Aeolian
0:25	Category 5 (WPAC): Nancy (Sep 1961, 185 kt) - tam-tam
0:26	Category 5 (WPAC): Violet (Oct 1961, 180 kt) - tam-tam
0:36	Eruption: Agung, Bali (Mar 1963) - gong; the drone darkens as the aerosols cool the planet
0:54	Strong El Niño (1965-11, ONI +1.8)
0:56	1966: routine satellite coverage begins - fewer storms before this is partly fewer observed
0:59	Category 5 (WPAC): Kit (Jun 1966, 170 kt) - tam-tam
1:15	Solar maximum (1968-11) - the shimmer is fullest
1:43	Strong El Niño (1972-11, ONI +1.8)
1:48	Category 5 (NEPAC): Ava (Jun 1973, 140 kt) - tam-tam
1:51	Strong La Niña (1973-12, ONI -2.0)
2:05	Strong La Niña (1975-12, ONI -1.6)
2:33	Solar maximum (1979-12) - the shimmer is fullest
2:34	PDO phase turns warm (1980) - Dorian
2:49	Eruption: El Chichon, Mexico (Apr 1982) - gong; the drone darkens as the aerosols cool the planet
2:55	Strong El Niño (1983-01, ONI +2.1)
3:28	CO₂ passes 350 ppm (1987-11) - the drone is louder and more open
3:36	Strong La Niña (1988-12, ONI -1.9)
3:37	PDO phase turns cool (1989) - Aeolian
3:39	Solar maximum (1989-06) - the shimmer is fullest
3:54	Eruption: Pinatubo, Philippines (Jun 1991) - gong; the drone darkens as the aerosols cool the planet
3:58	PDO phase turns warm (1992) - Dorian
3:58	Strong El Niño (1992-01, ONI +1.5)
4:15	Category 5 (NEPAC): Emilia (Jul 1994, 140 kt) - tam-tam
4:15	Category 5 (NEPAC): Gilma (Jul 1994, 140 kt) - tam-tam
4:16	Category 5 (NEPAC): John (Aug 1994, 150 kt) - tam-tam
4:35	Global SST anomaly passes +0.5 C (1997-06) - cello in three voices
4:36	Most active WPAC+NEPAC season by ACE: 1997 (840)
4:37	Category 5 (NEPAC): Guillermo (Aug 1997, 140 kt) - tam-tam
4:37	Category 5 (NEPAC): Oliwa (Sep 1997, 140 kt) - tam-tam
4:37	Category 5 (NEPAC): Linda (Sep 1997, 160 kt) - tam-tam
4:39	Strong El Niño (1997-12, ONI +2.4)
4:39	Category 5 (NEPAC): Paka (Dec 1997, 160 kt) - tam-tam
4:40	PDO phase turns cool (1998) - Aeolian
4:53	Strong La Niña (1999-12, ONI -1.5)
5:09	Solar maximum (2002-03) - the shimmer is fullest
5:11	Category 5 (NEPAC): Elida (Jul 2002, 140 kt) - tam-tam
5:12	Category 5 (NEPAC): Hernan (Sep 2002, 140 kt) - tam-tam
5:13	Category 5 (NEPAC): Kenna (Oct 2002, 145 kt) - tam-tam
5:40	Category 5 (NEPAC): Ioke (Aug 2006, 140 kt) - tam-tam
5:50	Strong La Niña (2008-01, ONI -1.8)
5:52	Ocean heat content (0-2000 m) past the halfway mark of its rise (2008-06) - the contrabass adds a fifth
6:02	Category 5 (NEPAC): Rick (Oct 2009, 155 kt) - tam-tam
6:03	Strong El Niño (2009-12, ONI +1.5)
6:07	Category 5 (NEPAC): Celia (Jun 2010, 140 kt) - tam-tam
6:09	Strong La Niña (2010-11, ONI -1.6)
6:30	Category 5 (WPAC): Haiyan (Nov 2013, 170 kt) - tam-tam
6:33	Solar maximum (2014-04) - the shimmer is fullest
6:36	Category 5 (NEPAC): Genevieve (Aug 2014, 140 kt) - tam-tam
6:36	Category 5 (NEPAC): Marie (Aug 2014, 140 kt) - tam-tam
6:37	Global temperature passes +0.8 C vs 1951-80 (2014-10) - a high violin enters
6:39	PDO phase turns warm (2015) - Dorian
6:40	CO₂ passes 400 ppm (2015-03) - the drone is louder and more open
6:44	Category 5 (NEPAC): Patricia (Oct 2015, 185 kt) - tam-tam
6:45	Strong El Niño (2015-12, ONI +2.6)
6:50	Category 5 (WPAC): Meranti (Sep 2016, 170 kt) - tam-tam
7:00	PDO phase turns cool (2018) - Aeolian
7:02	The Kettle’s biggest peak in the span: 1317 m3/s (2018-05-11) - the fullest harp bloom
7:04	Category 5 (NEPAC): Lane (Aug 2018, 140 kt) - tam-tam
7:05	Category 5 (NEPAC): Walaka (Oct 2018, 140 kt) - tam-tam
7:05	Category 5 (NEPAC): Willa (Oct 2018, 140 kt) - tam-tam
7:19	Category 5 (WPAC): Goni (Oct 2020, 170 kt) - tam-tam
7:23	Category 5 (WPAC): Surigae (Apr 2021, 170 kt) - tam-tam
7:28	Eruption: Hunga Tonga-Hunga Ha’apai (Jan 2022) - gong; mostly water vapour, not sulphate: the drone barely darkens
7:36	CO₂ passes 420 ppm (2023-04) - the drone is louder and more open
7:39	Category 5 (NEPAC): Jova (Sep 2023, 140 kt) - tam-tam
7:40	Category 5 (NEPAC): Otis (Oct 2023, 145 kt) - tam-tam
7:41	Strong El Niño (2023-12, ONI +2.0)
7:47	Category 5 (NEPAC): Kristy (Oct 2024, 140 kt) - tam-tam
7:49	WPAC best tracks end Dec 2024 (not yet in the archive) - that basin’s drums fall silent

NoteAtlantic: 71 events

Time	What you hear
0:00	Strong El Niño (1958-01, ONI +1.7)
0:26	Category 5 (ATL): Esther (Sep 1961, 140 kt) - tam-tam
0:26	Category 5 (ATL): Hattie (Oct 1961, 145 kt) - tam-tam
0:36	Eruption: Agung, Bali (Mar 1963) - gong; the drone darkens as the aerosols cool the planet
0:42	AMO phase turns cool (1964) - Aeolian
0:54	Strong El Niño (1965-11, ONI +1.8)
0:56	1966: routine satellite coverage begins - fewer storms before this is partly fewer observed
1:01	Category 5 (ATL): Inez (Sep 1966, 140 kt) - tam-tam
1:08	Category 5 (ATL): Beulah (Sep 1967, 140 kt) - tam-tam
1:15	Solar maximum (1968-11) - the shimmer is fullest
1:21	Category 5 (ATL): Camille (Aug 1969, 150 kt) - tam-tam
1:43	Strong El Niño (1972-11, ONI +1.8)
1:51	Strong La Niña (1973-12, ONI -2.0)
2:05	Strong La Niña (1975-12, ONI -1.6)
2:17	Category 5 (ATL): Anita (Sep 1977, 150 kt) - tam-tam
2:25	Satellites begin measuring Arctic sea ice (1978-11) - a high glass voice fades in
2:31	Category 5 (ATL): David (Aug 1979, 150 kt) - tam-tam
2:33	Solar maximum (1979-12) - the shimmer is fullest
2:38	Category 5 (ATL): Allen (Aug 1980, 165 kt) - tam-tam
2:49	Eruption: El Chichon, Mexico (Apr 1982) - gong; the drone darkens as the aerosols cool the planet
2:55	Strong El Niño (1983-01, ONI +2.1)
3:28	CO₂ passes 350 ppm (1987-11) - the drone is louder and more open
3:34	Category 5 (ATL): Gilbert (Sep 1988, 160 kt) - tam-tam
3:36	Strong La Niña (1988-12, ONI -1.9)
3:37	Strongest winter NAO (1989-02, +2.0) - the loudest wind
3:39	Solar maximum (1989-06) - the shimmer is fullest
3:41	Category 5 (ATL): Hugo (Sep 1989, 140 kt) - tam-tam
3:54	Eruption: Pinatubo, Philippines (Jun 1991) - gong; the drone darkens as the aerosols cool the planet
3:58	Strong El Niño (1992-01, ONI +1.5)
4:02	Category 5 (ATL): Andrew (Aug 1992, 150 kt) - tam-tam
4:33	AMO phase turns warm (1997) - Dorian
4:35	Global SST anomaly passes +0.5 C (1997-06) - cello in three voices
4:39	Strong El Niño (1997-12, ONI +2.4)
4:45	Category 5 (ATL): Mitch (Oct 1998, 155 kt) - tam-tam
4:53	Strong La Niña (1999-12, ONI -1.5)
5:09	Solar maximum (2002-03) - the shimmer is fullest
5:19	Category 5 (ATL): Isabel (Sep 2003, 145 kt) - tam-tam
5:26	Category 5 (ATL): Ivan (Sep 2004, 145 kt) - tam-tam
5:32	Most active ATL season by ACE: 2005 (244)
5:32	Category 5 (ATL): Emily (Jul 2005, 140 kt) - tam-tam
5:33	Category 5 (ATL): Katrina (Aug 2005, 150 kt) - tam-tam
5:34	Category 5 (ATL): Rita (Sep 2005, 155 kt) - tam-tam
5:34	Category 5 (ATL): Wilma (Oct 2005, 160 kt) - tam-tam
5:47	Category 5 (ATL): Dean (Aug 2007, 150 kt) - tam-tam
5:47	Category 5 (ATL): Felix (Sep 2007, 150 kt) - tam-tam
5:50	Strong La Niña (2008-01, ONI -1.8)
5:52	Ocean heat content (0-2000 m) past the halfway mark of its rise (2008-06) - the contrabass adds a fifth
6:03	Strong El Niño (2009-12, ONI +1.5)
6:09	Strong La Niña (2010-11, ONI -1.6)
6:22	Lowest September sea ice in the record: 2012 (3.57 million km2) - the glass voice at its thinnest
6:33	Solar maximum (2014-04) - the shimmer is fullest
6:37	Global temperature passes +0.8 C vs 1951-80 (2014-10) - a high violin enters
6:40	CO₂ passes 400 ppm (2015-03) - the drone is louder and more open
6:45	Strong El Niño (2015-12, ONI +2.6)
6:51	Category 5 (ATL): Matthew (Oct 2016, 145 kt) - tam-tam
6:57	Category 5 (ATL): Irma (Sep 2017, 155 kt) - tam-tam
6:58	Category 5 (ATL): Maria (Sep 2017, 150 kt) - tam-tam
7:05	Category 5 (ATL): Michael (Oct 2018, 140 kt) - tam-tam
7:11	Category 5 (ATL): Dorian (Sep 2019, 160 kt) - tam-tam
7:12	Category 5 (ATL): Lorenzo (Sep 2019, 140 kt) - tam-tam
7:28	Eruption: Hunga Tonga-Hunga Ha’apai (Jan 2022) - gong; mostly water vapour, not sulphate: the drone barely darkens
7:33	Category 5 (ATL): Ian (Sep 2022, 140 kt) - tam-tam
7:35	The AMO record ends (2023-01; the index is no longer updated) - its pad falls silent and the mode holds
7:36	CO₂ passes 420 ppm (2023-04) - the drone is louder and more open
7:39	Category 5 (ATL): Lee (Sep 2023, 145 kt) - tam-tam
7:41	Strong El Niño (2023-12, ONI +2.0)
7:45	Category 5 (ATL): Beryl (Jul 2024, 145 kt) - tam-tam
7:47	Category 5 (ATL): Milton (Oct 2024, 155 kt) - tam-tam
7:53	Category 5 (ATL): Erin (Aug 2025, 140 kt) - tam-tam
7:54	Category 5 (ATL): Humberto (Sep 2025, 140 kt) - tam-tam
7:54	Category 5 (ATL): Melissa (Oct 2025, 165 kt) - tam-tam

What to be careful about

• Early storms are undercounted. Before routine satellite coverage (about 1966), and in the Northeast Pacific into the 1970s, weaker storms were missed. The sparse early drums are partly fewer observed storms.
• Early West Pacific intensities are uncertain. The US Joint Typhoon Warning Center’s wind estimates from the 1950s–70s are thought to run high, so some early “Category 5” tam-tams may be overstated. The Pacific guide names only the eight strongest West Pacific Category 5 storms.
• Records end at different times. The AMO index stopped updating in January 2023, so its pad falls silent there. West Pacific best tracks currently end in December 2024. Both gaps are marked in the guides, not filled.
• Sea ice starts in November 1978, when the satellite record starts. Longer histories reconstructed from ship logs, ice charts and aircraft exist, but they are a different kind of record and are not used here.
• The ocean phases lag. The PDO and AMO are smoothed over three years, so their turns arrive about two years after the conventional dates (the 1976–77 Pacific shift sounds in 1980).
• Interpretive choices include the volcanic cooling shape (Hunga Tonga is treated as mostly water vapour and barely darkens the drone), the ENSO-to-chord table, and every loudness curve. None of them is a measurement.

Sources and licences

• CO₂: NOAA Global Monitoring Laboratory, Mauna Loa monthly mean (the Keeling curve).
• ENSO, NAO: NOAA Climate Prediction Center, Oceanic Niño Index and North Atlantic Oscillation.
• PDO, global ocean surface temperature: NOAA National Centers for Environmental Information (ERSSTv5; Climate at a Glance).
• AMO: NOAA Physical Sciences Laboratory.
• Global temperature: NASA GISS Surface Temperature Analysis (GISTEMP v4).
• Ocean heat content: Cheng, L. et al. (2024). IAPv4 ocean temperature and ocean heat content gridded dataset. Earth System Science Data 16, 3517–3546. doi:10.5194/essd-16-3517-2024. Institute of Atmospheric Physics, Chinese Academy of Sciences.
• Storms: NOAA National Hurricane Center HURDAT2 (Atlantic; Northeast and Central Pacific). NOAA NCEI International Best Track Archive for Climate Stewardship (IBTrACS) v04r01 (West Pacific).
• Arctic sea ice: National Snow and Ice Data Center, Sea Ice Index, version 4 (G02135).
• Sunspots: WDC-SILSO, Royal Observatory of Belgium, Brussels (CC BY-NC 4.0).
• Kettle River: U.S. Geological Survey, gauge 12404500, Kettle River near Laurier, WA.
• Volcanic eruptions: Smithsonian Institution Global Volcanism Program (dates).
• Instruments: VSCO 2 Community Edition, Versilian Studios (CC0), plus synthesizers written in Python.

The sunspot numbers are licensed for non-commercial use only, and the heat-content series’ redistribution terms are not yet confirmed. Those two series are therefore not stored in the repository; fetch_climate.py downloads them. Because the music is derived in part from non-commercial data, these two pieces are released under CC BY-NC-SA 4.0, not the CC BY-SA 4.0 of Water Year Rings.

Source code

Noteclimate_pair.py: the score for both pieces

"""
CLIMATE PAIR (2026-09-25): two pieces on one multidecadal grammar, 1958-2025.

    python climate_pair.py pacific
    python climate_pair.py atlantic

GW's brief: atmospheric CO2 (Keeling), SST, ENSO, PDO and storm cycles on a
multidecadal scale; the rising elements in the drone and cello range with
increasing loudness; cyclones and hurricanes in the deeper percussion; moody
without being overly dark. Both pieces share the global layers and differ in
their ocean: the Pacific piece (PDO, West Pacific typhoons and NE+Central
Pacific hurricanes, and the Kettle River's annual peak, a valley that PDO/ENSO
steer) and the Atlantic piece (AMO, NAO, Atlantic hurricanes, Arctic sea ice). ENSO runs through both: El Nino years
bring more Pacific storms and fewer Atlantic ones, and the pair lets you hear it.

TIME: 1 bar = 1 calendar year (68 bars). The bar is felt as four dotted-
quarter beats = the four seasons; the twelve eighth notes are the months.

LAYER          DATA (R retrieved, D derived, S simulated)            SOUND
drone          Mauna Loa CO2, deseasonalized (R)                      D/A drone + organ pedal; level and filter rise
                                                                      with ppm; the seasonal CO2 cycle (R) makes it breathe
cello          global ocean SST anomaly (R)                           cello section on chord tones: 1 -> 3 voices, louder
contrabass     ocean heat content 0-2000 m, IAP v4.2 monthly (R)      solo contrabass under the cellos: the deep ocean
                                                                      filling; root, then root + 5th, louder
chords         ONI (R): El Nino lifts to VI/III, La Nina sinks to iv   half-year chords; neutral years circle i VII i v
month pulse    ONI (R)                                                one note per month; louder in strong ENSO, higher
                                                                      in El Nino, lower in La Nina
mode + pad     PDO (Pacific) / AMO (Atlantic), 3-yr mean (R/D)         Dorian in the warm phase, Aeolian in the cool;
                                                                      synth + violin pad brightens with the index
storms         HURDAT2 / IBTrACS best tracks (R); Pacific: West       tropical storm: bass-drum rub; hurricane: tuned
               Pacific panned left (west), NE Pacific right (east)
                                                                      timpani; major: timpani + concert bass drum + sub
                                                                      boom; Cat 5: tam-tam; landfall: cymbal. Season ACE
                                                                      (D) -> low storm rumble
sunspots       SILSO monthly sunspot number (R)                       high shimmer + glock sparkles (11-year cycle)
heat           GISTEMP global anomaly (R)                             a held high violin once it passes +0.8 C
volcanoes      large eruptions (R dates) + cooling decay (D)          gong; the drone darkens for ~2 years
Pacific only:  Kettle River near Laurier annual peak flow, USGS (R)    spring harp/piano bloom, sized by rank
Atlantic only: NAO (R)                                                winter wind (band noise), stronger when positive
               Arctic sea-ice extent, NSIDC, Nov 1978 on (R)          a high glass voice: breathes with the season,
                                                                      thins (fewer voices) as September ice declines
"""
import json
import subprocess
import sys
from datetime import datetime
from pathlib import Path

import numpy as np

sys.path.insert(0, str(Path(__file__).resolve().parent))
import sampler as S
from sampler import Instrument, OneShots, Bus, make_ir, SR
import seasonal as W
import synth as Y
import records as R

BASIN = (sys.argv[1] if len(sys.argv) > 1 else "pacific").lower()
if BASIN not in ("pacific", "atlantic"):
    sys.exit("usage: python climate_pair.py pacific|atlantic")
CFG = {
    # storms: basin -> (pan centre, gain); west on the left, as on a map
    "pacific": dict(name="climate_pacific_1958_2025", tonic=50, key="D", mode_index="pdo",
                    storms={"wpac": (-0.5, 0.3), "nepac": (0.45, 0.75)}, seed=1958, pulse="piano"),
    "atlantic": dict(name="climate_atlantic_1958_2025", tonic=45, key="A", mode_index="amo",
                     storms={"atl": (0.0, 1.0)}, seed=1492, pulse="harp"),
}[BASIN]
NAME, TONIC = CFG["name"], CFG["tonic"]
OUT = W.LIGHT_OUT / "climate"
HEAVY = W.HEAVY_OUT / "climate"
OUT.mkdir(parents=True, exist_ok=True)
HEAVY.mkdir(parents=True, exist_ok=True)
rng = np.random.default_rng(CFG["seed"])
np.random.seed(CFG["seed"])
NOISE_SEED = 99  # the high noise band's own generator
print(f"seeds: {CFG['seed']}, {NOISE_SEED}")  # renders.py records these in the lock

J = R.climate_monthly()  # data/climate/parquet (or the parsed JSON it was built from)
J["volcanoes"] = R.volcanoes()
MONTHS = J["months"]
NM = len(MONTHS)
Y0 = int(MONTHS[0][:4])
NY = NM // 12
ser = {k: np.array([np.nan if v is None else v for v in a], dtype=float) for k, a in J["series"].items()}
# sunspots (CC BY-NC) and ocean heat (terms unconfirmed) live in their own, uncommitted file
FO = R.climate_monthly("fetched_only_monthly")
ser.update({k: np.array([np.nan if v is None else v for v in a], dtype=float) for k, a in FO["series"].items()})


def fill(x):
    """Linear interpolation over interior gaps; edges hold the nearest value."""
    x = x.copy()
    ok = np.isfinite(x)
    x[~ok] = np.interp(np.flatnonzero(~ok), np.flatnonzero(ok), x[ok])
    return x


def smooth(x, w):
    return W.smooth(x, w)


def n01(x, lo=None, hi=None):
    lo = np.nanpercentile(x, 2) if lo is None else lo
    hi = np.nanpercentile(x, 98) if hi is None else hi
    return np.clip((x - lo) / (hi - lo), 0, 1)


# ------------------------------------------------------------------ grid ----
BAR = 7.0                     # seconds per year
MONTH = BAR / 12
TAIL = 14.0
DUR = NY * BAR + TAIL
n = int(DUR * SR)
print(f"{BASIN}: {NY} years, {DUR / 60:.1f} min")


def t_of(year_idx, month_frac=0.0):
    return year_idx * BAR + month_frac * MONTH


def t_date(dt):
    y = dt.year - Y0
    frac = (dt.timetuple().tm_yday - 1) / (366 if dt.year % 4 == 0 else 365) * 12
    return t_of(y, frac)


def smp(sec):
    return int(max(0, sec) * SR)


t_month = np.arange(n) / SR / MONTH           # sample -> fractional month index
at_sample = lambda monthly: np.interp(t_month, np.arange(NM) + 0.5, monthly)
events = []


def mark(sec, text):
    events.append((round(sec, 1), text))


# -------------------------------------------------------------- features ----
co2 = fill(ser["co2_ppm"])
co2d = fill(ser["co2_deseason_ppm"])
co2_cycle = co2 - co2d                                    # the biosphere's annual breath (R)
oni = fill(ser["oni"])
sst = fill(ser["sst_global_anom"])
gis = fill(ser["gistemp_anom"])
sun = fill(ser["sunspots"])
idx_raw = ser[CFG["mode_index"]]
idx_ok = np.isfinite(idx_raw)
idx_last = int(np.flatnonzero(idx_ok)[-1])
idx = fill(idx_raw)
idx3 = smooth(idx, 36)
nao = fill(ser["nao"])
ohc = fill(ser["ohc_0_2000_1e22J"])
ice_raw = ser["seaice_extent_mkm2"]
ice_first = int(np.flatnonzero(np.isfinite(ice_raw))[0])     # satellite record begins (Nov 1978)
ice = fill(ice_raw)

co2_n = n01(co2d, 315, 428)                               # fixed endpoints: 1958 -> ~2025
sst_n = n01(sst, -0.25, 1.05)

# volcanic cooling envelope (D): e-folding ~1.2 years after each eruption
volc = np.zeros(NM)
for v in J["volcanoes"]:
    dt = datetime.fromisoformat(v["date"])
    m0 = (dt.year - Y0) * 12 + dt.month - 1
    k = np.arange(NM) - m0
    volc = np.maximum(volc, np.where(k >= 0, np.exp(-k / 14.0) * v.get("cooling", 0.4), 0))

# ---------------------------------------------------------------- harmony ---
# mode per year from the 3-year mean of PDO/AMO with hysteresis: Dorian (warm) / Aeolian (cool)
modes = np.zeros(NY, dtype=int)
HYST = 0.25 * float(np.nanstd(idx3))   # hysteresis scaled to the index (PDO and AMO differ ~4x in range)
cur = 2 if idx3[5] > 0 else 1
for y in range(NY):
    v = idx3[y * 12 + 6]
    if v > HYST:
        cur = 2
    elif v < -HYST:
        cur = 1
    modes[y] = cur
NEUTRAL = [0, 6, 0, 4]            # i VII i v
chords = []                       # per half-year
for h in range(NY * 2):
    o = float(np.mean(oni[h * 6:(h + 1) * 6]))
    if o >= 1.5:
        d = 2                     # III: a strong El Nino lifts furthest
    elif o >= 0.5:
        d = 5                     # VI
    elif o <= -0.5:
        d = 3                     # iv
    else:
        d = NEUTRAL[h % 4]
    chords.append(d)


def tones(h, lo, hi, prev=None, ext=(0, 2, 4, 8)):
    return W.voice_lead(prev, W.chord(chords[h], modes[h // 2], TONIC, ext), lo, hi)


def root(h, lo, hi):
    return W.voice_lead(None, [W.degree_to_midi(chords[h], modes[h // 2], TONIC)], lo, hi)[0]


# ------------------------------------------------------------ instruments ---
K = S.vsco_kit()
piano, harp, cello, vln, glock, cbass = K["piano"], K["harp"], K["cello"], K["vln"], K["glock"], K["bass"]
perc = OneShots()
timp = OneShots(S.VSCO / "Percussion" / "Timpani", glob="*Hit_*.wav")
TIMP_F0 = {1: 41.5, 2: 46.8, 3: 49.5, 4: 52.4, 5: 54.6}     # principal-mode MIDI, FFT-measured 2026-09-25
organ = Instrument("Keys/Organ/Quiet", glob="NT5_PedalQuiet_*.wav",
                   mapping=lambda nm: (int(nm.split("_")[2]) - 28, 1))
pulse_ins = piano if CFG["pulse"] == "piano" else harp

B = {k: Bus(k, n, **kw) for k, kw in {
    "drone":  dict(gain=1.8, hp=25, send=0.3),
    "cello":  dict(gain=1.1, hp=45, send=0.45),
    "deep":   dict(gain=1.8, hp=28, lp=900, send=0.3),       # ocean heat content: contrabass
    "pad":    dict(gain=0.7, hp=140, send=0.6),
    "pulse":  dict(gain=3.6 if CFG["pulse"] == "piano" else 1.7, hp=90, send=0.45),
    "storms": dict(gain=0.8, hp=25, send=0.35),
    "rumble": dict(gain=1.0, hp=30, lp=420, send=0.2),
    "shimmer": dict(gain=0.5, hp=1500, send=0.6),
    "heat":   dict(gain=2.0, hp=500, send=0.5),
    "events": dict(gain=1.3, hp=40, send=0.5),
    "local":  dict(gain=1.5, hp=80, send=0.5),     # Kettle bloom / NAO wind
}.items()}
if BASIN == "atlantic":
    B["ice"] = Bus("ice", n, gain=0.5, hp=1200, send=0.65)   # Arctic sea ice


def repitch(x, semis):
    r = 2 ** (semis / 12)
    idx_ = np.arange(int(len(x) / r)) * r
    return np.stack([np.interp(idx_, np.arange(len(x)), x[:, c]) for c in range(2)], axis=1)


def tuned_timp(target_midi, dyn):
    k = min(TIMP_F0, key=lambda i: abs(((target_midi - TIMP_F0[i] + 6) % 12) - 6))
    shift = ((target_midi - TIMP_F0[k] + 6) % 12) - 6
    try:
        x = timp.hit(f"Timpani{k}_Hit_v{dyn}", 1.0, rng)
    except FileNotFoundError:          # Timpani1 has only v1 and v3
        x = timp.hit(f"Timpani{k}_Hit_v3", 1.0, rng)
    return repitch(x, shift)


# ------------------------------------------------------------ drone (CO2) ---
lvl_m = 0.18 + 0.82 * co2_n ** 1.3                          # rising loudness
cut_m = (160 + 1500 * co2_n ** 1.5) * (1 - volc)            # opens with ppm, darkens after eruptions
breathe_m = n01(co2_cycle, -4, 4)                           # seasonal CO2 cycle, 0..1
B["drone"].add(Y.drone([TONIC - 12, TONIC - 5], n, at_sample(cut_m), at_sample(lvl_m), rng=rng,
                       breathe=at_sample(breathe_m)), 0)
for y in range(NY):                                          # sampled organ pedal joins as CO2 rises
    c = float(np.mean(co2_n[y * 12:(y + 1) * 12]))
    if c > 0.2:
        B["drone"].add(W.sustain(organ, TONIC - 12, 0.5, BAR, rng), smp(t_of(y)), gain=0.25 + 0.9 * (c - 0.2))
for ppm in (350, 400, 420):
    m = int(np.argmax(co2d >= ppm))
    if co2d[m] >= ppm:
        mark(t_of(m // 12, m % 12), f"CO2 passes {ppm} ppm ({MONTHS[m]}) - the drone is louder and more open")

# ------------------------------------------------------------ cello (SST) ---
prev = None
for h in range(NY * 2):
    s_ = float(np.mean(sst_n[h * 6:(h + 1) * 6]))
    nv = 1 + int(s_ > 0.35) + int(s_ > 0.7)
    full = tones(h, 38, 62, prev, ext=(0, 4, 8))
    prev = full
    voices = [full[0]] + ([full[1]] if nv > 1 else []) + ([full[2]] if nv > 2 else [])
    for i, m in enumerate(voices):
        B["cello"].add(W.sustain(cello, m, 0.3 + 0.55 * s_, BAR / 2, rng), smp(t_of(h / 2)),
                       gain=(0.3 + 0.9 * s_ ** 1.2) * (1.0 if i == 0 else 0.75), pan=(-0.2, 0.25, -0.05)[i])
m = int(np.argmax(sst >= 0.5))
mark(t_of(m // 12, m % 12), f"Global SST anomaly passes +0.5 C ({MONTHS[m]}) - cello in three voices")

# ------------------------------------------- contrabass (ocean heat 0-2000 m) --
# GW's note: a second, lower voice for the deep ocean, below the SST cellos.
ohc_n = n01(ohc, float(np.nanmin(ohc)), float(np.nanmax(ohc)))
for h in range(NY * 2):
    o_ = float(np.mean(ohc_n[h * 6:(h + 1) * 6]))
    r_ = root(h, 28, 40)
    notes = [r_] + ([r_ + 7] if o_ > 0.55 else [])
    for i, m in enumerate(notes):
        B["deep"].add(W.sustain(cbass, m, 0.3 + 0.5 * o_, BAR / 2, rng), smp(t_of(h / 2)),
                      gain=(0.15 + 1.0 * o_ ** 1.4) * (1.0 if i == 0 else 0.6), pan=(-0.1, 0.2)[i])
m = int(np.argmax(ohc_n >= 0.55))
mark(t_of(m // 12, m % 12), f"Ocean heat content (0-2000 m) past the halfway mark of its rise ({MONTHS[m]}) - "
                            "the contrabass adds a fifth")

# -------------------------------------------------------- pad (PDO / AMO) ---
prev = None
IDXN = n01(idx3, -1.5 if CFG["mode_index"] == "pdo" else -0.3, 1.5 if CFG["mode_index"] == "pdo" else 0.3)
for h in range(NY * 2):
    mi = h * 6 + 3
    if mi > idx_last + 3:                                    # index no longer updated: the pad falls silent
        break
    v = tones(h, 57, 81, prev, ext=(0, 2, 4, 8))
    prev = v
    b_ = float(IDXN[mi])
    for i, m in enumerate(v):
        B["pad"].add(Y.pad_note(m, BAR / 2, vel=0.35 + 0.2 * b_, cutoff=450 + 2200 * b_, a=1.6, r=2.4, rng=rng),
                     smp(t_of(h / 2)), pan=(-0.5, -0.15, 0.15, 0.5)[i % 4])
    if b_ > 0.45:
        for i, m in enumerate(v[1:]):
            B["pad"].add(W.sustain(vln, m + 12, 0.25 + 0.3 * b_, BAR / 2, rng), smp(t_of(h / 2)),
                         gain=0.35 * b_, pan=(-0.55, 0.0, 0.55)[i % 3])
if idx_last < NM - 1:
    mark(t_of(idx_last // 12, idx_last % 12),
         f"The {CFG['mode_index'].upper()} record ends ({MONTHS[idx_last]}; the index is no longer updated) - "
         "its pad falls silent and the mode holds")
for y in range(1, NY):
    if modes[y] != modes[y - 1]:
        mark(t_of(y), f"{CFG['mode_index'].upper()} phase turns {'warm' if modes[y] == 2 else 'cool'} "
                      f"({Y0 + y}) - {W.MODES[modes[y]][0]}")

# ------------------------------------------------------ month pulse (ENSO) --
PAT = [0, 2, 1, 3, 2, 1, 0, 2, 1, 3, 2, 1]
for mth in range(NM):
    h = mth // 6
    o = oni[mth]
    shift = 12 if o > 1.0 else (-12 if o < -1.0 else 0)
    pool = sorted(tones(h, 60 + shift, 79 + shift))
    note = pool[PAT[mth % 12] % len(pool)]
    vel = float(np.clip(0.22 + 0.28 * min(abs(o), 2.2) / 2.2 + (0.08 if mth % 3 == 0 else 0), 0.1, 0.8))
    at = smp(t_of(mth // 12, mth % 12) + rng.normal(0, 0.006))
    B["pulse"].add(pulse_ins.note(note, vel, dur=MONTH * 1.6, release=1.2, rng=rng), at,
                   pan=0.3 * np.sin(mth * 0.52))
    if abs(o) > 1.0:  # strong events: a quiet synth pluck doubles the pulse an octave up
        B["pulse"].add(Y.pluck_note(note + 12, vel=0.25 * vel, decay=0.35, cutoff=1800, rng=rng), at, pan=-0.3)
for thr, label in ((1.5, "El Nino"), (-1.5, "La Nina")):
    y_prev = -9
    for mth in range(NM):
        if (oni[mth] >= thr if thr > 0 else oni[mth] <= thr) and mth // 12 > y_prev + 1:
            seg = oni[mth:mth + 18]
            pk = mth + int(np.argmax(seg) if thr > 0 else np.argmin(seg))
            mark(t_of(pk // 12, pk % 12), f"Strong {label} ({MONTHS[pk]}, ONI {oni[pk]:+.1f})")
            y_prev = pk // 12

# ---------------------------------------------------------------- storms ----
rumble_env = np.zeros(NM * 8)                                # 1/8-month resolution
duck_hits, duck_depth = [], []
ace_year = np.zeros(NY)
cat_counts = {}
for basin, (pan0, bgain) in CFG["storms"].items():
    storms = R.storms(basin)
    cc = np.zeros(6, dtype=int)
    cat5 = []
    for s_ in storms:
        dt = datetime.fromisoformat(s_["peak_time"])
        if dt.year < Y0 or dt.year >= Y0 + NY:
            continue
        t = t_date(dt)
        h = min(int(t / (BAR / 2)), NY * 2 - 1)
        cat = s_["category"]
        cc[cat] += 1
        ace_year[dt.year - Y0] += s_["ace"]
        r_ = root(h, 38, 50)
        at = smp(t + rng.normal(0, 0.01))
        pan = float(np.clip(pan0 + rng.normal(0, 0.25), -0.85, 0.85))
        if cat == 0:
            B["storms"].add(perc.hit("bassdrum_rub", 0.35, rng), at, gain=0.5 * bgain, pan=pan)
        else:
            dyn = 1 if cat <= 1 else (3 if cat == 2 else 4)
            B["storms"].add(tuned_timp(r_ if cat != 2 else r_ + 7, dyn), at, gain=(0.5 + 0.12 * cat) * bgain, pan=pan)
        if cat >= 3:
            B["storms"].add(perc.hit(f"BDrumNewhit_v{min(7, 3 + cat)}", 1.0, rng), at, gain=(0.6 + 0.1 * cat) * bgain,
                            pan=pan * 0.5)
            B["storms"].add(Y.kick(0.6 + 0.08 * cat, f0=38, f_start=90, decay=0.9 + 0.15 * cat, click=0.05), at,
                            gain=0.7 * bgain)
            duck_hits.append(at)
        if cat == 5:
            B["storms"].add(perc.hit("gongHit_f", 1.0, rng), at + int(0.05 * SR), gain=0.55 * bgain, pan=pan * 0.6)
            cat5.append((s_["peak_kt"], t, s_["name"], dt))
        if s_["landfall"] and cat >= 1:
            B["storms"].add(perc.hit("susCymb1-hit_mp", 0.6, rng), at + int(0.03 * SR), gain=(0.25 + 0.05 * cat) * bgain,
                            pan=-pan if len(CFG["storms"]) == 1 else pan)
        k0 = int(t / MONTH * 8)
        w = np.arange(-8, 9)
        lo, hi = max(0, k0 - 8), min(len(rumble_env), k0 + 9)
        rumble_env[lo:hi] += np.sqrt(max(s_["ace"], 0.1)) * np.exp(-0.5 * (w[lo - k0 + 8:hi - k0 + 8] / 4.0) ** 2)
    # the guide names every Cat 5 in a basin with few of them, and the eight strongest where there are many
    for kt, t, nm, dt in (cat5 if len(cat5) <= 40 else sorted(cat5, reverse=True)[:8]):
        mark(t, f"Category 5 ({basin.upper()}): {nm} ({dt:%b %Y}, {kt} kt) - tam-tam")
    last = max(datetime.fromisoformat(s_["peak_time"]) for s_ in storms)
    if last.year < Y0 + NY - 1:
        mark(t_of(last.year + 1 - Y0), f"{basin.upper()} best tracks end {last:%b %Y} (not yet in the archive) - "
                                       "that basin's drums fall silent")
    cat_counts[basin] = cc.tolist()
    print(f"storms {basin} by category", cc.tolist(), f"({len(cat5)} Cat 5)")
rumble_s = np.interp(t_month * 8, np.arange(len(rumble_env)), rumble_env)
rumble_s = rumble_s / (np.percentile(rumble_env[rumble_env > 0], 98) + 1e-9)
B["rumble"].add(Y.noise_band(n, 35, 260, rng=rng) * np.clip(rumble_s, 0, 1.4)[:, None] * 0.35, 0)
top = int(np.argmax(ace_year))
basins = "+".join(b.upper() for b in CFG["storms"])
mark(t_of(top, 6), f"Most active {basins} season by ACE: {Y0 + top} ({ace_year[top]:.0f})")
mark(t_of(1966 - Y0), "1966: routine satellite coverage begins - fewer storms before this is partly fewer observed")

# ------------------------------------------------------ shimmer (sunspots) --
sun_n = n01(smooth(sun, 6), 0, 250)
for h in range(NY * 2):
    s_ = float(np.mean(sun_n[h * 6:(h + 1) * 6]))
    if s_ < 0.12:
        continue
    for m in tones(h, 84, 98, ext=(0, 4, 8))[:2]:
        B["shimmer"].add(Y.pad_note(m, BAR / 2, vel=0.2 + 0.5 * s_, cutoff=6000, a=2.0, r=2.5, rng=rng, voices=3,
                                    detune=6), smp(t_of(h / 2)), gain=s_)
    for _ in range(int(round(6 * s_))):
        m = int(rng.choice(tones(h, 86, 100)))
        B["shimmer"].add(glock.note(m, 0.18, dur=0.8, release=2.0, rng=rng),
                         smp(t_of(h / 2) + rng.uniform(0, BAR / 2)), gain=0.5, pan=rng.uniform(-0.8, 0.8))
for mth in range(24, NM - 24):
    s6 = smooth(sun, 13)
    if s6[mth] == s6[mth - 24:mth + 24].max() and s6[mth] > 80:
        mark(t_of(mth // 12, mth % 12), f"Solar maximum ({MONTHS[mth]}) - the shimmer is fullest")

# ------------------------------------------------------- heat (GISTEMP) -----
gis12 = smooth(gis, 12)
for h in range(NY * 2):
    a = float(np.mean(gis12[h * 6:(h + 1) * 6]))
    if a < 0.8:
        continue
    m = W.voice_lead(None, [W.degree_to_midi(8, modes[h // 2], TONIC)], 81, 93)[0]   # the 9th, held high
    B["heat"].add(W.sustain(vln, m, 0.3, BAR / 2, rng), smp(t_of(h / 2)), gain=min(1.0, (a - 0.8) / 0.5))
m = int(np.argmax(gis12 >= 0.8))
mark(t_of(m // 12, m % 12), f"Global temperature passes +0.8 C vs 1951-80 ({MONTHS[m]}) - a high violin enters")

# -------------------------------------------------------------- volcanoes ---
for v in J["volcanoes"]:
    dt = datetime.fromisoformat(v["date"])
    t = t_date(dt)
    B["events"].add(perc.hit("gongscrape_mf", 0.9, rng), smp(t - 1.2), gain=0.6)
    B["events"].add(perc.hit("gongHit_mf", 0.9, rng), smp(t), gain=0.8)
    mark(t, f"Eruption: {v['name']} ({dt:%b %Y}) - gong; " + v.get("note", "the drone darkens as the aerosols cool the planet"))

# ------------------------------------------------------ local layer ---------
if BASIN == "pacific":
    # GW's note: the Kettle at Laurier (USGS, complete since 1929) replaces the Granby, whose record has a 1958-66 gap
    gr = R.kettle_annual_peak()
    vals = np.array([v["max_daily_m3s"] for k_, v in gr.items() if Y0 <= int(k_) < Y0 + NY])
    for y in range(NY):
        rec = gr.get(str(Y0 + y))
        if rec is None:
            continue
        dt = datetime.fromisoformat(rec["date"])
        t = t_date(dt)
        p = float((vals < rec["max_daily_m3s"]).mean())
        h = min(int(t / (BAR / 2)), NY * 2 - 1)
        v_ = sorted(tones(h, 55, 86))
        k = 3 + int(round(5 * p))
        for i, m in enumerate((v_ + [x + 12 for x in v_])[:k]):
            B["local"].add(harp.note(m, 0.3 + 0.35 * p, rng=rng), smp(t + i * 0.07), pan=-0.5 + i * 0.12)
        B["local"].add(piano.note(v_[0], 0.3 + 0.3 * p, dur=2.0, release=1.5, rng=rng), smp(t), pan=0.1)
    missing = [Y0 + y for y in range(NY) if str(Y0 + y) not in gr and Y0 + y <= 2024]
    top_y = max(gr, key=lambda k_: gr[k_]["max_daily_m3s"] if Y0 <= int(k_) < Y0 + NY else -1)
    mark(t_date(datetime.fromisoformat(gr[top_y]["date"])),
         f"The Kettle's biggest peak in the span: {gr[top_y]['max_daily_m3s']:.0f} m3/s ({gr[top_y]['date']}) - "
         "the fullest harp bloom")
    if missing:
        mark(t_of(missing[0] - Y0), f"No Kettle record {missing[0]}-{missing[-1]} - the spring bloom is silent")
else:
    wind = np.zeros(NM)
    for mth in range(NM):
        if (mth % 12) in (11, 0, 1, 2):                      # Dec-Mar
            wind[mth] = np.clip(0.25 + 0.3 * nao[mth], 0, 1.2)
    w_s = at_sample(smooth(wind, 2))
    lo_n = Y.noise_band(n, 300, 1400, rng=rng)
    hi_n = Y.noise_band(n, 1400, 4200, rng=np.random.default_rng(NOISE_SEED))
    gust = 0.7 + 0.3 * np.sin(2 * np.pi * np.arange(n) / SR / 5.3) * np.sin(2 * np.pi * np.arange(n) / SR / 2.1)
    B["local"].add((lo_n * 0.6 + hi_n * 0.25) * (w_s * gust)[:, None] * 0.22, 0)
    top = int(np.argmax(np.where((np.arange(NM) % 12) < 3, nao, -9)))
    mark(t_of(top // 12, top % 12), f"Strongest winter NAO ({MONTHS[top]}, {nao[top]:+.1f}) - the loudest wind")

    # Arctic sea ice (GW's note: a careful bring-in when the satellite record starts). A glass voice of
    # sine partials on high chord tones: its level follows the monthly extent (the seasonal breath,
    # March max, September min); its number of voices follows that year's September extent, so it
    # thins over the decades. Silent before Nov 1978: no satellite record, and nothing is invented.
    sep = {y: ice[y * 12 + 8] for y in range(NY) if y * 12 + 8 >= ice_first}
    sep_lo, sep_hi = min(sep.values()), max(sep.values())
    lvl_ice = np.where(np.arange(NM) >= ice_first, n01(ice, 3.0, 16.5), 0.0)
    ramp = np.clip((np.arange(NM) - ice_first) / 3.0, 0, 1)            # a three-month fade-in, not a jump
    lvl_ice_s = at_sample(lvl_ice * ramp)
    for h in range(NY * 2):
        m0 = h * 6
        if m0 + 6 <= ice_first:
            continue
        y = h // 2
        sv = sep.get(y, sep.get(y - 1, sep_hi))
        nv = 1 + int(round(3 * (sv - sep_lo) / (sep_hi - sep_lo)))      # 4 voices at the most ice, 1 at the least
        pool = sorted(tones(h, 84, 100, ext=(0, 2, 4, 8)))[:nv]
        L = int((BAR / 2 + 2.5) * SR)
        tt = np.arange(L) / SR
        env = np.minimum(1, tt / 1.5) * np.clip((BAR / 2 + 2.5 - tt) / 2.5, 0, 1)
        for i, m in enumerate(pool):
            f = Y.midi_hz(m)
            trem = 1 + 0.25 * np.sin(2 * np.pi * (0.23 + 0.07 * i) * tt + rng.uniform(0, 6.28))
            x = (np.sin(2 * np.pi * f * tt) + 0.18 * np.sin(2 * np.pi * 2.005 * f * tt)) * env * trem * 0.08
            a0 = smp(t_of(h / 2))
            g_ = lvl_ice_s[a0:a0 + L]
            x = x[:len(g_)] * g_
            B["ice"].add(x[:, None].repeat(2, 1), a0, pan=(-0.6, 0.5, -0.2, 0.25)[i % 4])
    mark(t_of(ice_first // 12, ice_first % 12), f"Satellites begin measuring Arctic sea ice ({MONTHS[ice_first]}) - "
                                                  "a high glass voice fades in")
    y_min = min(sep, key=sep.get)
    mark(t_of(y_min, 8), f"Lowest September sea ice in the record: {Y0 + y_min} ({sep[y_min]:.2f} million km2) - "
                         "the glass voice at its thinnest")

# ----------------------------------------------------- sidechain + echoes ---
duck = Y.sidechain(n, duck_hits, depth=0.35, release=0.6).astype(np.float32)
for k in ("drone", "cello", "pad"):
    B[k].buf *= duck[:, None]
B["pulse"].buf += Y.pingpong(B["pulse"].buf, MONTH * 1.5, fb=0.35, mix=0.3).astype(np.float32)
B["shimmer"].buf += Y.pingpong(B["shimmer"].buf, MONTH, fb=0.5, mix=0.4).astype(np.float32)

# ---------------------------------------------------------------- master ----
irs = {"hall": make_ir(4.5, bright=0.45, seed=21)}
mix, stems = S.master(list(B.values()), irs, HEAVY / f"{NAME}.wav", stems_dir=HEAVY / f"{NAME}_stems", wet_gain=0.42)
S.to_mp3(HEAVY / f"{NAME}.wav", OUT / f"{NAME}.mp3")
AUDIO = W.SITE_AUDIO                                          # the site's copy: 160 kbps, encoded from the WAV
AUDIO.mkdir(parents=True, exist_ok=True)
subprocess.run(["ffmpeg", "-loglevel", "error", "-y", "-i", str(HEAVY / f"{NAME}.wav"), "-af",
                "loudnorm=I=-16:TP=-1.5:LRA=11", "-ar", "44100", "-b:a", "160k", str(AUDIO / f"{NAME}.mp3")], check=True)


def active_db(x, thr_db=-50):
    w = int(0.1 * SR)
    m_ = (x ** 2).mean(axis=1)[: len(x) // w * w].reshape(-1, w).mean(axis=1)
    db = 10 * np.log10(m_ + 1e-12)
    on = db > thr_db
    return float(np.median(db[on])) if on.any() else None


def seg_db(x, a, b):
    s_ = x[smp(a):smp(b)]
    return float(10 * np.log10((s_ ** 2).mean() + 1e-12))


events.sort()
level = {k: active_db(v) for k, v in stems.items()}
first10, last10 = seg_db(mix, 5, 75), seg_db(mix, NY * BAR - 70, NY * BAR)
SCORE = {"name": NAME, "basin": BASIN, "key": CFG["key"], "years": [Y0, Y0 + NY - 1], "bar_seconds": BAR,
           "duration_s": DUR, "events": events, "stem_active_db": level,
           "storms_by_category": cat_counts,
           "mix_db_first_decade": first10, "mix_db_last_decade": last10,
           "years_detail": [{"year": Y0 + y, "mode": W.MODES[modes[y]][0], "chords": chords[2 * y:2 * y + 2],
                             "co2_ppm": round(float(np.mean(co2d[y * 12:(y + 1) * 12])), 1),
                             "sst_anom": round(float(np.mean(sst[y * 12:(y + 1) * 12])), 2),
                             "ace": round(float(ace_year[y]), 1),
                             "ohc_1e22J": round(float(np.mean(ohc[y * 12:(y + 1) * 12])), 2),
                             "seaice_sep_mkm2": (round(float(ice_raw[y * 12 + 8]), 2)
                                                 if np.isfinite(ice_raw[y * 12 + 8]) else None)}
                            for y in range(NY)]}
with open(OUT / f"{NAME}_score.json", "w") as fh:
    json.dump(SCORE, fh, indent=1)
print("\n".join(f"{int(t // 60)}:{int(t % 60):02d}  {e}" for t, e in events))
print("stem active dBFS:", {k: (round(v, 1) if v is not None else None) for k, v in level.items()})
print(f"mix loudness first decade {first10:.1f} dB, last decade {last10:.1f} dB")
print("peak", float(np.abs(mix).max()), "nan", bool(np.isnan(mix).any()))

Notefetch_climate.py: data downloads and parsing

"""
Fetch and parse the records for the climate pair (2026-09-25):
monthly global climate indices 1958-2025, Atlantic and NE/Central Pacific
tropical-cyclone tracks (HURDAT2), and the Granby River's annual maximum
daily flow (the Pacific piece's link back to the valley).

Raw downloads are cached in data/climate/raw/ (skipped if present); parsed
output goes to data/climate/indices_monthly.json and storms_<basin>.json.

Sources (all public; attribution text in the pieces' credits):
  CO2      NOAA GML Mauna Loa monthly mean (Keeling curve)            co2_mm_mlo.txt
  ONI      NOAA CPC Oceanic Nino Index (3-month running Nino-3.4 anomaly, ERSSTv5)
  PDO      NOAA NCEI Pacific Decadal Oscillation (ERSSTv5)
  AMO      NOAA PSL Atlantic Multidecadal Oscillation, unsmoothed, long
  NAO      NOAA CPC North Atlantic Oscillation, monthly standardized
  SST      NOAA NCEI Climate at a Glance: global ocean surface temperature anomaly
  GISTEMP  NASA GISS global land-ocean temperature anomaly (v4)
  Sunspots SILSO (Royal Observatory of Belgium) monthly mean total sunspot number, v2
  HURDAT2  NOAA NHC best tracks: Atlantic and NE+Central Pacific
  IBTrACS  NOAA NCEI v04r01, West Pacific; USA_WIND (JTWC 1-minute) so categories
           match HURDAT2's Saffir-Simpson scale (added 2026-09-25, GW's note)
  Sea ice  NSIDC Sea Ice Index v4 (G02135), Arctic monthly extent, 1978-11 onward
  OHC      IAP (Institute of Atmospheric Physics, CAS) v4.2, global ocean heat content
           0-2000 m, monthly, 10^22 J (Cheng et al. 2024, ESSD 16:3517). Chosen over
           NOAA NCEI's pentadal series: monthly, 1940 on, no splice needed.
  Kettle   USGS 12404500 Kettle River near Laurier, WA: annual peak daily flow, from
           the daily record in ../kettle_laurier_12404500_lowflow.json (GW's note:
           the Kettle replaces the Granby, whose record has a 1958-66 gap)
  Granby   WSC 08NN002 annual statistics via ECCC OGC API (kept; no longer used)

Licence note: SILSO sunspot data are CC BY-NC 4.0, and the IAP heat-content terms
are not yet confirmed, so those two series go to their own file
(fetched_only_monthly.json), which is not committed to the repo.
"""
import json
import re
import urllib.request
from datetime import datetime
from pathlib import Path

import numpy as np

HERE = Path(__file__).resolve().parent.parent
RAW = HERE / "data" / "climate" / "raw"
OUT = HERE / "data" / "climate"
RAW.mkdir(parents=True, exist_ok=True)
UA = {"User-Agent": "curl/8 (sonification-workbench; personal research)"}
Y0, Y1 = 1958, 2025

SOURCES = {
    "co2_mm_mlo.txt": "https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_mm_mlo.txt",
    "oni.ascii.txt": "https://www.cpc.ncep.noaa.gov/data/indices/oni.ascii.txt",
    "ersst.v5.pdo.dat": "https://www.ncei.noaa.gov/pub/data/cmb/ersst/v5/index/ersst.v5.pdo.dat",
    "amon.us.long.data": "https://psl.noaa.gov/data/correlation/amon.us.long.data",
    "nao.monthly.ascii": "https://www.cpc.ncep.noaa.gov/products/precip/CWlink/pna/norm.nao.monthly.b5001.current.ascii",
    "global_ocean_sst.csv": "https://www.ncei.noaa.gov/access/monitoring/climate-at-a-glance/global/time-series/globe/ocean/tavg/1/0/1950-2025/data.csv",
    "gistemp_glb.csv": "https://data.giss.nasa.gov/gistemp/tabledata_v4/GLB.Ts+dSST.csv",
    "sunspots_silso.csv": "https://www.sidc.be/SILSO/INFO/snmtotcsv.php",
    "hurdat2_atl.txt": "https://www.nhc.noaa.gov/data/hurdat/hurdat2-1851-2025-091226.txt",
    "hurdat2_nepac.txt": "https://www.nhc.noaa.gov/data/hurdat/hurdat2-nepac-1949-2025-091426.txt",
    "ibtracs_wp.csv": "https://www.ncei.noaa.gov/data/international-best-track-archive-for-climate-stewardship-ibtracs/"
                      "v04r01/access/csv/ibtracs.WP.list.v04r01.csv",
    "iap_ohc_v4.2_monthly.txt": "http://www.ocean.iap.ac.cn/ftp/images_files/IAPv4.2_OHC_estimate_update.txt",
    **{f"seaice_N_{m:02d}_extent_v4.0.csv": f"https://noaadata.apps.nsidc.org/NOAA/G02135/north/monthly/data/"
                                           f"N_{m:02d}_extent_v4.0.csv" for m in range(1, 13)},
    "granby_annual_stats.json": "https://api.weather.gc.ca/collections/hydrometric-annual-statistics/items"
                                "?STATION_NUMBER=08NN002&limit=500&f=json",
}

# Large volcanic eruptions in the span, for the "cooling" event layer. Dates are
# the paroxysmal eruption dates as listed by the Smithsonian Global Volcanism
# Program [verify each against GVP before publication].
# cooling = relative weight of the drone-darkening envelope (D, interpretive):
# Pinatubo is the benchmark stratospheric-aerosol cooling of the record; Hunga
# Tonga injected mostly water vapour, with a small and debated net effect.
VOLCANOES = [
    ("1963-03-17", "Agung, Bali", 5, 0.4, None),
    ("1982-04-04", "El Chichon, Mexico", 5, 0.4, None),
    ("1991-06-15", "Pinatubo, Philippines", 6, 0.6, None),
    ("2022-01-15", "Hunga Tonga-Hunga Ha'apai", 5, 0.1,
     "mostly water vapour, not sulphate: the drone barely darkens"),
]


def get(name):
    p = RAW / name
    if not p.exists() or not p.stat().st_size:
        req = urllib.request.Request(SOURCES[name], headers=UA)
        p.write_bytes(urllib.request.urlopen(req, timeout=120).read())
        print("fetched", name)
    return p.read_text(errors="replace")


MONTHS = [(y, m) for y in range(Y0, Y1 + 1) for m in range(1, 13)]
IDX = {ym: i for i, ym in enumerate(MONTHS)}


def blank():
    return np.full(len(MONTHS), np.nan)


def put(arr, y, m, v, missing=None):
    if (y, m) in IDX and v is not None and np.isfinite(v) and (missing is None or abs(v - missing) > 1e-6):
        arr[IDX[(y, m)]] = v


def parse_co2():
    avg, des = blank(), blank()
    for ln in get("co2_mm_mlo.txt").splitlines():
        if ln.startswith("#") or not ln.strip():
            continue
        f = ln.split()
        y, m, a, d = int(f[0]), int(f[1]), float(f[3]), float(f[4])
        put(avg, y, m, a, -99.99)
        put(des, y, m, d, -99.99)
    return avg, des


def parse_oni():
    # rows: SEAS YR TOTAL ANOM; the 3-month season is assigned to its centre month
    centre = {"DJF": 1, "JFM": 2, "FMA": 3, "MAM": 4, "AMJ": 5, "MJJ": 6, "JJA": 7, "JAS": 8,
              "ASO": 9, "SON": 10, "OND": 11, "NDJ": 12}
    out = blank()
    for ln in get("oni.ascii.txt").splitlines()[1:]:
        f = ln.split()
        if len(f) == 4 and f[0] in centre:
            put(out, int(f[1]), centre[f[0]], float(f[3]))
    return out


def parse_year_by_12(text, missing):
    out = blank()
    for ln in text.splitlines():
        f = ln.split()
        if len(f) >= 13 and re.fullmatch(r"\d{4}", f[0]):
            for m in range(12):
                put(out, int(f[0]), m + 1, float(f[m + 1]), missing)
    return out


def parse_nao():
    out = blank()
    for ln in get("nao.monthly.ascii").splitlines():
        f = ln.split()
        if len(f) == 3 and f[0].isdigit():
            put(out, int(f[0]), int(f[1]), float(f[2]))
    return out


def parse_ncei_csv():
    out = blank()
    for ln in get("global_ocean_sst.csv").splitlines():
        f = ln.split(",")
        if len(f) >= 2 and re.fullmatch(r"\d{6}", f[0]):
            put(out, int(f[0][:4]), int(f[0][4:]), float(f[1]))
    return out


def parse_gistemp():
    out = blank()
    for ln in get("gistemp_glb.csv").splitlines():
        f = ln.split(",")
        if f and re.fullmatch(r"\d{4}", f[0]):
            for m in range(12):
                if f[m + 1].strip() not in ("***", ""):
                    put(out, int(f[0]), m + 1, float(f[m + 1]))
    return out


def parse_sunspots():
    out = blank()
    for ln in get("sunspots_silso.csv").splitlines():
        f = ln.split(";")
        if len(f) >= 4:
            put(out, int(f[0]), int(f[1]), float(f[3]), -1)
    return out


def saffir_simpson(kt):
    for cat, lo in ((5, 137), (4, 113), (3, 96), (2, 83), (1, 64)):
        if kt >= lo:
            return cat
    return 0


def parse_hurdat(name):
    """Per storm: name, peak wind (kt) and its date, min pressure, ACE, landfall, category."""
    storms, cur = [], None
    for ln in get(name).splitlines():
        f = [x.strip() for x in ln.split(",")]
        if re.fullmatch(r"(AL|EP|CP)\d{6}", f[0]):
            cur = {"id": f[0], "name": f[1].title(), "rows": []}
            storms.append(cur)
        elif cur is not None and re.fullmatch(r"\d{8}", f[0]):
            cur["rows"].append(f)
    out = []
    for s in storms:
        y = int(s["id"][-4:])
        if not (Y0 <= y <= Y1):
            continue
        best, ace, landfall, pmin = None, 0.0, False, None
        for f in s["rows"]:
            date, hhmm, rec, status = f[0], f[1], f[2], f[3]
            v = int(f[6]) if f[6] not in ("", "-99", "-999") else 0
            p = int(f[7]) if f[7] not in ("", "-999") else None
            if rec == "L":
                landfall = True
            if p is not None and p > 0:
                pmin = p if pmin is None else min(pmin, p)
            if hhmm in ("0000", "0600", "1200", "1800") and status in ("TS", "HU", "SS") and v >= 34:
                ace += v * v / 1e4
            if best is None or v > best[0]:
                best = (v, datetime.strptime(date + hhmm, "%Y%m%d%H%M"))
        if best is None:
            continue
        out.append({"id": s["id"], "name": s["name"], "peak_kt": best[0], "peak_time": best[1].isoformat(),
                    "category": saffir_simpson(best[0]), "ace": round(ace, 2), "landfall": landfall,
                    "min_pressure_mb": pmin})
    return out


def parse_ibtracs_wp():
    """West Pacific storms from IBTrACS, same record shape as parse_hurdat()."""
    import csv
    rows = csv.reader(get("ibtracs_wp.csv").splitlines())
    head = next(rows)
    next(rows)                                         # units row
    c = {k: i for i, k in enumerate(head)}
    storms = {}
    for r in rows:
        if r[c["BASIN"]] != "WP" or r[c["TRACK_TYPE"]] != "main":
            continue
        yr = int(r[c["SEASON"]])
        if not (Y0 <= yr <= Y1):
            continue
        s = storms.setdefault(r[c["SID"]], {"id": r[c["SID"]], "name": r[c["NAME"]].title(), "rows": []})
        w = r[c["USA_WIND"]].strip()
        s["rows"].append((r[c["ISO_TIME"]], int(w) if w else None, r[c["NATURE"]], r[c["LANDFALL"]].strip(),
                          r[c["USA_PRES"]].strip()))
    out = []
    for s in storms.values():
        winds = [(w, t) for t, w, *_ in s["rows"] if w is not None]
        if not winds:
            continue
        v, t = max(winds, key=lambda x: x[0])
        ace = sum(w * w / 1e4 for tt, w, nat, *_ in s["rows"]
                  if w is not None and w >= 34 and nat == "TS" and tt[11:16] in ("00:00", "06:00", "12:00", "18:00"))
        landfall = any(lf == "0" and w is not None and w >= 34 for _, w, _, lf, _ in s["rows"])
        pres = [int(p) for *_, p in s["rows"] if p]
        out.append({"id": s["id"], "name": s["name"], "peak_kt": v,
                    "peak_time": datetime.strptime(t, "%Y-%m-%d %H:%M:%S").isoformat(),
                    "category": saffir_simpson(v), "ace": round(ace, 2), "landfall": landfall,
                    "min_pressure_mb": min(pres) if pres else None})
    return sorted(out, key=lambda s: s["peak_time"])


def parse_seaice():
    """Arctic monthly sea-ice extent (million km2); -9999 = missing (Dec 1987-Jan 1988)."""
    out = blank()
    for m in range(1, 13):
        for ln in get(f"seaice_N_{m:02d}_extent_v4.0.csv").splitlines()[1:]:
            f = [x.strip() for x in ln.split(",")]
            if len(f) >= 5 and f[0].isdigit():
                put(out, int(f[0]), int(f[1]), float(f[4]), -9999)
    return out


def parse_ohc():
    """Global 0-2000 m ocean heat content anomaly (10^22 J, baseline 2006-15) and its
    95% half-width, monthly, from IAP v4.2 (columns 9 and 11)."""
    val, err = blank(), blank()
    for ln in get("iap_ohc_v4.2_monthly.txt").splitlines():
        f = ln.split()
        if len(f) >= 11 and f[0].isdigit():
            put(val, int(f[0]), int(f[1]), float(f[8]))
            put(err, int(f[0]), int(f[1]), float(f[10]))
    return val, err


def parse_kettle():
    """Kettle River near Laurier (USGS 12404500): each calendar year's peak daily flow and date."""
    d = json.loads((HERE / "data" / "kettle_laurier_12404500_lowflow.json").read_text())["daily"]
    out = {}
    for date, q in zip(d["date"], d["q"]):
        if q is None:
            continue
        y = int(date[:4])
        if y not in out or q > out[y]["max_daily_m3s"]:
            out[y] = {"max_daily_m3s": q, "date": date[:10]}
    return dict(sorted(out.items()))


def parse_granby():
    d = json.loads(get("granby_annual_stats.json"))
    out = {}
    for f in d["features"]:
        p = f["properties"]
        if p["DATA_TYPE_EN"] == "Discharge" and p["MAX_VALUE"] is not None:
            out[int(p["IDENTIFIER"].split(".")[1])] = {"max_daily_m3s": p["MAX_VALUE"], "date": p["MAX_DATE"],
                                                       "symbol": p["MAX_SYMBOL_EN"] or ""}
    return dict(sorted(out.items()))


def main():
    co2, co2_des = parse_co2()
    pdo = parse_year_by_12(get("ersst.v5.pdo.dat"), 99.99)
    amo = parse_year_by_12(get("amon.us.long.data"), -99.99)
    series = {
        "co2_ppm": co2, "co2_deseason_ppm": co2_des, "oni": parse_oni(), "pdo": pdo, "amo": amo,
        "nao": parse_nao(), "sst_global_anom": parse_ncei_csv(), "gistemp_anom": parse_gistemp(),
        "seaice_extent_mkm2": parse_seaice(),
    }
    ohc, ohc_err = parse_ohc()
    sunspots = parse_sunspots()
    months = [f"{y}-{m:02d}" for y, m in MONTHS]
    js = {"months": months,
          "series": {k: [None if not np.isfinite(v) else round(float(v), 4) for v in a] for k, a in series.items()},
          "volcanoes": [{"date": d, "name": n, "vei": v, "cooling": c, **({"note": t} if t else {})}
                        for d, n, v, c, t in VOLCANOES],
          "sources": {k: v for k, v in SOURCES.items()}}
    (OUT / "indices_monthly.json").write_text(json.dumps(js))
    rnd = lambda a, k: [None if not np.isfinite(v) else round(float(v), k) for v in a]
    (OUT / "fetched_only_monthly.json").write_text(json.dumps(
        {"months": months, "series": {"sunspots": rnd(sunspots, 1), "ohc_0_2000_1e22J": rnd(ohc, 2),
                                      "ohc_0_2000_err95": rnd(ohc_err, 2)},
         "terms": {"sunspots": "CC BY-NC 4.0, WDC-SILSO, Royal Observatory of Belgium, Brussels",
                   "ohc": "IAP v4.2 (Cheng et al. 2024); redistribution terms not yet confirmed"}}))
    series.update(sunspots=sunspots, ohc_0_2000_1e22J=ohc)
    for k, a in series.items():
        ok = np.isfinite(a)
        first = months[int(np.argmax(ok))] if ok.any() else None
        last = months[len(ok) - 1 - int(np.argmax(ok[::-1]))] if ok.any() else None
        print(f"{k:18s} {ok.sum():4d}/{len(a)} months  {first} .. {last}")
    for basin in ("atl", "nepac", "wpac"):
        st = parse_ibtracs_wp() if basin == "wpac" else parse_hurdat(f"hurdat2_{basin}.txt")
        (OUT / f"storms_{basin}.json").write_text(json.dumps(st, indent=0))
        cats = np.bincount([s["category"] for s in st], minlength=6)
        print(f"storms {basin}: {len(st)} systems, by category {cats.tolist()}, "
              f"max ACE {max(s['ace'] for s in st):.1f}")
    ke = parse_kettle()
    (OUT / "kettle_annual_max.json").write_text(json.dumps(ke, indent=0))
    print(f"kettle annual max: {len(ke)} years, {min(ke)}-{max(ke)}")
    gr = parse_granby()
    (OUT / "granby_annual_max.json").write_text(json.dumps(gr, indent=0))
    miss = [y for y in range(Y0, 2025) if y not in gr]
    print(f"granby annual max: {len(gr)} years, missing {Y0}-2024: {miss}")


if __name__ == "__main__":
    main()

Notesynth.py: the synthesizers

"""
Synthesis layer for the electronic-ambient pieces (2026-09-25): band-limited
oscillators, a slowly swept low-pass filter, electronic drums, filtered
noise, tempo-synced ping-pong delay and sidechain ducking. Everything is
vectorized numpy (no per-sample Python loops) so a 7-8 minute render stays
well under a minute. Outputs are (n, 2) float arrays at sampler.SR, the same
shape the sampled Instruments return, so synth and samples share buses.

Pieces:
  osc(kind, f, n)            - 'sine' | 'saw' | 'square' | 'tri', polyBLEP saws/squares,
                               f may be a per-sample array (glides, vibrato)
  supersaw(f, n, voices)     - detuned saw stack, stereo-spread
  swept_lowpass(x, cutoff)   - 2nd-order low-pass whose cutoff follows an array
  env_adsr(n, a, d, s, r)    - linear-attack, exponential-decay envelope
  pad_note / sub_note / pluck_note / drone
  kick / hat / clap / rim    - electronic kit voices
  noise_band(n, lo, hi)      - band-passed noise (wind, rain, storm rumble)
  pingpong(x, delay_s, fb)   - stereo feedback delay, each repeat darker
  sidechain(n, hits, depth)  - gain curve that ducks after each kick
"""
import numpy as np
from scipy.signal import butter, sosfilt, sosfilt_zi

from sampler import SR


def midi_hz(m):
    return 440.0 * 2 ** ((np.asarray(m, dtype=float) - 69) / 12)


def _as_arr(f, n):
    f = np.asarray(f, dtype=float)
    return np.full(n, float(f)) if f.ndim == 0 else f[:n]


def _polyblep(t, dt):
    """polyBLEP residual for a phase t in [0,1) with increment dt (vectorized)."""
    y = np.zeros_like(t)
    a = t < dt
    x = t[a] / dt[a]
    y[a] = x + x - x * x - 1
    b = t > 1 - dt
    x = (t[b] - 1) / dt[b]
    y[b] = x * x + x + x + 1
    return y


def osc(kind, f, n, phase0=0.0):
    f = _as_arr(f, n)
    dt = f / SR
    ph = (phase0 + np.cumsum(dt)) % 1.0
    if kind == "sine":
        return np.sin(2 * np.pi * ph)
    if kind == "saw":
        return 2 * ph - 1 - _polyblep(ph, dt)
    if kind == "square":
        sq = np.where(ph < 0.5, 1.0, -1.0)
        return sq + _polyblep(ph, dt) - _polyblep((ph + 0.5) % 1.0, dt)
    if kind == "tri":
        return 2 * np.abs(2 * ph - 1) - 1
    raise ValueError(kind)


def supersaw(f, n, voices=5, detune_cents=14, spread=0.8, rng=None):
    rng = rng or np.random.default_rng(0)
    out = np.zeros((n, 2))
    f = _as_arr(f, n)
    for i in range(voices):
        c = (i - (voices - 1) / 2) / max(1, (voices - 1) / 2) * detune_cents
        x = osc("saw", f * 2 ** (c / 1200), n, rng.uniform())
        pan = (i - (voices - 1) / 2) / max(1, (voices - 1) / 2) * spread
        a = (pan + 1) * np.pi / 4
        out[:, 0] += x * np.cos(a)
        out[:, 1] += x * np.sin(a)
    return out / np.sqrt(voices)


def swept_lowpass(x, cutoff, q_order=2, block=1024):
    """Low-pass with a time-varying cutoff (Hz, scalar or per-sample array).
    Coefficients are recomputed per block and the filter state carried over,
    which is click-free for the slow sweeps used here."""
    x = np.atleast_2d(x.T).T if x.ndim == 1 else x
    n = len(x)
    cut = _as_arr(cutoff, n)
    y = np.empty_like(x)
    zi = None
    for s in range(0, n, block):
        c = float(np.clip(cut[min(s + block // 2, n - 1)], 20, SR * 0.45))
        sos = butter(q_order, c, "lowpass", fs=SR, output="sos")
        if zi is None:
            zi = np.stack([sosfilt_zi(sos) * x[0, ch] for ch in range(x.shape[1])], axis=-1)
        seg = x[s:s + block]
        out, zi = sosfilt(sos, seg, axis=0, zi=zi)
        y[s:s + block] = out
    return y


def bandpass(x, lo, hi, order=2):
    sos = butter(order, [lo, hi], "bandpass", fs=SR, output="sos")
    return sosfilt(sos, x, axis=0)


def highpass(x, hz, order=2):
    return sosfilt(butter(order, hz, "highpass", fs=SR, output="sos"), x, axis=0)


def env_adsr(n, a=0.01, d=0.2, s=0.7, r=0.4, hold=None):
    """Attack (linear) -> decay (exp) to sustain -> release (exp) over n samples.
    `hold` = seconds before release starts (default: n minus release)."""
    t = np.arange(n) / SR
    hold = (n / SR - r) if hold is None else hold
    e = np.where(t < a, t / max(a, 1e-4), s + (1 - s) * np.exp(-(t - a) / max(d, 1e-4)))
    rel = t >= hold
    if rel.any():
        e0 = e[np.argmax(rel)]
        e[rel] = e0 * np.exp(-(t[rel] - hold) / max(r / 4, 1e-4))
    return e


def stereo(x, width=0.0, delay_ms=9):
    """mono -> (n,2); width>0 adds a short Haas offset on the right."""
    if width <= 0:
        return np.stack([x, x], axis=1)
    d = int(delay_ms * SR / 1000)
    r = np.concatenate([np.zeros(d), x[:-d]]) if d else x
    return np.stack([x, (1 - width) * x + width * r], axis=1)


# -------------------------------------------------------------- voices -----
def pad_note(midi, dur, vel=0.6, cutoff=1200, a=1.5, r=2.5, rng=None, voices=5, detune=12):
    """Warm analogue-style pad: supersaw through a low-pass, slow attack."""
    n = int((dur + r) * SR)
    x = supersaw(midi_hz(midi), n, voices=voices, detune_cents=detune, rng=rng)
    x = swept_lowpass(x, cutoff)
    return x * env_adsr(n, a=a, d=1.0, s=0.85, r=r, hold=dur)[:, None] * vel * 0.35


def sub_note(midi, dur, vel=0.7, r=0.25, drive=1.4):
    """Sine sub with a touch of saturation (adds the 2nd/3rd harmonics that
    make a sub audible on small speakers)."""
    n = int((dur + r) * SR)
    x = np.tanh(drive * osc("sine", midi_hz(midi), n)) / np.tanh(drive)
    return stereo(x * env_adsr(n, a=0.02, d=0.3, s=0.9, r=r, hold=dur) * vel * 0.6)


def pluck_note(midi, vel=0.6, decay=0.35, cutoff=2600, rng=None):
    """Short filtered saw pluck for arpeggios."""
    n = int((decay * 4 + 0.05) * SR)
    x = 0.6 * osc("saw", midi_hz(midi), n) + 0.4 * osc("square", midi_hz(midi) * 1.002, n)
    e = np.exp(-np.arange(n) / SR / decay)
    cut = 300 + cutoff * e ** 0.7  # filter closes with the note
    x = swept_lowpass((x * e)[:, None], cut, block=256)[:, 0]
    return stereo(x * vel * 0.4, width=0.3)


def drone(midis, n, cutoff, level, rng=None, breathe=None, voices=3, detune=7):
    """Long sustained drone over n samples. cutoff/level: per-sample arrays.
    breathe: optional per-sample 0..1 amplitude modulation (e.g. an annual cycle)."""
    rng = rng or np.random.default_rng(1)
    x = np.zeros((n, 2))
    for m in midis:
        x += supersaw(midi_hz(m), n, voices=voices, detune_cents=detune, spread=0.6, rng=rng)
        x += 0.5 * stereo(osc("sine", midi_hz(m - 12), n))  # sub-octave body
    x = swept_lowpass(x, cutoff, block=2048)
    amp = _as_arr(level, n)
    if breathe is not None:
        amp = amp * (0.75 + 0.25 * _as_arr(breathe, n))
    return x * amp[:, None] * 0.25 / max(1, len(midis))


# ---------------------------------------------------------------- drums ----
def kick(vel=0.8, f0=52, f_start=140, decay=0.45, click=0.25):
    n = int((decay * 3) * SR)
    t = np.arange(n) / SR
    f = f0 + (f_start - f0) * np.exp(-t / 0.035)
    x = np.sin(2 * np.pi * np.cumsum(f) / SR) * np.exp(-t / decay)
    x += click * np.random.default_rng(3).standard_normal(n) * np.exp(-t / 0.003)
    return stereo(np.tanh(1.5 * x) * vel * 0.8)


def hat(vel=0.3, open_=False, rng=None):
    rng = rng or np.random.default_rng(5)
    d = 0.22 if open_ else 0.045
    n = int(d * 4 * SR)
    x = highpass(rng.standard_normal(n), 7000, order=4) * np.exp(-np.arange(n) / SR / d)
    return stereo(x * vel * 0.5, width=0.4)


def clap(vel=0.4, rng=None):
    rng = rng or np.random.default_rng(6)
    n = int(0.35 * SR)
    t = np.arange(n) / SR
    env = sum(np.exp(-np.maximum(t - o, 0) / 0.006) * (t >= o) for o in (0, 0.011, 0.022))
    env = env + 0.6 * np.exp(-np.maximum(t - 0.03, 0) / 0.09) * (t >= 0.03)
    x = bandpass(rng.standard_normal(n), 900, 3200) * env
    return stereo(x * vel * 0.6, width=0.5)


def rim(vel=0.35):
    n = int(0.12 * SR)
    t = np.arange(n) / SR
    x = (np.sin(2 * np.pi * 1650 * t) + 0.5 * np.sin(2 * np.pi * 540 * t)) * np.exp(-t / 0.012)
    return stereo(x * vel * 0.5)


def noise_band(n, lo, hi, rng=None, order=2):
    """Stereo band-passed noise (independent channels -> wide)."""
    rng = rng or np.random.default_rng(8)
    x = rng.standard_normal((n, 2))
    return bandpass(x, lo, hi, order)


# --------------------------------------------------------------- effects ---
def pingpong(x, delay_s, fb=0.45, repeats=6, damp_hz=3500, mix=0.35):
    """Stereo ping-pong echo: repeats alternate sides and get darker. Returns
    only the wet signal (add it to the dry)."""
    n = len(x)
    d = int(delay_s * SR)
    wet = np.zeros_like(x)
    tap = x.mean(axis=1)
    sos = butter(1, damp_hz, "lowpass", fs=SR, output="sos")
    for k in range(1, repeats + 1):
        tap = sosfilt(sos, tap)
        off = k * d
        if off >= n:
            break
        side = (k + 1) % 2
        wet[off:, side] += tap[: n - off] * (fb ** k)
    return wet * mix


def sidechain(n, hit_samples, depth=0.45, release=0.18):
    """Gain curve: dips by `depth` at each kick and recovers exponentially."""
    g = np.zeros(n)
    k = np.exp(-np.arange(int(release * 5 * SR)) / SR / release)
    for h in hit_samples:
        h = int(h)
        if 0 <= h < n:
            m = min(len(k), n - h)
            g[h:h + m] = np.maximum(g[h:h + m], k[:m])
    return 1 - depth * g


# ------------------------------------------------ warm voices (2026-09-25) --
def saturate(x, drive=1.3):
    """Soft tape-style saturation: tanh with unity gain at low level."""
    return np.tanh(drive * x) / np.tanh(drive)


def ep_note(midi, dur, vel=0.5, bright=0.5, r=0.6):
    """Electric-piano (Rhodes-ish) voice: 2-operator FM, modulator at 1x with
    a decaying index (the bell-like 'tine' attack), plus a quiet 14x tine."""
    n = int((dur + r) * SR)
    t = np.arange(n) / SR
    f = float(midi_hz(midi))
    idx = (0.6 + 2.2 * bright * vel) * np.exp(-t / 0.35)
    mod = np.sin(2 * np.pi * f * t)
    x = np.sin(2 * np.pi * f * t + idx * mod)
    x += 0.08 * vel * np.sin(2 * np.pi * f * 14 * t) * np.exp(-t / 0.05)
    amp = env_adsr(n, a=0.004, d=0.9, s=0.35, r=r, hold=dur)
    trem = 1 + 0.12 * np.sin(2 * np.pi * 4.5 * t)                 # stereo tremolo
    y = x * amp * vel * 0.45
    return np.stack([y * trem, y * (2 - trem)], axis=1)


VOWELS = {"oh": (450, 800, 2830), "ah": (730, 1090, 2440), "oo": (325, 700, 2530),
          "eh": (530, 1840, 2480), "ee": (300, 2200, 3000)}


def vox_chop(midi, dur, vowel="oh", vel=0.5, bend=0.0, rng=None):
    """Synthesized vocal chop: a glottal-ish pulse (saw with vibrato) through
    three formant band-passes. `bend` semitones glide into the note, like a
    chopped and repitched vocal. Not a sample of any real voice."""
    rng = rng or np.random.default_rng(12)
    n = int((dur + 0.08) * SR)
    t = np.arange(n) / SR
    f = float(midi_hz(midi)) * 2 ** ((bend * np.exp(-t / 0.05)) / 12) * (1 + 0.006 * np.sin(2 * np.pi * 5.5 * t))
    src = osc("saw", f, n) + 0.05 * rng.standard_normal(n)
    y = np.zeros(n)
    for fc, g in zip(VOWELS[vowel], (1.0, 0.7, 0.35)):
        bw = 0.12 * fc + 60
        y += g * bandpass(src, max(40, fc - bw), fc + bw, order=2)
    e = env_adsr(n, a=0.008, d=0.12, s=0.8, r=0.06, hold=dur)
    return stereo(y * e * vel * 0.9, width=0.25)


def lead_note(midi, dur, vel=0.5, cutoff=3000, voices=4, detune=10, rng=None):
    """Bright-but-warm detuned lead for arpeggios: supersaw, fixed low-pass,
    quick attack, medium decay."""
    n = int((dur + 0.5) * SR)
    x = supersaw(midi_hz(midi), n, voices=voices, detune_cents=detune, spread=0.7, rng=rng)
    x = sosfilt(butter(2, cutoff, "lowpass", fs=SR, output="sos"), x, axis=0)
    e = env_adsr(n, a=0.005, d=0.25, s=0.45, r=0.45, hold=dur)
    return x * e[:, None] * vel * 0.4


def split_normal_cdf(x, mode, sd_lo, sd_hi):
    """CDF of a two-piece normal (different spread either side of the mode),
    used to approximate an asymmetric posterior from its median and 95% interval."""
    from scipy.stats import norm
    x = np.asarray(x, dtype=float)
    w_lo = sd_lo / (sd_lo + sd_hi)
    return np.where(x < mode, 2 * w_lo * norm.cdf(x, mode, sd_lo),
                    w_lo + (1 - w_lo) * (2 * norm.cdf(x, mode, sd_hi) - 1))

To rebuild: from scripts/, run python fetch_climate.py, then python climate_pair.py pacific and python climate_pair.py atlantic, then python climate_check_plot.py pacific (and atlantic) and python export_climate_guide.py, then python export_climate_viz.py for the visualizer data. The instrument samples come from fetch_samples.py, as for Water Year Rings.