#!/usr/bin/env python3
"""Sample Experiment: frequency response of a first-order RC low-pass filter.

Computes |H(f)| and phase for H(f) = 1 / (1 + j*2*pi*f*R*C) over 10 Hz - 1 MHz,
writes the dataset as JSONL plus a sidecar metadata JSON, and renders a Bode plot.
Fully deterministic: re-running this script reproduces the dataset bit-for-bit.
"""
import json, math, hashlib, datetime, pathlib

R_OHMS = 1000.0        # 1 kOhm  (design ref: Yageo RC0603FR-071KL)
C_FARADS = 100e-9      # 100 nF  (design ref: Murata GRM188R71H104KA93D)
F_CUTOFF = 1.0 / (2.0 * math.pi * R_OHMS * C_FARADS)   # ~1591.55 Hz

POINTS_PER_DECADE = 20
F_START_EXP, F_STOP_EXP = 1, 6   # 10 Hz .. 1 MHz

root = pathlib.Path(__file__).resolve().parent.parent
data_path = root / "data" / "rc_lowpass_response.jsonl"
rows = []
n = (F_STOP_EXP - F_START_EXP) * POINTS_PER_DECADE + 1
for i in range(n):
    f = 10 ** (F_START_EXP + i / POINTS_PER_DECADE)
    w_rc = 2.0 * math.pi * f * R_OHMS * C_FARADS
    gain = 1.0 / math.sqrt(1.0 + w_rc * w_rc)
    rows.append({
        "point": i,
        "frequency_hz": round(f, 6),
        "gain_linear": round(gain, 9),
        "gain_db": round(20.0 * math.log10(gain), 6),
        "phase_deg": round(-math.degrees(math.atan(w_rc)), 6),
    })

with open(data_path, "w") as fh:
    for r in rows:
        fh.write(json.dumps(r) + "\n")

sha = hashlib.sha256(data_path.read_bytes()).hexdigest()
sidecar = {
    "dataset": "rc_lowpass_response.jsonl",
    "title": "RC low-pass filter frequency response (computed)",
    "description": "Magnitude and phase of a first-order RC low-pass filter (R=1 kOhm, C=100 nF, fc~=1591.5 Hz), 20 points/decade from 10 Hz to 1 MHz.",
    "created_utc": "2026-07-10T00:00:00Z",
    "generator": {"script": "analysis/run_experiment.py", "runtime": "python3", "deterministic": True},
    "method": "closed-form transfer function H(f) = 1/(1 + j*2*pi*f*R*C); no physical measurement",
    "parameters": {"R_ohms": R_OHMS, "C_farads": C_FARADS, "f_cutoff_hz": round(F_CUTOFF, 3),
                   "points_per_decade": POINTS_PER_DECADE, "f_range_hz": [10, 1000000]},
    "schema": {
        "format": "jsonl",
        "fields": [
            {"name": "point", "type": "int", "unit": None, "description": "0-based sample index"},
            {"name": "frequency_hz", "type": "float", "unit": "Hz", "description": "stimulus frequency"},
            {"name": "gain_linear", "type": "float", "unit": "V/V", "description": "|H(f)|"},
            {"name": "gain_db", "type": "float", "unit": "dB", "description": "20*log10(|H(f)|)"},
            {"name": "phase_deg", "type": "float", "unit": "degrees", "description": "arg(H(f)), negative = lag"},
        ],
    },
    "rows": len(rows),
    "sha256": sha,
    "license": "CC-BY-4.0",
}
(root / "data" / "rc_lowpass_response.sidecar.json").write_text(json.dumps(sidecar, indent=2) + "\n")
print(f"wrote {len(rows)} rows, fc={F_CUTOFF:.2f} Hz, sha256={sha[:16]}...")

# Bode plot
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

freqs = [r["frequency_hz"] for r in rows]
gains = [r["gain_db"] for r in rows]
phases = [r["phase_deg"] for r in rows]

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(9, 6.5), sharex=True)
ax1.semilogx(freqs, gains, color="#1a4480", lw=2)
ax1.axvline(F_CUTOFF, color="#c05621", ls="--", lw=1, label=f"fc = {F_CUTOFF:.1f} Hz (-3 dB)")
ax1.axhline(-3.0103, color="#c05621", ls=":", lw=1)
ax1.set_ylabel("Gain (dB)"); ax1.grid(True, which="both", alpha=0.3); ax1.legend()
ax1.set_title("RC low-pass frequency response - R=1 k$\\Omega$, C=100 nF")
ax2.semilogx(freqs, phases, color="#1a4480", lw=2)
ax2.axvline(F_CUTOFF, color="#c05621", ls="--", lw=1)
ax2.set_ylabel("Phase (deg)"); ax2.set_xlabel("Frequency (Hz)"); ax2.grid(True, which="both", alpha=0.3)
fig.tight_layout()
fig.savefig(root / "docs" / "bode-plot.png", dpi=110)
print("plot saved")