#!/usr/bin/env python3
"""
phi_chladni.py
Quick synthesizer: generate Chladni-like patterns using φ-recursive modal selection.
Requires: numpy, scipy, matplotlib
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jv, jn_zeros
from math import log
from typing import List, Tuple

phi = (1 + 5**0.5) / 2.0

def log_phi(x: float) -> float:
    return np.log(x) / np.log(phi)

def choose_modes(f: float, f0: float, alpha=1.0, beta=1.0/phi, gamma=0.0, max_n=6, max_m=6, top_k=6):
    E = log_phi(max(f, 1e-12) / f0)
    candidates = []
    for n in range(0, max_n+1):
        for m in range(0, max_m+1):
            score = abs(alpha*n + beta*m + gamma - E)
            candidates.append(((n,m), score))
    candidates.sort(key=lambda x: x[1])
    return [c[0] for c in candidates[:top_k]]

def generate_pattern(f: float = 440.0,
                     f0: float = 110.0,
                     plate_radius: float = 1.0,
                     grid_n: int = 800,
                     modes: List[Tuple[int,int]] = None,
                     radial_zeros_per_mode: int = 3,
                     eta: float = 0.9, zeta: float = 0.6):
    if modes is None:
        modes = choose_modes(f, f0)
    # grid in polar coords
    xs = np.linspace(-plate_radius, plate_radius, grid_n)
    ys = np.linspace(-plate_radius, plate_radius, grid_n)
    X, Y = np.meshgrid(xs, ys)
    R = np.sqrt(X**2 + Y**2)
    TH = np.arctan2(Y, X)
    U = np.zeros_like(R)
    # mask outside circle
    mask = R <= plate_radius

    for (n, m) in modes:
        # pick Bessel zeros (order m)
        # get a few roots and scale them by phi for radial spacing assumption
        # using scipy's jn_zeros to get true zeros, then warp them by phi-decay
        roots = jn_zeros(m, radial_zeros_per_mode)
        for j_idx, root in enumerate(roots):
            # scale radius by phi recursion hypothesis
            scale = phi**(-j_idx)  # r_j = r_plate * phi^{-j}
            k = root / (plate_radius * scale + 1e-12)
            amp = phi**(-(eta * n + zeta * m + 0.5*j_idx))
            phase = (n + m + j_idx) * 0.37  # deterministic-ish phase
            # generate mode contribution, zero outside plate radius
            contribution = amp * jv(m, k * R) * np.cos(m * TH + phase)
            # attenuate outside
            contribution = contribution * mask
            U += contribution

    # clip and normalize for visualization
    U = U * mask
    return X, Y, U

def plot_chladni(X, Y, U, levels=0, title="φ-Chladni pattern"):
    plt.figure(figsize=(6,6))
    # show nodal lines (contour at 0)
    if levels == 0:
        cs = plt.contour(X, Y, U, levels=[0], linewidths=1.0)
    else:
        cs = plt.contour(X, Y, U, levels=levels, linewidths=0.6)
    plt.gca().set_aspect('equal')
    plt.title(title)
    plt.axis('off')
    plt.show()

if __name__ == "__main__":
    f = 440.0    # driving frequency (Hz)
    f0 = 110.0   # base freq
    modes = choose_modes(f, f0, alpha=1.0, beta=1.0/phi, gamma=0.0, max_n=7, max_m=7, top_k=8)
    X, Y, U = generate_pattern(f=f, f0=f0, plate_radius=1.0, grid_n=800, modes=modes,
                               radial_zeros_per_mode=4, eta=0.9, zeta=0.6)
    plot_chladni(X, Y, U, title=f"φ-Chladni f={f}Hz modes={modes}")
