import numpy as np
import matplotlib.pyplot as plt

# ---------------------------
# CONSTANTS
# ---------------------------

R = 8.314        # J/mol/K
T_fixed = 308.15  # K

# ---------------------------
# PAPER PARAMETERS (CABEQ)
# ---------------------------

Ea = 35300      # J/mol
T_star = 414.6  # K

KM = 3.21e-3    # mol/L
KES1 = 7.57e-7  # mol/L
KES2 = 1.27e-8  # mol/L
KP = 1.22e-2    # mol/L

# ---------------------------
# ENZYME
# ---------------------------

# 1 g/L enzyme in same units as model
urease_1gL = 1 

# ---------------------------
# REACTION RATE FUNCTION
# ---------------------------

def reaction_rate(pH, S_sat=0.1, P_prod=0.0, urease=urease_1gL):
    """
    Reaction rate at given pH under substrate-saturated conditions
    and no product yet
    """
    H = 10**(-pH)
    S = max(S_sat, 1e-16)
    P = max(P_prod, 1e-16)

    kT = np.exp(-Ea/R * (1/T_fixed - 1/T_star))

    den_substrate = KM + S
    den_product = 1 + P / KP
    den_pH = 1 + (10**(-pH))/KES1 + KES2/(10**(-pH))

    ve = kT * urease * S / (den_substrate * den_product * den_pH)
    return ve

# ---------------------------
# CALCULATE RATE VS pH
# ---------------------------

pH_values = np.linspace(4, 9, 300)
rates = np.array([reaction_rate(pH) for pH in pH_values])

# ---------------------------
# PLOT (styled like your other figure)
# ---------------------------

plt.figure(figsize=(4.7,4))
plt.plot(pH_values, rates, color='black', lw=2)
plt.xlabel('pH / -', fontsize=14)
plt.ylabel('$r$ / mol l$^{-1}$ min$^{-1}$', fontsize=14)
plt.grid(False)

# Make full box (all spines visible)
ax = plt.gca()
for spine in ax.spines.values():
    spine.set_visible(True)
    spine.set_linewidth(1)
    spine.set_color('black')

plt.xticks([4, 5, 6, 7, 8, 9], fontsize = 12 )  # only these pH values on x-axis
plt.yticks([0, 0.005, 0.010, 0.015, 0.020, 0.025], fontsize = 12 )  # rate axis
plt.tight_layout()

# Make full box (all spines visible)
ax = plt.gca()
for spine in ax.spines.values():
    spine.set_visible(True)
    spine.set_linewidth(1)
    spine.set_color('black')

# ---------------------------
# ADD VERTICAL DASHED LINES
# ---------------------------

ax.axvline(
    x=6.8,
    color='black',
    linestyle='--',
    linewidth=1.2
)

ax.axvline(
    x=7.2,
    color='black',
    linestyle='--',
    linewidth=1.2
)
ymax = ax.get_ylim()[1]

ax.text(6.8, -ymax*0.05, '6.8', ha='right', va='top', fontsize=9)
ax.text(7.2, -ymax*0.05, '7.2', ha='left', va='top', fontsize=9)
# ---------------------------
# TICKS
# ---------------------------

plt.xticks([4, 5, 6, 7, 8, 9], fontsize=12)
plt.yticks([0, 0.005, 0.010, 0.015, 0.020, 0.025], fontsize=12)

plt.tight_layout()

plt.savefig("responseph.pdf", format='pdf', bbox_inches='tight')

# ---------------------------
# FIND MAXIMUM RATE
# ---------------------------

max_index = np.argmax(rates)           # index of max rate
pH_max = pH_values[max_index]          # corresponding pH
rate_max = rates[max_index]            # maximum rate

print(f"Maximum reaction rate: {rate_max:.6f} mol L^-1 min^-1 at pH = {pH_max:.2f}")



plt.show()