{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "from collections import defaultdict\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import chempy\n",
    "from chempy.chemistry import Species, Equilibrium\n",
    "from chempy.equilibria import EqSystem\n",
    "print('ChemPy: %s' % chempy.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "substances = list(map(Species.from_formula,\n",
    "    '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()))        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_conc = defaultdict(lambda: 1e-50, {'H+': .1, 'OH-': 1e-7, 'HSCN': 1e-3, 'Fe+3': 1e-3, 'H2O': 55.4})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H2O_c = init_conc['H2O']\n",
    "w_autop = Equilibrium({'H2O': 1}, {'H+': 1, 'OH-': 1}, 10**-14/H2O_c)\n",
    "HSCN_pr = Equilibrium({'HSCN': 1}, {'H+': 1, 'SCN-': 1}, 10**1.28, ref={'doi': '10.1139/v00-152'})\n",
    "FeOH = Equilibrium({'Fe+3': 1, 'H2O': 1}, {'FeOH+2': 1, 'H+': 1}, 10**-2.774 / H2O_c, ref={'doi': '10.1039/B001811M'})\n",
    "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'})\n",
    "FeSCN_B1 = Equilibrium({'Fe+3': 1, 'SCN-': 1}, {'FeSCN+2': 1}, 10**2.065, ref={'doi': '10.1039/B001811M'})\n",
    "FeSCN_B2 = Equilibrium({'Fe+3': 1, 'SCN-': 2}, {'Fe(SCN)2+': 1}, 10**3.34, ref={'doi': '10.1351/pac199769071489'})\n",
    "FeSCN_B3 = Equilibrium({'Fe+3': 1, 'SCN-': 3}, {'Fe(SCN)3': 1}, 10**3.82, ref={'doi': '10.1351/pac199769071489'})\n",
    "FeSCN_B4 = Equilibrium({'Fe+3': 1, 'SCN-': 4}, {'Fe(SCN)4-': 1}, 10**3.6, ref={'doi': '10.1351/pac199769071489'})\n",
    "FeSCN_B5 = Equilibrium({'Fe+3': 1, 'SCN-': 5}, {'Fe(SCN)5-2': 1}, 10**4.3, ref={'doi': '10.1351/pac199769071489'})\n",
    "FeSCN_B6 = Equilibrium({'Fe+3': 1, 'SCN-': 6}, {'Fe(SCN)6-3': 1}, 10**5, ref={'doi': '10.1351/pac199769071489'})\n",
    "FeOOH_S = Equilibrium({'FeOOH(s)': 1, 'H+': 3}, {'Fe+3': 1, 'H2O': 2}, 10**0.4*H2O_c**2, ref='Jmvkt tabell')\n",
    "equilibria = w_autop, HSCN_pr, FeOH, Fe2OH2, FeSCN_B1, FeSCN_B2, FeSCN_B3, FeSCN_B4, FeSCN_B5, FeSCN_B6, FeOOH_S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eqsys = chempy.equilibria.EqSystem(equilibria, substances)\n",
    "eqsys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SCN_varied = np.logspace(-6, 0)\n",
    "plt.figure(figsize=(16, 6))\n",
    "from chempy.equilibria import NumSysLog, NumSysLin\n",
    "NumSysLog.small = 1e-80\n",
    "xout, info, sanity = eqsys.roots(init_conc, SCN_varied, 'SCN-', plot_kwargs={'latex_names': True}, NumSys=(NumSysLog, NumSysLin))\n",
    "plt.gca().set_xscale('log')\n",
    "plt.gca().set_yscale('log')\n",
    "plt.gca().set_ylim([1e-5, 1e-2])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}