"""Reproduce a hypothetical uncoupled oscillator-bank magnitude illustration.

Run: python reproduce.py. Requires existing numpy and matplotlib installations.
No device properties, microphone measurements, coupling or electrical model.
"""
from pathlib import Path
import csv
import json
import platform
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

OUT = Path(__file__).resolve().parent
NATURAL_HZ = np.array([250., 500., 1000.])
DAMPING = [0.1, 0.5]
FREQUENCY_HZ = np.geomspace(50., 2000., 1601)

def magnitude(f, fn, zeta):
    r = np.asarray(f) / np.asarray(fn)
    return 1 / np.sqrt((1-r*r)**2 + (2*zeta*r)**2)

def main():
    plt.rcParams.update({"font.size": 10, "svg.fonttype": "none", "pdf.fonttype": 42})
    fig, axes = plt.subplots(1, 2, figsize=(11.5, 6), sharey=True)
    rows = []
    colors = ["#0072B2", "#D55E00", "#009E73"]
    for ax, zeta in zip(axes, DAMPING):
        h = magnitude(FREQUENCY_HZ[:, None], NATURAL_HZ[None, :], zeta)
        rms = np.sqrt(np.mean(h**2, axis=1))
        assert np.isfinite(h).all() and np.isfinite(rms).all()
        assert np.allclose(magnitude(NATURAL_HZ, NATURAL_HZ, zeta), 1/(2*zeta))
        assert np.allclose(magnitude(0., NATURAL_HZ, zeta), 1.)
        for i, (fn, color) in enumerate(zip(NATURAL_HZ, colors)):
            ax.plot(FREQUENCY_HZ, h[:, i], color=color, lw=1.7, label=f"Channel: $f_n$ = {fn:g} Hz")
        ax.plot(FREQUENCY_HZ, rms, color="#202124", lw=2.2, ls="--", label="RMS across three channels")
        ax.axhline(1., color="#777777", lw=.8, ls=":")
        ax.set(xscale="log", xlim=(50,2000), ylim=(0,5.4), xlabel="Forcing frequency, f (Hz; logarithmic scale)", title=fr"Damping ratio $\zeta$ = {zeta:g}")
        ax.set_xticks([50,100,250,500,1000,2000], ["50","100","250","500","1000","2000"])
        ax.grid(alpha=.2)
        for f, hs, a in zip(FREQUENCY_HZ, h, rms):
            rows.append([f,zeta,*hs,a])
    fig.legend(*axes[0].get_legend_handles_labels(), loc="upper center", bbox_to_anchor=(.53,.91), ncol=4, frameon=False, fontsize=9)
    axes[0].set_ylabel("Displacement magnitude / each channel's static response")
    fig.suptitle("Illustrative response of three uncoupled damped oscillators", fontsize=15, y=.96)
    fig.subplots_adjust(left=.08, right=.98, top=.77, bottom=.28, wspace=.09)
    fig.text(.08,.18, r"Model: $H_i=[(1-r_i^2)^2+(2\zeta r_i)^2]^{-1/2}$, $r_i=f/f_{n,i}$; aggregate: $\sqrt{(H_1^2+H_2^2+H_3^2)/3}$.", fontsize=10)
    fig.text(.08,.115,"Horizontal dotted reference: H = 1 (static-response magnitude). All parameters are hypothetical.\nRMS is a descriptive channel-magnitude statistic; it is not a coherent voltage sum or measured microphone performance.", fontsize=9, linespacing=1.6)
    fig.text(.08,.055,"Assumptions: linear steady state, independent harmonic forcing, no coupling, no transduction or acoustic loading.", fontsize=9)
    for ext in ["svg", "pdf", "png"]:
        fig.savefig(OUT / f"resonator-bank.{ext}", dpi=180, facecolor="white")
    plt.close(fig)
    with (OUT/"resonator-bank.csv").open("w",newline="",encoding="utf-8") as stream:
        writer=csv.writer(stream)
        writer.writerow(["forcing_frequency_hz","damping_ratio","H_fn250hz","H_fn500hz","H_fn1000hz","rms_channel_magnitude"])
        writer.writerows(rows)
    methods={
        "title":"Illustrative uncoupled damped oscillator bank",
        "status":"Hypothetical mathematical illustration, not measured or predicted device performance",
        "equation_of_motion":"m*x'' + c*x' + k*x = F0*cos(2*pi*f*t)",
        "definitions":{"natural_frequency_hz":"sqrt(k/m)/(2*pi)","damping_ratio":"c/(2*sqrt(k*m))","static_response":"F0/k for each channel","frequency_ratio":"r_i=f/fn_i"},
        "magnitude_formula":"H_i=1/sqrt((1-r_i^2)^2+(2*zeta*r_i)^2)",
        "aggregate_formula":"sqrt(mean(H_i^2)) across the three channels at each frequency",
        "natural_frequencies_hz":NATURAL_HZ.tolist(),
        "damping_ratios":DAMPING,
        "grid":{"spacing":"logarithmic (numpy.geomspace)","start_hz":50,"stop_hz":2000,"points":1601,"includes_endpoints":True},
        "csv_rows":len(rows),
        "normalization":"Each displacement magnitude is divided by that channel's own static displacement F0/k; equal normalized forcing is not equal absolute displacement or voltage.",
        "assumptions":["linear time-invariant steady state","three independent uncoupled oscillators","same damping ratio assigned to all channels in each panel","harmonic forcing","no noise or transients"],
        "excluded":["coherent sum and phase interference","electrical transduction","coupling","microphone sensitivity","acoustic loading","actual source parameters","actual device geometry or performance"],
        "checks":{"all_outputs_finite":True,"H_at_zero_equals_one":True,"H_at_fn_equals_one_over_2zeta":True},
        "software":{"python":platform.python_version(),"numpy":np.__version__,"matplotlib":matplotlib.__version__},
        "reproduction":"python reproduce.py",
        "figure_alt":"Two panels compare dimensionless response versus frequency for three uncoupled channels at 250, 500 and 1000 Hz. Damping 0.1 produces narrow peaks near 5; damping 0.5 produces broad low peaks. The dashed black curve is the RMS of the three normalized magnitudes, not an electrical sum."
    }
    (OUT/"methods.json").write_text(json.dumps(methods,indent=2)+"\n",encoding="utf-8")
    print(f"Wrote {len(rows)} CSV rows and SVG/PDF/PNG figures to {OUT}")

if __name__ == "__main__":
    main()