"""Reproduce a bounded synthetic Gaussian chromatography teaching figure.

Run: python reproduce.py. Requires numpy and matplotlib. Outputs beside script.
No experimental chromatogram or chemical identification is represented.
"""
from pathlib import Path
import csv
import json
import math
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

OUT = Path(__file__).resolve().parent
SIGMA = 10.0
MU1 = 120.0
TIME = np.linspace(0.0, 240.0, 2401)
DELTAS = (20.0, 40.0, 60.0)
BLUE, MAGENTA, NAVY = "#0072B2", "#C23A91", "#142C52"

def density(t, mu):
    return np.exp(-0.5 * ((t - mu) / SIGMA) ** 2) / (SIGMA * np.sqrt(2 * np.pi))

def main():
    plt.rcParams.update({"font.size": 12, "axes.titlesize": 14,
                         "axes.labelsize": 12, "svg.fonttype": "none",
                         "pdf.fonttype": 42})
    fig, axes = plt.subplots(3, 1, figsize=(10, 11), sharex=True, sharey=True)
    records, checks = [], []
    for ax, delta in zip(axes, DELTAS):
        mu2 = MU1 + delta
        first, second = density(TIME, MU1), density(TIME, mu2)
        total = first + second
        rs = delta / (4 * SIGMA)
        a1, a2 = float(np.trapezoid(first, TIME)), float(np.trapezoid(second, TIME))
        assert abs(a1 - 1) < 2e-9 and abs(a2 - 1) < 2e-9
        assert rs == 2 * (mu2 - MU1) / (4 * SIGMA + 4 * SIGMA)
        exact_window2 = 0.5 * (math.erf((240-mu2)/(SIGMA*math.sqrt(2))) - math.erf((0-mu2)/(SIGMA*math.sqrt(2))))
        checks.append({"separation_s": delta, "resolution_Rs": rs,
                       "peak_1_window_area_trapezoid": a1,
                       "peak_2_window_area_trapezoid": a2,
                       "peak_2_window_area_analytic": exact_window2,
                       "sum_window_area_trapezoid": float(np.trapezoid(total, TIME))})
        ax.plot(TIME, total, color=NAVY, lw=2.6, label="Sum of both peaks")
        ax.plot(TIME, first, color=BLUE, lw=2, ls="--", label="Peak 1: center 120 s")
        ax.plot(TIME, second, color=MAGENTA, lw=2, ls=":", label=f"Peak 2: center {mu2:g} s")
        ax.set_title(f"Center separation {delta:g} s  |  Resolution Rs = {rs:g}", loc="left", pad=10)
        ax.set_ylabel("Normalized density (1/s)")
        ax.set_ylim(0, 0.056)
        ax.set_yticks([0, 0.02, 0.04])
        ax.set_xlim(0, 240)
        ax.grid(alpha=0.22)
        ax.spines[["top", "right"]].set_visible(False)
        ax.legend(loc="upper left", fontsize=10, frameon=False)
        for t, p1, p2, summed in zip(TIME, first, second, total):
            records.append((delta, SIGMA, MU1, mu2, rs, float(t), float(p1), float(p2), float(summed)))
    axes[-1].set_xlabel("Time (s)")
    axes[-1].set_xticks(np.arange(0, 241, 40))
    fig.suptitle("Two synthetic Gaussian chromatography peaks", fontsize=19, x=0.11, ha="left", y=0.985)
    fig.text(0.11, 0.949, "Each individual peak has unit area; both have standard deviation 10 s.", fontsize=12)
    fig.text(0.11, 0.035, "Rs = separation / (4 × standard deviation). Width 4σ uses extrapolated tangents,\n"
             "not finite peak support. Ideal curves only: no measured noise, asymmetry, baseline,\n"
             "drift or chemical identification. Density is not concentration or detector voltage.", fontsize=10, linespacing=1.5)
    fig.subplots_adjust(top=0.90, bottom=0.14, left=0.11, right=0.97, hspace=0.40)
    fig.savefig(OUT / "chromatography-peaks.svg")
    fig.savefig(OUT / "chromatography-peaks.pdf")
    fig.savefig(OUT / "chromatography-peaks.png", dpi=160)
    plt.close(fig)
    with (OUT / "chromatography-peaks.csv").open("w", newline="", encoding="utf-8") as f:
        writer = csv.writer(f)
        writer.writerow(["separation_s", "sigma_s", "peak_1_center_s", "peak_2_center_s", "resolution_Rs", "time_s", "peak_1_density_per_s", "peak_2_density_per_s", "sum_density_per_s"])
        writer.writerows(records)
    methods = {
        "title": "Synthetic equal-area Gaussian chromatography peaks",
        "status": "Constructed teaching model, not an experimental chromatogram",
        "equations": {"individual_density": "p(t; mu,sigma)=exp(-(t-mu)^2/(2*sigma^2))/(sigma*sqrt(2*pi))",
                      "sum": "p1(t)+p2(t)", "resolution": "Rs=2*(mu2-mu1)/(wb1+wb2)=delta/(4*sigma)",
                      "baseline_width_convention": "wb=4*sigma: tangents at Gaussian inflection points intercept the zero baseline at mu +/- 2*sigma; Gaussian support remains infinite"},
        "inputs": {"peak_1_center_s": MU1, "peak_2_center_s": [MU1+d for d in DELTAS], "sigma_s": SIGMA,
                   "individual_peak_area_on_real_line": 1, "sum_area_on_real_line": 2,
                   "separations_s": list(DELTAS), "baseline_width_each_s": 40},
        "sampling": {"start_s": 0, "end_s": 240, "step_s": 0.1, "points_per_case": len(TIME), "csv_data_rows": len(records)},
        "units": "Time in seconds; normalized density in inverse seconds. Neither concentration nor detector voltage.",
        "limits": ["No actual chromatograph is modeled.", "No measured noise, asymmetry, baseline or drift is included.",
                   "No chemical identity, selectivity, retention chemistry or real-instrument performance is inferred.",
                   "Individual Gaussian areas are one on the entire real line; the finite plot omits negligible tails.",
                   "The tangent baseline width is a convention, not finite support.",
                   "Rs=1.5 is a geometric value here; it does not imply zero overlap or validated separation of real analytes."],
        "validation": {"area_absolute_tolerance": 2e-9, "integration": "NumPy trapezoidal rule on the plotted grid; analytic Gaussian finite-window area also reported", "cases": checks},
        "software": {"numpy": np.__version__, "matplotlib": matplotlib.__version__},
        "reproduction": "Run python reproduce.py using installed numpy and matplotlib. Outputs are written beside the script."
    }
    (OUT / "methods.json").write_text(json.dumps(methods, indent=2) + "\n", encoding="utf-8")
    print(json.dumps({"rows": len(records), "area_checks": checks}, indent=2))

if __name__ == "__main__":
    main()