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

# In[ ]:


from __future__ import print_function, division
from collections import defaultdict
import numpy as np
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')
import chempy
from chempy.chemistry import Species, Equilibrium
from chempy.equilibria import EqSystem
print('ChemPy: %s' % chempy.__version__)


# In[ ]:


substances = list(map(Species.from_formula,
    'H+ OH- HSCN SCN- H2O Fe+3 FeSCN+2 Fe(SCN)2+ Fe(SCN)3 Fe(SCN)4- Fe(SCN)5-2 Fe(SCN)6-3 Fe2(OH)2+4 FeOH+2 FeOOH(s)'.split()))        


# In[ ]:


init_conc = defaultdict(lambda: 1e-50, {'H+': .1, 'OH-': 1e-7, 'HSCN': 1e-3, 'Fe+3': 1e-3, 'H2O': 55.4})


# In[ ]:


H2O_c = init_conc['H2O']
w_autop = Equilibrium({'H2O': 1}, {'H+': 1, 'OH-': 1}, 10**-14/H2O_c)
HSCN_pr = Equilibrium({'HSCN': 1}, {'H+': 1, 'SCN-': 1}, 10**1.28, ref={'doi': '10.1139/v00-152'})
FeOH = Equilibrium({'Fe+3': 1, 'H2O': 1}, {'FeOH+2': 1, 'H+': 1}, 10**-2.774 / H2O_c, ref={'doi': '10.1039/B001811M'})
Fe2OH2 = Equilibrium({'Fe+3': 2, 'H2O': 2}, {'Fe2(OH)2+4': 1, 'H+': 2}, 10**-2.812 / H2O_c**2, ref={'doi': '10.1039/B001811M'})
FeSCN_B1 = Equilibrium({'Fe+3': 1, 'SCN-': 1}, {'FeSCN+2': 1}, 10**2.065, ref={'doi': '10.1039/B001811M'})
FeSCN_B2 = Equilibrium({'Fe+3': 1, 'SCN-': 2}, {'Fe(SCN)2+': 1}, 10**3.34, ref={'doi': '10.1351/pac199769071489'})
FeSCN_B3 = Equilibrium({'Fe+3': 1, 'SCN-': 3}, {'Fe(SCN)3': 1}, 10**3.82, ref={'doi': '10.1351/pac199769071489'})
FeSCN_B4 = Equilibrium({'Fe+3': 1, 'SCN-': 4}, {'Fe(SCN)4-': 1}, 10**3.6, ref={'doi': '10.1351/pac199769071489'})
FeSCN_B5 = Equilibrium({'Fe+3': 1, 'SCN-': 5}, {'Fe(SCN)5-2': 1}, 10**4.3, ref={'doi': '10.1351/pac199769071489'})
FeSCN_B6 = Equilibrium({'Fe+3': 1, 'SCN-': 6}, {'Fe(SCN)6-3': 1}, 10**5, ref={'doi': '10.1351/pac199769071489'})
FeOOH_S = Equilibrium({'FeOOH(s)': 1, 'H+': 3}, {'Fe+3': 1, 'H2O': 2}, 10**0.4*H2O_c**2, ref='Jmvkt tabell')
equilibria = w_autop, HSCN_pr, FeOH, Fe2OH2, FeSCN_B1, FeSCN_B2, FeSCN_B3, FeSCN_B4, FeSCN_B5, FeSCN_B6, FeOOH_S


# In[ ]:


eqsys = chempy.equilibria.EqSystem(equilibria, substances)
eqsys


# In[ ]:


SCN_varied = np.logspace(-6, 0)
plt.figure(figsize=(16, 6))
from chempy.equilibria import NumSysLog, NumSysLin
NumSysLog.small = 1e-80
xout, info, sanity = eqsys.roots(init_conc, SCN_varied, 'SCN-', plot_kwargs={'latex_names': True}, NumSys=(NumSysLog, NumSysLin))
plt.gca().set_xscale('log')
plt.gca().set_yscale('log')
plt.gca().set_ylim([1e-5, 1e-2])