"""
Isomerisierung von EZ-Kongorot (Z) zu EE-Kongorot (E), pH 10

Mechanismus
  (1) Z + OH⁻  ⇌  Z⁻ + H₂O      k₁, k₋₁   (schnelles Vorgleichgewicht, K₁ = k₁/k₋₁)
  (2) Z⁻  →  E⁻                   k₂        (geschwindigkeitsbestimmend)
  (3) E⁻ + H₂O  →  E + OH⁻        k₃        (schnell, E⁻ quasistationär)
  (4) Z  →  E                     k₀        (unkatalysierte Parallelreaktion)

Ergebnis (Vorgleichgewicht + Quasistationarität von E⁻):
  k_obs = (k₀ + K₁·k₂·[OH]) / (1 + K₁·[OH])
  [Z]   = [Z₀]·e^(−k_obs·t) / (1 + K₁[OH])
  [Z⁻]  = [Z₀]·K₁[OH]·e^(−k_obs·t) / (1 + K₁[OH])
  [E]   = [Z₀]·(1 − e^(−k_obs·t))
"""


import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# ---------------------------------------------------------------- Parameter
pH = 10
OH = 10.0 ** (pH - 14)     # mol/L  (25 °C, pKw = 14)  ->  1e-4 mol/L
Z0 = 1.0e-5                # mol/L  Anfangskonzentration EZ-Kongorot

k0 = 2.0e-4                # 1/s         unkatalysierte Isomerisierung (langsam)
k1 = 1.0e8                 # L/(mol s)   Deprotonierung (schnell)
K1 = 3.0e4                 # L/mol       Gleichgewichtskonstante der Deprotonierung
km1 = k1 / K1              # 1/s         Rückreaktion k₋₁ = k₁/K₁
k2 = 2.0e-3                # 1/s         Isomerisierung des Anions (RDS)

SHOW_NUMERIC_CHECK = False  # True: exakte numerische Lösung gestrichelt einzeichnen

# ------------------------------------------------- analytische Lösung
a = K1 * OH                                   # K₁·[OH]
k_obs = (k0 + K1 * k2 * OH) / (1 + a)

t = np.linspace(0, 5 / k_obs, 1000)           # s, ca. 5 Zeitkonstanten
decay = np.exp(-k_obs * t)

Z = Z0 * decay / (1 + a)                      # Edukt (protoniert)
Zm = Z0 * a * decay / (1 + a)                 # Intermediat (Anion)
E = Z0 * (1 - decay)                          # Produkt

# ---------------------------- Kontrolle: volles System ohne Näherung
def rhs(_t, y):
    z, zm, e = y
    dz = -(k0 + k1 * OH) * z + km1 * zm
    dzm = k1 * OH * z - (km1 + k2) * zm
    de = k0 * z + k2 * zm
    return [dz, dzm, de]

# Start im Gleichgewicht von Z und Z⁻ (Vorgleichgewicht ist in µs eingestellt)
y0 = [Z0 / (1 + a), Z0 * a / (1 + a), 0.0]
sol = solve_ivp(rhs, (t[0], t[-1]), y0, t_eval=t,
                method="Radau", rtol=1e-8, atol=1e-14)

dev = max(np.max(np.abs(sol.y[0] - Z)),
          np.max(np.abs(sol.y[1] - Zm)),
          np.max(np.abs(sol.y[2] - E))) / Z0
print(f"[OH] = {OH:.1e} mol/L,  K₁[OH] = {a:.2f},  k_obs = {k_obs:.3e} 1/s")
print(f"Halbwertszeit t½ = {np.log(2) / k_obs:.0f} s")
print(f"Max. Abweichung Näherung vs. exakte Lösung: {dev:.2e} (relativ zu [Z₀])")

# ----------------------------------------------------------------- Plot
sc = 1e6  # mol/L -> µmol/L
fig, ax = plt.subplots(figsize=(9, 5))

ax.plot(t, Z * sc, color="tab:blue", lw=2, label="[Z] EZ-Kongorot (Edukt)")
ax.plot(t, Zm * sc, color="tab:orange", lw=2, label="[Z⁻] Anion (Intermediat)")
ax.plot(t, E * sc, color="tab:green", lw=2, label="[E] EE-Kongorot (Produkt)")

if SHOW_NUMERIC_CHECK:
    for y, c in zip(sol.y, ("tab:blue", "tab:orange", "tab:green")):
        ax.plot(t, y * sc, ls="--", color="k", lw=1)

ax.set_xlabel("Zeit t / s")
ax.set_ylabel("Konzentration c / µmol L⁻¹")
ax.set_title(f"Isomerisierung EZ- → EE-Kongorot, pH {pH}")
ax.set_xlim(t[0], t[-1])
ax.set_ylim(0, Z0 * sc * 1.05)
ax.grid(alpha=0.3)

# Legende außerhalb des Diagramms -> schneidet garantiert keine Kurve
ax.legend(loc="upper left", bbox_to_anchor=(1.02, 1.0), frameon=True)
fig.tight_layout()

fig.savefig("kongorot_kinetik.png", dpi=200, bbox_inches="tight")
plt.show()