#!/usr/bin/env python
# coding: utf-8

# In[ ]:


import matplotlib.pyplot as plt
import sympy as sp
from chempy import Substance, Equilibrium, ReactionSystem
from chempy.kinetics.ode import get_odesys
get_ipython().run_line_magic('matplotlib', 'inline')


# In[ ]:


K1, K3, R1f, R3f = sp.symbols('K_1 K_3 R_1f R_3f')

def mono_rev_single():
    '''Reversible mechanism with one central complex (pg 36)
    '''
    E, A, P, X = map(Substance, 'EAPX')
    equilibria = e1, e3 = [
        Equilibrium({'E', 'A'}, {'X'}, param=K1),
        Equilibrium({'X'}, {'E', 'P'}, param=K3)
    ]
    reactions = e1.as_reactions(kf=R1f) + e3.as_reactions(kf=R3f)
    return ReactionSystem(reactions, (E, A, P, X))


# In[ ]:


rsys = mono_rev_single()


# In[ ]:


rsys


# In[ ]:


type(rsys.rxns[1].param)


# In[ ]:


odesys = get_odesys(rsys, include_params=False, params=True)[0]


# In[ ]:


assert len(odesys.params) == 4
odesys.params


# In[ ]:


odesys.get_jac()


# In[ ]:


c0 = {'E': 5, 'A': 2, 'P': 0, 'X': 0}
params = {K1: 17, K3: 23, R1f: 63, R3f: 43}


# In[ ]:


odesys.f_cb(42, [c0[k] for k in rsys.substances], [params[k] for k in odesys.params])


# In[ ]:


odesys.names, odesys.param_names  # get_odesys or SymbolicSys could determine names from Symbol.name


# In[ ]:


result = odesys.integrate(42, c0, params, integrator='cvode')


# In[ ]:


result.plot(xscale='log')
plt.gca().set_xlim([1e-6, plt.gca().get_xlim()[1]])
_ = plt.legend(loc='best')