{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Debye-Hückel theory\n",
    "The [Debye-Hückel equation](https://en.wikipedia.org/wiki/Debye%E2%80%93H%C3%BCckel_equation) describes the non-ideality of dilute electrolytes. In this notebook we will showcase how ChemPy can work *both* with floating points representations (and then optionally with units) and with symbolic input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import log as ln\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from chempy.units import default_constants as consts, default_units as u\n",
    "from chempy.properties.water_density_tanaka_2001 import water_density\n",
    "from chempy.properties.water_permittivity_bradley_pitzer_1979 import water_permittivity\n",
    "import chempy.electrolytes as dh\n",
    "import chempy.symbolic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def DH_A_B(T):\n",
    "    eps_r = water_permittivity(T, 1*u.bar, units=u)\n",
    "    rho = water_density(T, units=u)\n",
    "    b0 = 1 * u.mol/u.kg\n",
    "    A = dh.A(eps_r, T, rho, b0=b0, constants=consts, units=u).simplified/ln(10)\n",
    "    B = dh.B(eps_r, T, rho, b0=b0, constants=consts, units=u).rescale(1/u.angstrom)\n",
    "    return A, B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DH_A_B(298.15*u.kelvin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = (np.linspace(0, 40)+273.15)*u.Kelvin\n",
    "A, B = DH_A_B(T)\n",
    "%matplotlib inline\n",
    "plt.figure(figsize=(18,4))\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.plot(T-273.15*u.K, A)\n",
    "plt.xlabel(r\"Temperature / $^\\circ C$\")\n",
    "plt.ylabel(\"A\")\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.plot(T-273.15*u.K, B)\n",
    "plt.xlabel(r\"Temperature / $^\\circ C$\")\n",
    "plt.ylabel(\"B / $\\AA^{-1}$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "sp.init_printing()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "constants = chempy.symbolic.get_constant_symbols()\n",
    "constants.Avogadro_constant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var_symbs = A, B, eps_r, T, rho, b0, I, z, a, I0, C, gamma = sp.symbols('A B epsilon_r T rho b^o I z a I^o C gamma')\n",
    "one = sp.S(1)\n",
    "var_symbs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A_expr = dh.A(eps_r, T, rho, b0, constants, backend='sympy')\n",
    "sp.Eq(A, A_expr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "B_expr = dh.B(eps_r, T, rho, b0, constants, backend=sp)\n",
    "sp.Eq(B, B_expr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.Eq(sp.log(gamma), dh.limiting_log_gamma(I, z, A, I0=I0, backend=sp))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.Eq(sp.log(gamma), dh.extended_log_gamma(I, z, a, A, B, C=C, I0=I0, backend=sp))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One motivation for having these expressions available as functions in ChemPy is that the more equations we take from (well tested) libraries the smaller is the risk of introducing typos (which can take a long time to spot)."
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}