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

# Work in progress, it is far from trivial to make ``quantities`` and ``sympy`` play nice together

# In[ ]:


import math
from chempy import Equilibrium, ReactionSystem
from chempy.units import SI_base_registry, Backend, default_units as u, default_constants as c
from chempy.thermodynamics import GibbsEqConst
from chempy.kinetics.rates import MassAction, Arrhenius
from chempy.kinetics.ode import get_odesys


# In[ ]:


DH, DS, R = -20e3 * u.J/u.mol, -30 * u.J/u.mol/u.K, c.molar_gas_constant.definition
ref = 'reinterpreted from: M. W. Lister & D. E. Rivington: Can. J. Chem., 1955, 33(10): 1572-1590'
eq = Equilibrium({'Fe+3', 'SCN-'}, {'FeSCN+2'}, GibbsEqConst([DH/R, DS/R]))


# In[ ]:


be = Backend()
eq.equilibrium_constant({'temperature': 298.15*u.K}, backend=be)


# In[ ]:


A, Ea = math.exp(35.5)/u.molar/u.s, 72.2e3 * u.J/u.mol
fw, bw = eq.as_reactions(kf=MassAction(Arrhenius([A, Ea/R])))
fw.param.dedimensionalisation(SI_base_registry)


# In[ ]:


params = {'temperature': 298.15*u.K, 'Fe+3': 1e-2*u.molar, 'SCN-': 1e-3*u.molar, 'FeSCN+2': 1e-6*u.molar}
fw.rate(params, backend=be)


# In[ ]:


bw.param.dedimensionalisation(SI_base_registry)


# In[ ]:


bw.rate(params, backend=be)  # note incorrect formula in bw case (need to divide by c0)


# In[ ]:


rsys = ReactionSystem([fw, bw])
rsys


# In[ ]:


# Still failing
#odesys = get_odesys(rsys, units=SI_base_registry)[0]