#!/usr/bin/env python3
"""
The Sound of the Sri Yantra -- rendered literally.

From the essay's verdict (ayurastro.com/writings/the-sound-of-the-sri-yantra.html):
"one tone, struck twenty-seven times, from twenty-seven bearings,
 each strike placing one line by phase alone. Not a chord but a litany."

Mapping (all from the solved classical figure, Huet bases):
  * ONE frequency for every strike: 136.10 Hz (the 'Om' tone of the legend in section 1).
  * 27 strikes = 9 bases + 18 legs, in drawing order t1..t9 (base, left leg, right leg).
  * Stereo pan  = edge bearing:  pan = sin(2*theta)   (bases centre, diagonals wide).
  * Loudness and ring time = the edge's FULL CHORD across the unit plate (rule one:
    every line prints rim to rim).
  * Starting phase of each strike = 2*pi * (signed distance of the line from centre)
    -- the line really is "placed by phase alone".
  * The bindu is a point, not a line: no wave can draw it. It stays silent.

Output: 44.1 kHz stereo WAV (convert to MP3 with ffmpeg).
Tested with numpy only; Python 3.9+ compatible.
"""
from __future__ import annotations
import numpy as np
import wave

# ----------------------------------------------------------------------
# 1. The solved classical figure (section 4 table): base y, half-width, apex y
# ----------------------------------------------------------------------
TRIANGLES = [
    # (name, base_y, half_width, apex_y)   t1-t5 point down (shakti), t6-t9 up (shiva)
    ("t1", +0.774, 0.520, -0.074),
    ("t2", +0.538, 0.649, -0.670),
    ("t3", +0.336, 0.942, -1.000),
    ("t4", +0.208, 0.325, -0.470),
    ("t5", +0.103, 0.228, -0.204),
    ("t6", -0.074, 0.320, +0.538),
    ("t7", -0.204, 0.979, +1.000),
    ("t8", -0.470, 0.714, +0.774),
    ("t9", -0.670, 0.466, +0.336),
]

def edge_params(p1, p2):
    """Bearing (mod pi), full-chord length across unit circle, signed offset from centre."""
    p1 = np.asarray(p1, float); p2 = np.asarray(p2, float)
    u = p2 - p1
    u = u / np.hypot(*u)
    theta = np.arctan2(u[1], u[0]) % np.pi          # direction, mod 180 degrees
    s = p1[0] * u[1] - p1[1] * u[0]                  # signed distance of line from origin
    chord = 2.0 * np.sqrt(max(0.0, 1.0 - s * s))     # rule one: the full chord prints
    return theta, chord, s

EDGES = []
for name, yb, w, ya in TRIANGLES:
    EDGES.append((f"{name} base",      edge_params((-w, yb), (+w, yb))))
    EDGES.append((f"{name} left leg",  edge_params((-w, yb), (0.0, ya))))
    EDGES.append((f"{name} right leg", edge_params((+w, yb), (0.0, ya))))
assert len(EDGES) == 27

# ----------------------------------------------------------------------
# 2. Synthesis
# ----------------------------------------------------------------------
FS      = 44100
F0      = 136.10          # the one tone (change to taste: 108.0, 220.0, ...)
GAP     = 2.4             # seconds between strikes -- litany pace
LEAD    = 0.5
TAIL    = 8.0
rng     = np.random.default_rng(43)   # 43 small triangles

n_total = int(FS * (LEAD + GAP * len(EDGES) + TAIL))
L = np.zeros(n_total)
R = np.zeros(n_total)

for k, (label, (theta, chord, s)) in enumerate(EDGES):
    tau  = 1.1 + 1.5 * (chord / 2.0)              # long chords ring longer
    dur  = min(9.0, 5.0 * tau)
    n    = int(FS * dur)
    t    = np.arange(n) / FS

    env  = (1.0 - np.exp(-t / 0.004)) * np.exp(-t / tau)
    ph   = 2.0 * np.pi * s                        # position carried by phase alone

    # the pair of waves at +/-beta: same frequency, slightly different phase L/R
    beta = 0.35 * np.sin(theta)
    toneL = np.sin(2 * np.pi * F0 * t + ph - beta)
    toneR = np.sin(2 * np.pi * F0 * t + ph + beta)

    # soft mallet transient (the strike itself, not a tone)
    nzlen = int(FS * 0.06)
    nz = rng.standard_normal(nzlen) * np.exp(-np.arange(nzlen) / (FS * 0.012))
    nz = np.convolve(nz, np.hanning(96) / np.hanning(96).sum(), mode="same")
    mallet = np.zeros(n); mallet[:nzlen] = nz * 0.35

    amp  = 0.30 + 0.70 * (chord / 2.0)            # loudness = chord length
    pan  = np.sin(2.0 * theta)                    # bearing -> stereo field
    gl, gr = np.sqrt((1 - pan) / 2), np.sqrt((1 + pan) / 2)

    sigL = amp * (env * toneL + mallet)
    sigR = amp * (env * toneR + mallet)

    i0 = int(FS * (LEAD + k * GAP))
    i1 = min(i0 + n, n_total)
    L[i0:i1] += gl * sigL[: i1 - i0]
    R[i0:i1] += gr * sigR[: i1 - i0]

# gentle master fade and normalise to -1 dBFS
fade = int(FS * 1.5)
ramp = np.linspace(1, 0, fade)
L[-fade:] *= ramp; R[-fade:] *= ramp
peak = max(np.abs(L).max(), np.abs(R).max())
L *= 0.891 / peak; R *= 0.891 / peak

stereo = np.empty(2 * n_total, dtype=np.int16)
stereo[0::2] = (L * 32767).astype(np.int16)
stereo[1::2] = (R * 32767).astype(np.int16)

with wave.open("sri_yantra.wav", "wb") as f:
    f.setnchannels(2); f.setsampwidth(2); f.setframerate(FS)
    f.writeframes(stereo.tobytes())

print(f"{len(EDGES)} strikes, one tone at {F0} Hz, "
      f"{n_total / FS:.1f} s -> sri_yantra.wav")
for label, (theta, chord, s) in EDGES:
    print(f"  {label:14s} bearing {np.degrees(theta):6.1f}  chord {chord:.3f}  offset {s:+.3f}")