import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Lecture du fichier (séparateur tabulation, certaines lignes ont moins de colonnes)
data = pd.read_csv("BINOME10.txt", sep="\t", skiprows=1, header=None)

# Les 6 premières colonnes sont inutiles.
# Ensuite : paires (EA2=x, courbe=y) dans l'ordre :
# cols 6-7 : 0_0
# cols 8-9 : 10_10
# cols 10-11 : 15_10
# cols 12-13 : 20_10
# cols 14-15 : 25_10
# cols 16-17 : 30_10 (mal étiqueté 25_10{2} dans le fichier)

def concentration(v_mis):
    """Concentration en Fe(CN)6 en mmol/L. v_mis en mL."""
    v_total = v_mis + 80 + 10
    return 0.1 * v_mis / v_total * 1000  # en mmol/L

curves = [
    (0,                   6,  7),
    (concentration(10),   8,  9),
    (concentration(15),  10, 11),
    (concentration(20),  12, 13),
    (concentration(25),  14, 15),
    (concentration(30),  16, 17),
]

def to_potential(v):
    """Sortie potentiostat (V) -> Potentiel réel (mV). 0V<-> -2000mV, 2.5V<->0mV, 5V<-> +2000mV"""
    return (v - 2.5) * 800

def to_current(v):
    """Sortie potentiostat (V) -> Courant réel (mA). 0V<-> -20mA, 2.5V<->0mA, 5V<-> +20mA"""
    return (v - 2.5) * 8

# --- Graphe 1 : toutes les courbes ---
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

concentrations = []
i_diffusion = []
i_sigma = []

for conc, col_x, col_y in curves:
    x = to_potential(pd.to_numeric(data.iloc[:, col_x], errors="coerce"))
    y = to_current(pd.to_numeric(data.iloc[:, col_y], errors="coerce"))

    ax1.plot(x, y, label=f"{conc:.1f} mmol/L")

    # Moyenne et écart-type du courant sur le palier de diffusion [500, 1000] mV
    mask = (x >= 500) & (x <= 1000)
    i_mean = y[mask].mean()
    i_std  = y[mask].std(ddof=1)
    concentrations.append(conc)
    i_diffusion.append(i_mean)
    i_sigma.append(i_std)

ax1.set_xlabel("Potentiel (mV)")
ax1.set_ylabel("Courant (mA)")
ax1.set_title("Courbes Intensité Potentiel — \n(légende : [Fe(CN)₆] en mmol/L)")
ax1.legend()
ax1.grid(True, alpha=0.3)

# --- Graphe 2 : courant de diffusion vs concentration ---
concentrations = np.array(concentrations)
i_diffusion = np.array(i_diffusion)
i_sigma = np.array(i_sigma)

# Régression linéaire forcée à l'origine I = a·c, pondérée par 1/sigma^2
# Solution analytique : a = sum(I_i * c_i / sigma_i^2) / sum(c_i^2 / sigma_i^2)
w = 1.0 / i_sigma**2
a = np.sum(w * i_diffusion * concentrations) / np.sum(w * concentrations**2)
sigma_a = 1.0 / np.sqrt(np.sum(w * concentrations**2))

# Chi2 réduit : N points, 1 paramètre libre → ddl = N - 1
chi2 = np.sum(w * (i_diffusion - a * concentrations)**2)
ddl  = len(concentrations) - 1
chi2_red = chi2 / ddl

c_fit = np.linspace(0, concentrations.max(), 200)

ax2.errorbar(concentrations, i_diffusion, yerr=i_sigma,
             fmt="o", color="tab:blue", capsize=4, zorder=5, label="Données ± σ")
ax2.plot(c_fit, a * c_fit, color="tab:red", linestyle="--",
         label=f"I = ({a:.4f} ± {sigma_a:.4f})·c")
ax2.text(0.05, 0.95,
         f"χ²_red = {chi2_red:.2f}  (ddl = {ddl})",
         transform=ax2.transAxes, va="top",
         bbox=dict(boxstyle="round", facecolor="wheat", alpha=0.5))
ax2.set_xlabel("[Fe(CN)₆] (mmol/L)")
ax2.set_ylabel("Courant de diffusion (mA)")
ax2.set_title("Courant de diffusion vs concentration\n(palier 500–1000 mV, modèle I = a·c)")
ax2.legend()
ax2.grid(True, alpha=0.3)

print(f"Pente  a  = {a:.4f} ± {sigma_a:.4f} mA·L/mmol")
print(f"χ²_red    = {chi2_red:.3f}  (ddl = {ddl})")
if chi2_red < 2:
    print("→ Modèle linéaire valide (χ²_red ≈ 1).")
else:
    print("→ χ²_red élevé : le modèle ou les incertitudes sont à revoir.")

plt.tight_layout()
plt.show()