diff --git a/.gitignore b/.gitignore
index 52b635e92ad15fd50b706c63232d10df9e005dab..1c5289f640ee61bb56cdcd74eec4be819a4bae5b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -8,3 +8,5 @@
/teqp
/dev/docker/valgrind/output
/dev/docker/gcov/output
+/**/.ipynb_checkpoints
+/doc/build
\ No newline at end of file
diff --git a/.readthedocs.yml b/.readthedocs.yml
new file mode 100644
index 0000000000000000000000000000000000000000..268e96a2b03c462fde88f9d0cae930094199f1aa
--- /dev/null
+++ b/.readthedocs.yml
@@ -0,0 +1,29 @@
+# .readthedocs.yaml
+# Read the Docs configuration file
+# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details
+
+# Required
+version: 2
+
+# Set the version of Python and other tools you might need
+build:
+ os: ubuntu-20.04
+ tools:
+ python: "3.10"
+ # You can also specify other tool versions:
+ # nodejs: "16"
+ # rust: "1.55"
+ # golang: "1.17"
+
+# Build documentation in the docs/ directory with Sphinx
+sphinx:
+ configuration: doc/source/conf.py
+
+# If using Sphinx, optionally build your docs in additional formats such as PDF
+# formats:
+# - pdf
+
+# Optionally declare the Python requirements required to build your docs
+python:
+ install:
+ - requirements: doc/requirements.txt
\ No newline at end of file
diff --git a/doc/Makefile b/doc/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..d0c3cbf1020d5c292abdedf27627c6abe25e2293
--- /dev/null
+++ b/doc/Makefile
@@ -0,0 +1,20 @@
+# Minimal makefile for Sphinx documentation
+#
+
+# You can set these variables from the command line, and also
+# from the environment for the first two.
+SPHINXOPTS ?=
+SPHINXBUILD ?= sphinx-build
+SOURCEDIR = source
+BUILDDIR = build
+
+# Put it first so that "make" without argument is like "make help".
+help:
+ @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
+
+.PHONY: help Makefile
+
+# Catch-all target: route all unknown targets to Sphinx using the new
+# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS).
+%: Makefile
+ @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
diff --git a/doc/make.bat b/doc/make.bat
new file mode 100644
index 0000000000000000000000000000000000000000..dc1312ab09ca6fb0267dee6b28a38e69c253631a
--- /dev/null
+++ b/doc/make.bat
@@ -0,0 +1,35 @@
+@ECHO OFF
+
+pushd %~dp0
+
+REM Command file for Sphinx documentation
+
+if "%SPHINXBUILD%" == "" (
+ set SPHINXBUILD=sphinx-build
+)
+set SOURCEDIR=source
+set BUILDDIR=build
+
+%SPHINXBUILD% >NUL 2>NUL
+if errorlevel 9009 (
+ echo.
+ echo.The 'sphinx-build' command was not found. Make sure you have Sphinx
+ echo.installed, then set the SPHINXBUILD environment variable to point
+ echo.to the full path of the 'sphinx-build' executable. Alternatively you
+ echo.may add the Sphinx directory to PATH.
+ echo.
+ echo.If you don't have Sphinx installed, grab it from
+ echo.https://www.sphinx-doc.org/
+ exit /b 1
+)
+
+if "%1" == "" goto help
+
+%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O%
+goto end
+
+:help
+%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O%
+
+:end
+popd
diff --git a/doc/requirements.txt b/doc/requirements.txt
new file mode 100644
index 0000000000000000000000000000000000000000..6061e1b6f5e8d4d6fb5224f43510fcfc38a18509
--- /dev/null
+++ b/doc/requirements.txt
@@ -0,0 +1,8 @@
+breathe
+teqp
+numpy
+jupyter
+nbconvert
+nbsphinx
+pandas
+matplotlib
\ No newline at end of file
diff --git a/doc/source/algorithms/VLE.ipynb b/doc/source/algorithms/VLE.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..85dc68c2fcbc01d10f13702bfeb6ccece9c9ba6c
--- /dev/null
+++ b/doc/source/algorithms/VLE.ipynb
@@ -0,0 +1,324 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "8218498b",
+ "metadata": {},
+ "source": [
+ "# Phase equilibria\n",
+ "\n",
+ "Two basic approaches are implemented in teqp:\n",
+ "\n",
+ "* Iterative calculations given guess values\n",
+ "* Tracing along iso-curves (constant temperature, etc.) powered by the isochoric thermodynamics formalism"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "c0b4e863",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import teqp\n",
+ "import numpy as np\n",
+ "import pandas\n",
+ "import matplotlib.pyplot as plt\n",
+ "teqp.__version__"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e57e532b",
+ "metadata": {},
+ "source": [
+ "## Iterative Phase Equilibria"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1baf38e1",
+ "metadata": {},
+ "source": [
+ "### Pure fluid\n",
+ "\n",
+ "For a pure fluid, phase equilibrium between two phases is defined by equating the pressures and Gibbs energies in the two phases. This represents a 2D non-linear rootfinding problem. Newton's method can be used for the rootfinding, and in teqp, automatic differentiation is used to obtain the necessary Jacobian matrix so the implementation is quite efficient.\n",
+ "\n",
+ "The method requires guess values, which are the densities of the liquid and vapor densities. In some cases, ancillary or superancillary equations have been developed which provide curves of guess densities as a function of temperature.\n",
+ "\n",
+ "For a pure fluid, you can use the ``pure_VLE_T`` method to carry out the iteration."
+ ]
+ },
+ {
+ "cell_type": "raw",
+ "id": "f0ca0b22",
+ "metadata": {
+ "raw_mimetype": "text/restructuredtext"
+ },
+ "source": [
+ "The Python method is here: :py:meth:`pure_VLE_T <teqp.teqp.pure_VLE_T>`"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "2674227c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Instantiate the model\n",
+ "model = teqp.canonical_PR([300], [4e6], [0.1])\n",
+ "\n",
+ "T = 250 # [K], Temperature to be used\n",
+ "\n",
+ "# Here we use the superancillary to get guess values (actually these are more \n",
+ "# accurate than the results we will obtain from iteration!)\n",
+ "rhoL0, rhoV0 = model.superanc_rhoLV(T)\n",
+ "display('guess:', [rhoL0, rhoV0])\n",
+ "\n",
+ "# Carry out the iteration, return the liquid and vapor densities\n",
+ "# The guess values are perturbed to make sure the iteration is actually\n",
+ "# changing the values\n",
+ "teqp.pure_VLE_T(model, T, rhoL0*0.98, rhoV0*1.02, 10)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f8805ae1",
+ "metadata": {},
+ "source": [
+ "### Binary Mixture"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "76bccf19",
+ "metadata": {},
+ "source": [
+ "For a binary mixture, the approach is roughly similar to that of a pure fluid. The pressure is equated between phases, and the chemical potentials of each component in each phase are forced to be the same. \n",
+ "\n",
+ "Again, the user is required to provide guess values, in this case molar concentrations in each phase, and a Newton method is implemented to solve for the phase equilibrium. The analytical Jacobian is obtained from automatic differentiation.\n",
+ "\n",
+ "The ``mix_VLE_Tx`` function is the binary mixture analog to ``pure_VLE_T`` for pure fluids."
+ ]
+ },
+ {
+ "cell_type": "raw",
+ "id": "eef189fd",
+ "metadata": {
+ "raw_mimetype": "text/restructuredtext"
+ },
+ "source": [
+ "The Python method is here: :py:meth:`mix_VLE_Tx <teqp.teqp.mix_VLE_Tx>`"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b12bd318",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "zA = np.array([0.01, 0.99])\n",
+ "model = teqp.canonical_PR([300,310], [4e6,4.5e6], [0.1, 0.2])\n",
+ "model1 = teqp.canonical_PR([300], [4e6], [0.1])\n",
+ "T = 273.0 # [K]\n",
+ "# start off at pure of the first component\n",
+ "rhoL0, rhoV0 = model1.superanc_rhoLV(T)\n",
+ "\n",
+ "# then we shift to the given composition in the first phase\n",
+ "# to get guess values\n",
+ "rhovecA0 = rhoL0*zA\n",
+ "rhovecB0 = rhoV0*zA\n",
+ "\n",
+ "# carry out the iteration\n",
+ "code, rhovecA, rhovecB = teqp.mix_VLE_Tx(model, T, rhovecA0, rhovecB0, zA, \n",
+ " 1e-10, 1e-10, 1e-10, 1e-10, # stopping conditions\n",
+ " 10 # maximum number of iterations\n",
+ " )\n",
+ "code, rhovecA, rhovecB"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d72ef08e",
+ "metadata": {},
+ "source": [
+ "You can (and should) check the value of the return code to make sure the iteration succeeded. Do not rely on the numerical value of the enumerated return codes!"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f4e3f914",
+ "metadata": {},
+ "source": [
+ "# Tracing (isobars and isotherms)\n",
+ "\n",
+ "When it comes to mixture thermodynamics, as soon as you add another component to a pure component to form a binary mixture, the complexity of the thermodynamics entirely changes. For that reason, mixture iterative calculations for mixtures are orders of magnitude more difficult to carry out. Asymmetric mixtures can do all sorts of interesting things that are entirely unlike those of pure fluids, and the algorithms are therefore much, much more complicated. Formulating phase equilibrium problems is not much more complicated than for pure fluids, but the most challenging aspect is to obtain good guess values from which to start an iterative routine, and the difficulty of this problem increases with the complexity of the mixture thermodynamics.\n",
+ "\n",
+ "Ulrich Deiters and Ian Bell have developed a number of algorithms for tracing phase equilibrium solutions as the solution of ordinary differential equations rather than carrying out iterative routines for a given state point. The advantage of the tracing calculations is that they can often be initiated at a state point that is entirely known, for instance the pure fluid endpoint for a subcritical isotherm or isobar."
+ ]
+ },
+ {
+ "cell_type": "raw",
+ "id": "e0097771",
+ "metadata": {
+ "raw_mimetype": "text/restructuredtext"
+ },
+ "source": [
+ "The Python method is here: :py:meth:`trace_VLE_isotherm_binary <teqp.teqp.trace_VLE_isotherm_binary>`"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "63902dba",
+ "metadata": {},
+ "source": [
+ "The C++ implementation returns a string in JSON format, which can be conveniently operated upon, for instance after converting the returned data structure to a ``pandas.DataFrame``. A simple example of plotting a subcritical isotherm for a \"boring\" mixture is presented here:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "49dcba2b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "model = teqp.canonical_PR([300,310], [4e6,4.5e6], [0.1, 0.2])\n",
+ "model1 = teqp.canonical_PR([300], [4e6], [0.1])\n",
+ "T = 273.0 # [K]\n",
+ "rhoL0, rhoV0 = model1.superanc_rhoLV(T) # start off at pure of the first component\n",
+ "j = teqp.trace_VLE_isotherm_binary(model, T, np.array([rhoL0, 0]), np.array([rhoV0, 0]))\n",
+ "display(str(j)[0:100]+'...') # The first few bits of the data\n",
+ "df = pandas.DataFrame(j) # Now as a data frame\n",
+ "df.head(3)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9aecca78",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(df['xL_0 / mole frac.'], df['pL / Pa']/1e6)\n",
+ "plt.plot(df['xV_0 / mole frac.'], df['pL / Pa']/1e6)\n",
+ "plt.gca().set(xlabel='$x_1,y_1$ / mole frac.', ylabel='p / MPa')\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "a9c9fe55",
+ "metadata": {},
+ "source": [
+ "Isn't that exciting!\n",
+ "\n",
+ "You can also provide an optional set of flags to the function to control other behaviors of the function, and switch between simple Euler and adaptive RK45 integration (the default)"
+ ]
+ },
+ {
+ "cell_type": "raw",
+ "id": "264c5123",
+ "metadata": {
+ "raw_mimetype": "text/restructuredtext"
+ },
+ "source": [
+ "The options class is here: :py:meth:`TVLEOptions <teqp.teqp.TVLEOptions>`"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d4193110",
+ "metadata": {},
+ "source": [
+ "Supercritical isotherms work approximately in the same manner"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "c5b925ad",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "Tc_K = [190.564, 154.581]\n",
+ "pc_Pa = [4599200, 5042800]\n",
+ "acentric = [0.011, 0.022]\n",
+ "model = teqp.canonical_PR(Tc_K, pc_Pa, acentric)\n",
+ "model1 = teqp.canonical_PR([Tc_K[0]], [pc_Pa[0]], [acentric[0]])\n",
+ "T = 170.0 # [K] # Note: above Tc of the second component\n",
+ "rhoL0, rhoV0 = model1.superanc_rhoLV(T) # start off at pure of the first component\n",
+ "j = teqp.trace_VLE_isotherm_binary(model, T, np.array([rhoL0, 0]), np.array([rhoV0, 0]))\n",
+ "df = pandas.DataFrame(j) # Now as a data frame\n",
+ "plt.plot(df['xL_0 / mole frac.'], df['pL / Pa']/1e6)\n",
+ "plt.plot(df['xV_0 / mole frac.'], df['pL / Pa']/1e6)\n",
+ "plt.gca().set(xlabel='$x_1,y_1$ / mole frac.', ylabel='p / MPa')\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "936c5f55",
+ "metadata": {},
+ "source": [
+ "As of version 0.10.0, isobar tracing has been added to ``teqp``. It operates in fundamentally the same fashion as the isotherm tracing and the same recommendations about starting at a pure fluid apply"
+ ]
+ },
+ {
+ "cell_type": "raw",
+ "id": "81eacaf0",
+ "metadata": {
+ "raw_mimetype": "text/restructuredtext"
+ },
+ "source": [
+ "The tracer function class is here: :py:meth:`trace_VLE_isobar_binary <teqp.teqp.trace_VLE_isobar_binary>`"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "109bc65b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "Tc_K = [190.564, 154.581]\n",
+ "pc_Pa = [4599200, 5042800]\n",
+ "acentric = [0.011, 0.022]\n",
+ "model = teqp.canonical_PR(Tc_K, pc_Pa, acentric)\n",
+ "model1 = teqp.canonical_PR([Tc_K[0]], [pc_Pa[0]], [acentric[0]])\n",
+ "T = 170.0 # [K] # Note: above Tc of the second component\n",
+ "rhoL0, rhoV0 = model1.superanc_rhoLV(T) # start off at pure of the first component\n",
+ "p0 = rhoL0*model1.get_R(np.array([1.0]))*T*(1+model1.get_Ar01(T, rhoL0, np.array([1.0])))\n",
+ "j = teqp.trace_VLE_isobar_binary(model, p0, T, np.array([rhoL0, 0]), np.array([rhoV0, 0]))\n",
+ "df = pandas.DataFrame(j) # Now as a data frame\n",
+ "plt.plot(df['xL_0 / mole frac.'], df['T / K'])\n",
+ "plt.plot(df['xV_0 / mole frac.'], df['T / K'])\n",
+ "plt.gca().set(xlabel='$x_1,y_1$ / mole frac.', ylabel='T / K')\n",
+ "plt.show()"
+ ]
+ }
+ ],
+ "metadata": {
+ "celltoolbar": "Raw Cell Format",
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "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.10.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/doc/source/algorithms/critical_curves.ipynb b/doc/source/algorithms/critical_curves.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..3f4d16ac2886bcb374b01c0a004dc83aa0396ba3
--- /dev/null
+++ b/doc/source/algorithms/critical_curves.ipynb
@@ -0,0 +1,227 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "d88bffed",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "'0.9.4.dev0'"
+ ]
+ },
+ "execution_count": 2,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "import teqp\n",
+ "teqp.__version__"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "8218498b",
+ "metadata": {},
+ "source": [
+ "# Critical curves & points\n",
+ "\n",
+ "\n",
+ "## Pure Fluids\n",
+ "\n",
+ "Solving for the critical point involves finding the temperature and density that make\n",
+ "$$\n",
+ "\\left(\\frac{\\partial p}{\\partial \\rho}\\right)_T = \\left(\\frac{\\partial^2 p}{\\partial \\rho^2}\\right)_T = 0\n",
+ "$$\n",
+ "by 2D non-linear rootfinding. Newton steps are taken, and the analytic Jacobian is used (thanks to the ability to do derivatives with automatic differentiation). This is all handily wrapped up in the\n",
+ "``solve_pure_critical`` method which requires the user to provide guess values for temperature and density"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "46657a96",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Values taken from http://dx.doi.org/10.6028/jres.121.011\n",
+ "modelPR = teqp.canonical_PR([190.564], [4599200], [0.011])\n",
+ "\n",
+ "# Solve for the critical point from a point close to the critical point\n",
+ "T0 = 192.0\n",
+ "# Critical compressibility factor of P-R is 0.307401308698.. (see https://doi.org/10.1021/acs.iecr.1c00847)\n",
+ "rhoc = (4599200/(8.31446261815324*190.564))/0.3074\n",
+ "rho0 = rhoc*1.2345 # Perturb to make sure we are doing something in the solver\n",
+ "teqp.solve_pure_critical(modelPR, T0, rho0)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e15eeded",
+ "metadata": {},
+ "source": [
+ "## Mixtures\n",
+ "\n",
+ "A pure fluid has a single vapor-liquid critical point, but mixtures are different:\n",
+ "\n",
+ "* They may have multiple (or zero!) critical points for a given mixture composition\n",
+ "* The critical curves may not emanate from the pure fluid endpoints\n",
+ "\n",
+ "When it comes to critical points, intuition from pure fluids is not helpful, or sometimes even counter-productive. \n",
+ "\n",
+ "teqp has methods for working with the critical loci of binary mixtures (only binary mixtures, for now) and especially, methods for tracing the critical curves emanating from the pure fluid endpoints.\n",
+ "\n",
+ "The tracing method in teqp is based explicitly on the isochoric thermodynamics formalism introduced by Ulrich Deiters and Sergio Quinones-Cisneros. It uses the Helmholtz energy density as the fundamental potential and all other properties are derived from it. For critical curves it is based upon the integration of sets of ordinary differential equations; the differential equations are in the form of derivatives of the molar concentrations of each component in the mixture with respect to an integration variable. The set of ODE is then integrated.\n",
+ "\n",
+ "Here is an example of the construction of the critical curves emanating from the pure fluid endpoints for the mixture nitrogen + ethane."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "id": "81619ae2",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgYAAAG8CAYAAAClhm0uAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA4zElEQVR4nO3dd3gc5b328e9PvVqyJBfZlhtyxTa2kU2vIWA6hBKMIbRASELqCYG095DknJxDck5OQggBBwi9hYRiQugYMMW9927LTZZsSVaXdp/3j5U3wkjCZaXZnb0/16VrV7PS+vZ4rLn1zMwz5pxDREREBCDB6wAiIiISPVQMREREJEzFQERERMJUDERERCRMxUBERETCVAxEREQkTMVAREREwlQMREREJCzqi4GZjTKz+83seTP7utd5RERE/MyTYmBmD5tZmZktO2D5FDNbbWbrzOxOAOfcSufcrcCVQIkXeUVEROKFVyMGjwBT2i4ws0Tgj8C5wGhgqpmNbn3tImAW8Hb3xhQREYkvnhQD59z7wJ4DFk8G1jnnNjjnmoBngItbv/5l59yJwLTuTSoiIhJfkrwO0EZ/YGubz0uB48zsdOBLQCrwakffbGa3ALcAZGZmHjty5MguCyoi3a++KcC+xhZqWz8cYEBGShJZaUlkpSaRnpKIeR1UxCPz588vd871OtL3iaZi0N7/Z+ecmwnM/Lxvds5NB6YDlJSUuHnz5kU0nIh0H+ccmyvq+GBdOR+uLeej9eXUNrSQAEwu7MHJxfmcVFzA5CF5ZKRE048xEe+Y2eZIvE80/Y8qBYrafD4A2O5RFhHpZuU1jXy0voIP15Yza1052yrrAeiXk8aUMX05qbiAE48qoFd2qsdJRfwtmorBXGCYmQ0BtgFXAVd7G0lEukpdUwtzNu7hw3XlzFpXwcod1QD0SEvixKMKuPX0ozi5uIDB+RmY6QCBSHfxpBiY2dPA6UCBmZUC/+6ce8jMbgNeBxKBh51zy73IJyJdY2N5LW+v3MXbK8uYv3kvTYEgKYkJHDuoJ7efM4KTigsY2z+HxAQVARGveFIMnHNTO1j+Kp2cYCgisaUlEGTe5r3hMrChvBaAEX2yuf6kwZxcXMCkwXmkpyR6nFRE9oumQwlHzMwuBC4sLi72OopI3Kqqa2bmmjLeWVXGzNW7qapvJjnROH5oPtedOJgzR/amKC/D65gi0gFfFQPn3AxgRklJyc1eZxGJJxt21/DOqjLeWrmLuZv2Egg68jNT+OLoPnxhZG9OGd6LrFRf/bgR8S39TxWRQ9bZIYKvnTqUL4zqw/iiXJ0rIBKDVAxE5KDsP0Tw9soyZq4uo7qhhZTEBI4bmqdDBCI+omIgIh3asLuGt1eGDhHM2/yvQwRnH92Xs0b15uRhOkQg4jf6Hy0iYS2BIHM37eWdVZ8+RDCybza3njaUM0fqEIGI3/mqGOiqBJFD19EhguOPyuf6kwZzxggdIhCJJ+ac8zpDxOleCSKd21vbxGvLdzJj8XZmb9xDIOgoyErhjBG9+YIOEYjEJDOb75wrOdL30f98kThR09jCWyt28fLi7by/ZjctQcfQgkxuPa31KoIBuSToEIFI3FMxEPGxhuYAM1fvZsbi7by9ahcNzUH65aRx08lDuPCYfhzdr4fuQyAin6JiIOIzLYEgH66v4OVF23lj+U72NbaQn5nClSVFXHRMPyYO7KmRARHpkIqBiA8Eg455m/fy8uJtvLp0J3tqm8hOS2LKmL5cNL4fJwzNJykxweuYIhIDVAxEYpRzjmXbqnl58TZeWbKDHVUNpCUncNaoPlx4TD9OH9GL1CTdnEhEDo2KgUiMWVe2j5cXbWfGkh1sLK8lOdE4bXgv7jx3JGeN6kOmriYQkSPgq58gmsdA/GrrnjpmLNnOjMU7WLmjmgSDE47K52unDmXKmL7kZqR4HVFEfELzGIhEqbJ9DfxjyQ5mLN7Ogi2VAEwcmMuFx/Tj/HGF9M5O8zagiEQVzWMg4kOBoOPdVWU8OXsz763ZTdDBqMIe/HDKCC4c108zEIpIl1MxEIkCZfsaeG7uVp6es5VtlfX0zk7lG6cXc/H4fgzrk+11PBGJIyoGIh5xzvHJhj08MXszry/bSUvQcVJxPj89fxRnje5Dsi4vFBEPqBiIdLOq+mb+Nr+UJ2dvZv3uWnLSk7nuxMFMO24gQ3tleR1PROKcioFIN1lSWskTn2zm5cXbaWgOMr4ol/+54hguGFdIWrLmGxCR6KBiINKF6psCzFi8nSdmb2ZJaRXpyYlcOqE/044bxJj+OV7HExH5DF8VA81jINFiXdk+nvhkC39bUMq+hhaG9c7i5xcdzaUT+9MjLdnreCIiHfJVMXDOzQBmlJSU3Ox1Fok/TS1B3lixkyc+2cwnG/aQnGicO6aQa44fxKTBPXUXQxGJCb4qBiJeKN1bxzNztvLM3K2U1zQyoGc6P5wygitLiijISvU6nojIIVExEDkMgaDj/TW7eeKTzby7ugwHnDmiN9ccP4hTh/ciUbc1FpEYpWIgcgjKaxp5du5Wnp6zhdK99RRkhSYiumpyEQN6alZCEYl9KgYin8M5x9xNe3n8k828tmwHzQHH8UPzuPPckZw9ui8pSZqISET8Q8VApAPOOd5bs5vfv72WhVsqyU5LYtpxg7jm+IEU99Y0xSLiTyoGIgdwzjFzdagQLNpaSf/cdH55yRgunziA9BRNRCQi/qZiINLKOcc7q8r4/dtrWVJaRf/cdH516VguP3aADheISNxQMZC455zjrZVl3PP2WpZuq6IoL527LxvLpRNUCEQk/qgYSNwKBh1vrNjFPW+vZcWOagblZ/Dry8dx6YT+urOhiMQtXxUDTYksByMYdLy+fCe/f3stq3buY3B+Bv9zxTFcMr4fSSoEIhLnzDnndYaIKykpcfPmzfM6hkSZYNDxz2U7+cM7oUIwtCCT284s5qJjVAhEJPaZ2XznXMmRvo+vRgxE2hMIOl5duoM/vLOWNbtqOKpXJr+/ajwXjOunGQpFRA6gYiC+FQg6XlmynT+8s451ZTUM653FPVMncP7YQhUCEZEOqBiI77QEgryyJDRCsH53LcP7ZHHv1RM4b0whCSoEIiKdUjEQ32gJBHlp0XbufXcdG8trGdk3m/umTWTK0X1VCEREDpKKgfjCR+vL+dmLy1i/u5ZRhT24/5qJnD1ahUBE5FCpGEhM272vkf/8xwpeXLSdgXkZ3H/NsZxzdB/MVAhERA6HioHEpEDQ8dScLfz6tVU0Ngf59pnFfOOMYtKSdS8DEZEjoWIgMWfZtip+8sJSFpdWcVJxPr+4eAxH9cryOpaIiC+oGEjMqG5o5rdvrOGxjzeRl5nK768az0XH9NNhAxGRCFIxkKjnnGPGkh388pUVlNc08pXjB/H9s0eQk57sdTQREd/xVTHQvRL8Z2N5LT97cRmz1pUztn8OD11XwrgBuV7HEhHxLd0rQaJSQ3OA+2au5/6Z60lNSuD2KSOYdtwgzVgoItIB3StBfOu9Nbv5fy8tY3NFHReP78dPzh9F7+w0r2OJiMQFFQOJGjurGvjlKyv4x9IdDC3I5MmvHsdJxQVexxIRiSsqBuK5lkCQxz7ezG/fXENzIMi/fXE4t5w2lNQkzUkgItLdVAzEUwu37OUnLyxjxY5qTh/Ri19cNIaB+RlexxIRiVsqBuKJqrpm7n59FU/P2UKf7DT+NG0iU8b01ZwEIiIeUzGQbuWc44WF2/jPf6yksr6Zm04awne/OJysVG2KIiLRQD+NpdvUNbXw478v5cVF25k4MJcnLh3LqMIeXscSEZE2VAykW6zfXcPXn5jP2rIafnD2cL5xerFuiSwiEoVUDKTLvbp0Bz98fgkpSQk8fuNxnDxMlyCKiEQrFQPpMs2BIHf/cxUPztrIhIG53DdtIoU56V7HEhGRTqgYSJfYVd3AbU8tYO6mvVx/4mB+fN4oUpISvI4lIiKfQ8VAIu6TDRXc9tRC6ppauGfqBC46pp/XkURE5CCpGEjEOOd44P0N/Ob11QzOz+Dpm49jWJ9sr2OJiMghUDGQiKhuaOYHzy3mjRW7OH9cIXdfNk5zE4iIxCBf/eQ2swuBC4uLi72OEldW7qjm60/Mp3RvPf/vgtHccNJgzWAoIhKjfHU2mHNuhnPulpycHK+jxI2/zS/l0vs+pL45wDO3HM+NJw9RKRARiWG+GjGQ7tPQHODnM1bw9JwtnDA0nz9cPYGCrFSvY4mIyBFSMZBDVravga8+Oo8lpVV84/Sj+P4Xh5OU6KvBJxGRuKViIIdke2U90x6cza7qBv78lRK+OLqP15FERCSCVAzkoG2uqOXqP8+muqGZx2+azLGD8ryOJCIiEaZiIAdl7a59THtwNs2BIE/ffDxj+usETxERP1IxkM+1fHsV1z40h8QE49mvncBwTVokIuJbOmNMOrVgy16mTv+E9ORE/qpSICLiexoxkA59vL6Crz46l17ZqTzx1eMY0DPD60giItLFVAykXe+uLuPWx+czMC+DJ796HL17pHkdSUREuoGKgXzGa8t28K2nFzK8TzaP33QceZkpXkcSEZFuomIgn/LCwlJ+8NclHDMgh7/cMJmc9GSvI4mISDdSMZCwFxaW8v3nFnP8kHwevK6ETN0dUUQk7ugnvwAwd9Mefvj8Ek4Yms/D108iLTnR60giIuIBXa4obN1Tx9cen09Rzwz+NO1YlQIRkTimYhDn9jU0c9OjcwkEHQ9dP4mcDJ1TICISz1QM4lgg6Pj20wvZsLuWP02byJCCTK8jiYiIx3SOQRz7r1dX8u7q3fznpWM4sbjA6zgiIhIFNGIQp56Zs4UHZ23k+hMHM+24QV7HERGRKKFiEIc+Xl/BT19cxqnDe/HT80d5HUdERKKIikGc2VxRy9efnM/ggkzuvXoCSYnaBERE5F+0V4gj1Q3N3PjIXAAeuq6EHmm6AkFERD7NV8XAzC40s+lVVVVeR4k6LYEgtz21kM0Vddx/zbEMytcVCCIi8lm+KgbOuRnOuVtycnK8jhJ1fvP6at5fE7oC4fih+V7HERGRKOWrYiDtW7BlL9M/2MDUyQP58qSBXscREZEopmLgc40tAe54fgmFPdL48XkjvY4jIiJRThMc+dwf313P2rIa/nLDJLJ1sqGIiHwOjRj42Mod1dz37jq+NKE/Z4zo7XUcERGJASoGPtUSCHLH35aQk57Mzy4Y7XUcERGJETqU4FMPzdrIktIq7r16Aj0zU7yOIyIiMUIjBj60sbyW3765hrNH9+H8sYVexxERkRiiYuAzwaDjjr8tISUpgV9eMgYz8zqSiIjEEBUDn3lqzhbmbNzDT88fRZ8eaV7HERGRGKNi4COVdU3c/c9VnFScz5UlRV7HERGRGKRi4CMPf7iJfY0t/OyC0TqEICIih0XFwCeqG5p55MONnHN0H0b27eF1HBERiVEqBj7x+MebqW5o4VtnDvM6ioiIxDAVAx+oa2rhwQ82cMaIXozprztLiojI4VMx8IEnP9nC3rpmbtNogYiIHCEVgxjX0Bxg+gcbOKk4n2MH9fQ6joiIxDgVgxj37Nyt7N7XyG1naLRARESOnIpBDGtqCXL/e+uZNLgnxw/N8zqOiIj4gIpBDPv7glJ2VDVw25nDNG+BiIhEhIpBjGoJBLlv5nqOGZDDqcMKvI4jIiI+oWIQo95csYste+r45hnFGi0QEZGIUTGIUTOWbKcgK5UvjOrjdRQREfERFYMYVNvYwjuryjhvbF8SEzRaICIikaNiEIPeXlVGQ3OQC8b18zqKiIj4jIpBDHpl8Xb69EilRBMaiYhIhKkYxJh9Dc3MXLOb88YWkqDDCCIiEmEqBjHmzRW7aGrRYQQREekaKgYx5pUlO+ifm87EgbleRxERER9SMYghVXXNfLB2N+ePK9TcBSIi0iVUDGLI68t30hxwXDCu0OsoIiLiUyoGMeT15TspyktnbP8cr6OIiIhPqRjECOcc87fs5cShBTqMICIiXUbFIEZsqqijsq6ZCTrpUEREulBMFAMzu8TM/mxmL5nZ2V7n8cLCLXsBGK9iICIiXcizYmBmD5tZmZktO2D5FDNbbWbrzOxOAOfci865m4HrgS97ENdzi7ZWkpmSyLDe2V5HERERH/NyxOARYErbBWaWCPwROBcYDUw1s9FtvuSnra/HnYVbKhk3IFc3TRIRkS7lWTFwzr0P7Dlg8WRgnXNug3OuCXgGuNhC7gb+6Zxb0N1ZvdbQHGDljmqdXyAiIl0u2s4x6A9sbfN5aeuybwFnAZeb2a3tfaOZ3WJm88xs3u7du7s+aTdatq2KlqBjfFGu11FERMTnkrwOcID2xsmdc+4e4J7OvtE5Nx2YDlBSUuK6IJtnFm2tBHTioYiIdL1oGzEoBYrafD4A2O5RlqixcEsl/XPT6Z2d5nUUERHxuWgrBnOBYWY2xMxSgKuAlz3O5LkVO6oZN0CzHYqISNfz8nLFp4GPgRFmVmpmNznnWoDbgNeBlcBzzrnlXmWMBs45tlfWM6BnutdRREQkDnh2joFzbmoHy18FXu3mOFGrqr6ZxpYgfXroMIKIiHS9aDuUcETM7EIzm15VVeV1lIjZWd0AoGIgIiLdwlfFwDk3wzl3S06Of47H76wKFYO+OSoGIiLS9XxVDPyorLoRgL4aMRARkW6gYhDl9h9K6N0j1eMkIiISD1QMotzO6gbyMlNITUr0OoqIiMQBFYMot6uqQSceiohIt1ExiHI7qxvoq8MIIiLSTXxVDPx4ueKu6gZdkSAiIt3GV8XAb5crOucor2miIEsjBiIi0j18VQz8Jth6j8jkRP0ziYhI99AeJ4oFXagZJLR3M2oREZEuoGIQxfYXAzM1AxER6R4qBlGstReQoGIgIiLdRMUgiulQgoiIdDdfFQO/Xa4Y1IiBiIh0M18VA79drvivcww8DiIiInHDV8XAb1ww9KgRAxER6S4qBlHMoREDERHpXioGUSwlKfTP09gS9DiJiIjECxWDKJaenEhKUgJ765q8jiIiInFCxSCKmRm56clU1TV7HUVEROKEikGUy81IplLFQEREuomKQZTLTU/RoQQREek2vioGfpvgCEIjBlX1GjEQEZHu4ati4LcJjkCHEkREpHv5qhj4UW5GCpX1OpQgIiLdQ8UgyuVmJNPQHKShOeB1FBERiQMqBlEuLyMFgPKaRo+TiIhIPFAxiHJFeRkAbNlT53ESERGJByoGUW5QfqgYbK5QMRARka6nYhDlCnPSSUlMUDEQEZFuoWIQ5RITjAF56WyuqPU6ioiIxAEVgxgwOD+TTRoxEBGRbuCrYuDHmQ8hdJ7B5opanHNeRxEREZ/zVTHw48yHEBoxqGsKUF6jiY5ERKRr+aoY+NX+KxM26TwDERHpYioGMWB4n2wAVu6o9jiJiIj4nYpBDCjMSaMgK5XFW/117oSIiEQfFYMYYGaMG5DD0m2VXkcRERGfUzGIEWP757CurIbaxhavo4iIiI+pGMSIY4pyCDpYvl3nGYiISNdRMYgRY/vnArCktNLTHCIi4m9Jh/NNZtYTGAak7V/mnHs/UqHks3plp9IvJ40lpToBUUREus4hFwMz+yrwHWAAsAg4HvgYODOiyeQzxg7I0YiBiIh0qcM5lPAdYBKw2Tl3BjAB2B3RVNKucQNy2VRRR1Vds9dRRETEpw6nGDQ45xoAzCzVObcKGBHZWNKecQNCUz0v0WWLIiLSRQ6nGJSaWS7wIvCmmb0EbI9kqMPl15so7TdhYE+SEoyP1ld4HUVERHzqkIuBc+5S51ylc+4u4GfAQ8AlEc51WPx6E6X9slKTmDAwl1lry72OIiIiPnXQxcDM0szsu2Z2r5l9zcySnHPvOededs7ptn/d5OTiXizbXsXeWq1yERGJvEMZMXgUKAGWAucC/9sliaRTJw8rwDn4cL1GDUREJPIOpRiMds5d45x7ALgcOKWLMkknjhmQQ3Zakg4niIhIlziUYhC+Rs45pwn7PZKUmMAJQ/P5YG05zjmv44iIiM8cSjE4xsyqzWyfme0Dxu1/bmaawL8bnTKsgG2V9WyqqPM6ioiI+MxBz3zonEvsyiBy8E4e1guAWWt3M6Qg0+M0IiLiJwddDMzs5c5ed85ddORx5GAMzs+gf246768t59oTBnsdR0REfORQ7pVwAlAKPAXMBqxLEsnnMjO+MKo3z83bSl1TCxkph3UvLBERkc84lHMM+gI/AsYAvwe+CJS3zmXwXleEk46dO6aQhuYgM1frNhUiIhI5B10MnHMB59xrzrnrCN1RcR0w08y+1WXppEOTh+SRn5nCq0t3eB1FRER85JDGoM0sFTgfmAoMBu4B/h75WPJ5EhOMc8b05cWF22hoDpCWrHNDRUTkyB3KlMiPAh8BE4GfO+cmOed+6Zzb1mXppFPnjSmkrinAe2t0OEFERCLjUM4xuBYYDnwH+Kh1ToNqzWPgneOG5tEzI5l/6nCCiIhEyKHMY3A4t2iWLpScmMDZo/vyj6U7aGwJkJqkwwkiInJkfLWzN7MLzWx6VVWV11G6zXnjCqlpbOHtlWVeRxERER/wVTFwzs1wzt2Sk5PjdZRuc3JxAf1z03lq9havo4iIiA/4qhjEo8QE48uTipi1rpzNFbVexxERkRh3KFclnGBmmu0wCl1ZUkRigvH0nK1eRxERkRh3KCMG1wHzzewZM7vezPp2VSg5NH1z0jhzZG+en7+Vppag13FERCSGHcrMh7c65yYCdwE9gUfM7GMz+5WZnWpmOiXeQ1cfN5DymibeWLHT6ygiIhLDDvkcA+fcKufc/znnpgBnArOAKwjdWEk8cuqwXjoJUUREjtgRnXzonKt3zr3qnPuWc64kUqHk0CUmGFMnF/HR+go2luskRBEROTy6KsFHriwpIinBeGaORg1EROTwqBj4SO8eaZw1qg9/nV9KY0vA6zgiIhKDDrkYmNm5ZjbbzFab2XNmdkJXBJPDc/VxA9lT28Try3d5HUVERGLQ4YwY3Ad8HzgemA78xsymRjSVHLaTiwsYmJfBIx9uxDnndRwREYkxh1MMdjnnPnTO7XXOvQWcA/wkwrnkMCUkGDefOpQFWyr5YG2513FERCTGHE4x2GRm/2FmKa2fNwP7IphJjtCVJQPol5PG795ao1EDERE5JIdTDBzwJWCrmc0C1gEzzWxYRJPJYUtNSuSbZxazYEsl72vUQEREDsHhTHA01Tk3GhgEfBf4OZAJPGhmmqw/SlxxbBH9c9P5vzc1aiAiIgfvsC9XdM41OOfmOececs592zl3mnOuKJLh5PClJCVw25nFLNpaycw1u72OIyIiMULzGPjY5ccOYEDPdH6nUQMRETlIKgY+lpyYwLfOLGZxaRXvri7zOo6IiMQAFQOf+9LEAQzMy+B3b63VqIGIiHwuFQOfS04MnWuwpLSKN1ZoNkQREemcikEc+NKE/gzrncUvZqygtrHF6zgiIhLFfFUMzOxCM5teVVXldZSokpSYwK++NJZtlfX835trvI4jIiJRzFfFwDk3wzl3S05OjtdRos6kwXlMnTyQhz/cyLJtKk4iItI+XxUD6dydU0aSl5nKj/6+lJZA0Os4IiIShVQM4khORjL/fuFolm6r4rGPN3sdR0REopCKQZy5YFwhp4/oxf++sZrtlfVexxERkSijYhBnzIxfXjyGgHP8v5eWa24DERH5FBWDOFSUl8H3zhrOWyt38fz8Uq/jiIhIFFExiFM3nTyEE4bm89MXl+kqBRERCVMxiFNJiQn84eoJ5GWmcOsT89lb2+R1JBERiQIqBnGsICuV+6ZNpKy6ke88u4hAUOcbiIjEOxWDODdhYE/uuuho3l+zW7MiioiIioHA1MlFXFkygHvfXcfry3d6HUdERDykYiCYGb+4eAzjBuRw21MLmLF4u9eRRETEIyoGAkBaciKP33gcE4p68q2nF/LwrI1eRxIREQ+oGEhYTkYyj900mXOO7sMvXlnBf/1zJUGdkCgiEldUDORT0pITuW/asVxz/EAeeG8D//bXxTS16IZLIiLxIsnrABJ9EhNC0yb37ZHG/7yxhvKaRv50zbFkpWpzERHxO/2kl3aZGbedOYxe2an8+IVlTJ3+CQ9fP4le2aleR4s7gaCjtqmFmoYWahpDH80tQQLOEQxCwDkCwSCBYOhrg87REnQEg45A0LW+HlqelGAkJiSQlGAkJVr485SkBNKTE0lL3v+YSHpK62NyIokJ5vVqEJFuomIgnfrypIEUZKXyzacWcPn9H/HYjZMZlJ/pdayYEwg69tQ2Ubavgd37Ginb10h5TSPV9S3UNv5rh1/T0BIuAfsaQ6/VNQW8jk9qUgLZacn0SE+iR1oyPdKT6ZGWRI/0ZLLTQsvyMlMoyEolPyuFgszQY6ZGmURijvnx7nolJSVu3rx5XsfwlQVb9nLTI3NJMOOOc0dy2cQB+i0SqG8KtO7o/7XDb+/zitqmdmeWTE40slKTyEpLIjMliey0JDJTk8hKbX2eEnotq3XZ/tdSkhJIMCMxwUhMgMSEBBLNSEiApIQEEhNo83roI8EsNIIQdDQHggSCoZGFQNDR2BKgoTlIfVOA+uYADa0foedBahtbqG5oobqhmer6ZqobWtjX0Ex1fQvV9c00Bdo/DyU9OZH8rBTys1LpnZ1Kv5w0CnPTKcxJo39uOoW56fTJTiUpUac7iRwpM5vvnCs54vdRMZCDtX53Dd97dhFLSqsY1juLH5wzgrNH98HMvwWhoTlA6d46tu6tp3RP6HHrnjpK99azdW8dlXXNn/meBAtNN927Ryq9slLpnZ1Gr+zQ572zU0PPs9MoyEolPSXRg79V5DU0B9hT20R5TSMVNaHH8pomKmpCpai8ppFd1Q3sqGxgX2PLp743waB3dhqFuaGyMDg/k8EFmQzOz2BwQSb5mSm+3sZEIkXFoBMqBl3HOcdry3bymzdWs2F3LeOLcrljykhOOCrf62iHraaxhU3ltWyuqGNTRS2bymtDjxV17N7X+KmvTUlKYEDPdAb0zKCoZzr9ctM/tbPvlZ1KXmaKRlM6Ud3QzI7KBrZX1bOjsoEdVfVsb33cVllP6d76T42uZKcmhYpCQSZDWstCce8shvXO9k2xEokEFYNOqBh0vZZAkOfnl/K7t9ays7qBU4f34vazRzB2QI7X0dq1f+e/qSJUADaW17K5opaN5XWU13x65987OzX8G+vAvIxQCchLp6hnBgVZqSRop9+lmgNBtu6pC/87baqoDT9u21vP/s5gBoPzMxnRJ5sRfbMZ2Tf0OCg/U8VM4pKKQSdUDLpPQ3OAxz7exB/fXU9VfTP9c9OZPCSPyUPymDQ4j6N6ZXbbMHDbnX/osS78eODOv0+PVAblZzIkP5NBBRmhx/xMBhdkkJGiE+aiVVNLkC176lhXto9VO/exemfocVNFLft/lKUlJzCsdzajC3swfmAu44tyGd4nW2VBfE/FoBMqBt2vqr6ZFxduY/bGCuZs3EN5TRMA+Zkp4ZIweUgeowp7HPYP6OZAkL11TZRVN7az868N/5n7td35tz1mPShfO3+/qW8KsLZNWVi9cx/LtleFzwHJTElk7IAcxhf1ZMLAXCYU5dK7R5rHqUUiS8WgEyoG3nLOsbG8ljkb94Q+Nu2hdG89EDpeXNwni7Sk0DXzqQc+JicSDDr21DWxt7aJvXXN7K1rYk9tE/saWj7zZ/XpkRo6WU07fzmAc45NFXUs2rqXRVsqWbi1khXbq2lpPRbRLyeNCQN7Mr4ol/EDcxnbP4e0ZJ2zILFLxaATKgbRZ3tlPXM37WH2xj1srqilsTlIY0swfJlcY0uAxpYgDc0BEszomZFCXmYKuRmh6+N7ZqS0LkumV3ZoJEA7fzlUDc0Blm+vYuGWShZtrWThlkq2VYZKa0pSAiWDenLysAJOLi7g6H45OvwgMUXFoBMqBiJysMr2NbBoSyWzN+7hw3XlrNq5D4DcjGROPCqfk4t7ccqwAoryMjxOKtK5SBUD/bolInGtd3YaZx/dl7OP7guEisJH6yqYta6cWWvLeXXpTgAG5mWERxNOPCqf3IwUL2OLdBmNGIiIdMA5x/rdtXy4rpwP1pbzyYYKahpbMIMJRbmcO6aQKWP6ajRBooIOJXRCxUBEukJLIMji0ko+WFvOWyt3sWxbNQBj+vfg3DGFnDumL0N7ZXmcUuKVikEnVAxEpDtsqajjteU7eHXpThZtrQRgRJ9szh3bl3PHFDK8T5amc5Zuo2LQCRUDEelu2yvreW3ZTl5btpO5m/fgHAwtyGTKmL6cN7aQo/v1UEmQLqVi0AkVAxHxUtm+Bl5fvovXlu3gkw17CAQdowp7MHVyERcf05+cjGSvI4oPqRh0QsVARKLFntom/rFkO8/M3cry7dWkJiVw3thCrppUxOQheRpFkIhRMeiEioGIRKNl26p4Zu4WXlq4nX2NLQwpyOTLk4q4bOIAemWneh1PYlzcFAMzGwr8BMhxzl1+MN+jYiAi0ayuqYVXl+7k2blbmLtpL0kJxlmj+vDlyUWcOqyXZlyUwxLTxcDMHgYuAMqcc2PaLJ8C/B5IBB50zv13m9eeVzEQEb9ZV1bDs3O38LcF29hT20S/nDSuOWEQ0yYP0rkIckhivRicCtQAj+0vBmaWCKwBvgiUAnOBqc65Fa2vqxiIiG81tQR5a+Uunpy9mQ/XVZCRksiXJxVx40lDNIGSHJSYnhLZOfe+mQ0+YPFkYJ1zbgOAmT0DXAysOJj3NLNbgFsABg4cGLmwIiLdIKX1pMTzxhayYns1D36wgcc/3syjH23ivLGF3HLqUMYNyPU6psSBBK8DtNEf2Nrm81Kgv5nlm9n9wAQz+1FH3+ycm+6cK3HOlfTq1aurs4qIdJnR/Xrw2y+P54M7zuDmU4fy3urdXHTvh1z5wMe8t2Y30X5umMS2aLqJUntn2zjnXAVwa3eHERHxWmFOOj86dxS3nVHMs3O38vCsjVz38BwmDszlu2cN55RhBbrcUSIumkYMSoGiNp8PALZ7lEVEJGpkpyXz1VOGMvP2M/jVpWPZVd3IVx6ew+X3f8wHazWCIJEVTcVgLjDMzIaYWQpwFfCyx5lERKJGSlICVx83kHd+cBr/cckYtlfWc+1Dc7ji/o+ZtbZcBUEiwpNiYGZPAx8DI8ys1Mxucs61ALcBrwMrgeecc8u9yCciEs1SkxK55vhBzLz9dH55yRi2VdZzzUOzmfrnT1hSWul1PIlxUT/B0aEwswuBC4uLi29eu3at13FERLpFY0uAp2dv4Z531rGntomLjunH7eeM0GWOcSam5zHoaprHQETi0b6GZu5/bz0PfrAR5+C6Ewdx2xnDNFFSnIhUMYimcwxEROQIZKclc/s5I5l5++lcPL4fD87ayKm/eZeHZ22kORD0Op7ECBUDERGfKcxJ5zdXHMOr3z6FcQNy+MUrKzj/ng/4aH2519EkBqgYiIj41KjCHjx242SmX3ssdU0Brv7zbG57agE7quq9jiZRTMVARMTHzIyzj+7LW98/je+dNZw3V+zizP95jz++u47GloDX8SQKqRiIiMSBtOREvnPWMN76/mmcOryA37y+mim/0+EF+SxfFQMzu9DMpldVVXkdRUQkKhXlZfDAtSU8euNkgs5x9Z9n84O/LmZvbZPX0SRK6HJFEZE41dAc4J631zL9/Q30SE/mp+eP4tIJ/XX/hRilyxVFROSIpCUn8sMpI3nl2yczOD+D7z+3mGsems3milqvo4mHVAxEROLcyL49eP7WE/nlJWNYsrWKc373Pg/N2kgg6L8RZfl8KgYiIkJCgnHt8YN48/unceJRBfzylRVc+cDHrCur8TqadDMVAxERCeubk8ZD15Xwuy+PZ/3uGs675wP+NHM9LZo5MW6oGIiIyKeYGZdM6M8b3zuVM0f05u7XVnHZ/R+zfrdGD+KBr4qBLlcUEYmc3tlp3H/tsdx79QQ2V9Ry/j0f8MiHGwnq3ANf0+WKIiLyucqqG7jjb0t4d/VuTi4u4NeXj6NfbrrXsaQNXa4oIiLdpnePNB6+fhK/unQsC7bs5Zzfvc+LC7fhx18u452KgYiIHBQz4+rjBvLP75zCiD7ZfPfZRXzzqQXs0ayJvqJiICIih2RQfibPfu0E7pgykjdX7OKc373PO6t2eR1LIkTFQEREDlligvH104/ipW+eTH5mCjc+Mo8f/X0p9U26Y2OsUzEQEZHDNrpfD1667SS+dupQnp6zhYvuncXqnfu8jiVHQMVARESOSGpSIj86bxSP3TiZvXVNXHTvLJ6avUUnJsYoFQMREYmIU4f34tXvnMLkIXn8+IWl3PbUQqrqm72OJYfIV8VAExyJiHird3Yaj94wmTvPHcnry3dy/j0fsHDLXq9jySHwVTFwzs1wzt2Sk5PjdRQRkbiVkGDcetpRPHfrCQBc+cDHPDl7s8ep5GD5qhiIiEj0mDiwJ//41imcVFzAT15Yxo/+vpTGFl21EO1UDEREpMvkZCTz0HWT+MbpR/H0nC1Mnf4JZdUNXseSTqgYiIhIl0pMMH44ZSR/vHoiK3fs44I/zGLZNp0LFq1UDEREpFucP66QF755IsmJCVw1/RM+Xl/hdSRph4pBvGnRnOYi4p2RfXvw/NdPoDAnjev+Moe3Vmgq5WijYhBPGmvgsYtg5n97nURE4lhhTjrPfe0ERvXN5htPLdDIQZRRMYgXjTXw5BWwdQ4UDPc6jYjEuZ6ZKTx642QG5WXwtcfn8cHa3Zz9f++xZpemU/aaikE8aKiGJy6DrbPhsj/DmC95nUhEhNyMFB68roSWQJCbHp3H2rIabvjLXOqaWryOFtd8VQw082E76ivh8Uth2zy4/GEYc5nXiUREwgblZzKkIIumliDOQXlNIz98fonXseKar4qBZj48QG05PHoB7FgMVz4GR1/idSIRkU95bu5WNpTXhj9vbAny9soynpu71cNU8c1XxUDaqN4Bj5wP5Wvh6mdg5PleJxIR+Yy7X1tFffOnZ0Osbw5w92urPEokKgZ+tGcjPHwOVJXCtOeh+CyvE4mItOuOKSNJT0781LL05ETuPHekR4lExcBvylbBX86Fhir4yssw5BSvE4mIdOjKSUWcOao3qUmh3VFqUgJfGNWbK0qKPE4Wv1QM/GTb/FApcEG44VUYcKzXiUREPtdvLh9HQVYKBhRkpfLry8d5HSmuqRj4xcYP4NGLIDULbvgn9Dna60QiIgclIyWJv9wwmWF9svjLDZPISEnyOlJc09r3g1Wvwl+vh7whcO0L0KOf14lERA7J8D7ZvPG907yOIWjEIPYtfhaevQb6jgmNFKgUiIjIEVAxiGWzH4AXboHBJ8FXXoKMPK8TiYhIjNOhhFjkHLz3a5j5Kxh5AVz2ECSneZ1KRER8QMUg1gSD8MZP4JP74Jir4aI/QKL+GUVEJDK0R4klwQC8/G1Y9AQc93U451eQoKNBIiISOb7aq/j6JkotTfD8jaFScNqdMOW/VApERCTifLVn8e1NlJrr4dlpsOJFOPs/4YwfgZnXqURExId0KCHaNdXC01eFJjC64HdQcoPXiURExMdUDKJZ4z548grYOhsufQCO+bLXiURExOdUDKJVQzU8cVno/geXPQRjvuR1IhERiQO+Oscg5pWthD8eD6XzQqVg+wK44hGVAhER6TYaMYgWTbWhwwZVpfCX8yDQDFc+AqMv8jqZiIjEEY0YRIuXvgm1ZYCDQGPolsmjL/Y6lYiIxBkVg2iw4AlY8zq0NP5r2a7loeUiIiLdSMUgGrx9FzTXfXpZc11ouYiISDdSMYgGX7gLkjM+vSw5A876uSdxREQkfqkYRIOJ18DwcyCp9Q6JSWkwfApMmOZtLhERiTsqBtHi4j9CZi/AQo8X3+t1IhERiUMqBtEiJROm/RV6jQw9pmR6nUhEROKQ5jGIJr1HwTc/8TqFiIjEMY0YiIiISJivioGZXWhm06uqqryOIiIiEpN8VQycczOcc7fk5OR4HUVERCQm+aoYiIiIyJFRMRAREZEwFQMREREJUzEQERGRMBUDERERCVMxEBERkTAVAxEREQlTMRAREZEwFQMREREJUzEQERGRMBUDERERCVMxEBERkTAVAxEREQlTMRAREZEwFQMREREJUzEQERGRMBUDERERCVMxEBERkTAVAxEREQlTMRAREZEwFQMREREJ81UxMLMLzWx6VVWV11FERERikq+KgXNuhnPulpycHK+jiIiIxCRfFQMRERE5MioGIiIiEqZiICIiImEqBiIiIhKmYiAiIiJhKgYiIiISpmIgIiIiYSoGIiIiEqZiICIiImEqBiIiIhKmYiAiIiJhKgYiIiISpmIgIiIiYSoGIiIiEqZiICIiImEqBiIiIhKmYiAiIiJhKgYiIiISpmIgIiIiYSoGIiIiEqZiICIiImEqBiIiIhKmYiAiIiJhKgYiIiISpmIgIiIiYSoGIiIiEqZiICIiImEqBiIiIhKmYiAiIiJhKgYiIiISpmIgIiIiYSoGIiIiEqZiICIiImFJXgf4PGaWCdwHNAEznXNPehxJRETEtzwZMTCzh82szMyWHbB8ipmtNrN1ZnZn6+IvAc87524GLur2sCIiInHEq0MJjwBT2i4ws0Tgj8C5wGhgqpmNBgYAW1u/LNCNGUVEROKOJ8XAOfc+sOeAxZOBdc65Dc65JuAZ4GKglFA5AJ0TISIi0qWi6RyD/vxrZABCheA44B7gXjM7H5jR0Teb2S3ALa2fNh54mEIirgAo9zpEHNB67npax11P67h7jIjEm0RTMbB2ljnnXC1ww+d9s3NuOjAdwMzmOedKIpxP2tA67h5az11P67jraR13DzObF4n3iaah+VKgqM3nA4DtHmURERGJS9FUDOYCw8xsiJmlAFcBL3ucSUREJK54dbni08DHwAgzKzWzm5xzLcBtwOvASuA559zyw/wjpkcoqnRM67h7aD13Pa3jrqd13D0isp7NOReJ9xEREREfiKZDCSIiIuKxmCwG7c2caGZ5Zvamma1tfezZ5rUftc6muNrMzvEmdWzpYB3fZWbbzGxR68d5bV7TOj5EZlZkZu+a2UozW25m32ldrm05QjpZx9qWI8TM0sxsjpktbl3HP29dru04gjpZz5Hflp1zMfcBnApMBJa1WfZr4M7W53cCd7c+Hw0sBlKBIcB6INHrv0O0f3Swju8CftDO12odH946LgQmtj7PBta0rktty12/jrUtR24dG5DV+jwZmA0cr+2429ZzxLflmBwxcO3PnHgx8Gjr80eBS9osf8Y51+ic2wisIzTLonSig3XcEa3jw+Cc2+GcW9D6fB+hk277o205YjpZxx3ROj5ELqSm9dPk1g+HtuOI6mQ9d+Sw13NMFoMO9HHO7YDQDwOgd+vy9mZU7OwHg3TuNjNb0nqoYf/QoNbxETKzwcAEQr8FaFvuAgesY9C2HDFmlmhmi4Ay4E3nnLbjLtDBeoYIb8t+KgYdaXdGxW5P4Q9/Ao4CxgM7gP9tXa51fATMLAv4G/Bd51x1Z1/azjKt54PQzjrWthxBzrmAc248oYnpJpvZmE6+XOv4MHWwniO+LfupGOwys0KA1sey1uWaUTFCnHO7WjfMIPBn/jUspXV8mMwsmdAO60nn3N9bF2tbjqD21rG25a7hnKsEZhK6e6624y7Sdj13xbbsp2LwMnBd6/PrgJfaLL/KzFLNbAgwDJjjQb6Yt/8/eatLgf1XLGgdHwYzM+AhYKVz7rdtXtK2HCEdrWNty5FjZr3MLLf1eTpwFrAKbccR1dF67optOZpuonTQLDRz4ulAgZmVAv8O/DfwnJndBGwBrgBwzi03s+eAFUAL8E3nXMCT4DGkg3V8upmNJzQctQn4GmgdH4GTgGuBpa3HDQF+jLblSOpoHU/VthwxhcCjZpZI6JfN55xzr5jZx2g7jqSO1vPjkd6WNfOhiIiIhPnpUIKIiIgcIRUDERERCVMxEBERkTAVAxEREQlTMRAREZEwFQMREREJUzEQERGRMBUDEemUmd3c5l7vwTbPf9vO1z5gZicdsKymzfPzzGytmQ3sjuwicug0wZGIHBQz6w985Jwb1MnXLAKObTvDmpnVOOeyzOwLwHTgbOfc+i4PLCKHJSanRBYRT4wBlnb0opmNAta0N+2qmZ1C6AYv56kUiEQ3FQMROVhj+dcNWtpzLvBaO8tTCd1A53Tn3KquCCYikaNzDETkYHU6YgCcQ/vFoBn4CLipK0KJSGSpGIjIwepwxMDMMoBc51x793sPAlcCk8zsx12YT0QiQIcSRORzmVkCofu5d3Qo4Azg3Y6+3zlXZ2YXAB+Y2S7n3ENdEFNEIqDTEQMzu771PzNm9gszSz+SP8zM8szskdb7R7d3qdPtZnaPmX3vSP4cEYm4YqDUOdfYwesdnV8Q5pzbA0wBfmpmF0c4n0hEtN3vRfA9/83M7m29nNfaLO9jZs+Z2X1mdms731doZhvMbEwH73u6md12JBna83kjBicDGa3vMRBINLNHgC3AKOD91uVJzrnvmdnXgeFAT+CnwJmEfpi8A+EfDNe3hvyrmSU454Ktn08ATgJWAzsO/AsBg4Aq59zPDmYliEjkOOfWAKM7+ZITgXYLvXMuq83zrcCQyKYTiai2+71G4HwgHfgboX3fXcAaYLJzboqZrQAeIHSo7TvAZbTZ75lZCjDROTetdSd+MvBBmz9rhnPucTN73swecs41t8lyO/DXAwO2/vI8CMgB5pvZIODfAAPWA3uBPc65GWb2l9ZcHWX4jM87x2AW8JRz7pUDlt8P/AcwxDl3O1BkZlnAV4Cq1lATnXOP7V85B/ylTgFW7S8FrUYAK51zdwDnHzA60ReYB9zzOXlFxAPOuYkH/EATiVVt93vfBioJ/bI6GbgZuAP4BZDc+vXbnXO/B/4BXNTOfi8f2N36fDMwoM1rrwITzex/Cf1Cnb//BTO7gVAZqW8n46nOue/yr1G6b7R+XQWhgvI34FIzywZagMxOMnzG540YBDtYXg30aH0M/z2Abc65uzp7QzM7HbgQ+MEBL5USGn0AqCN0idP+FXIHMAn4i5ld7ZyrRkREJPLa7vcSgP9wzrUAtB4Cd60f++3fjyYfsHy/CqCg9flAYMn+F5xz9bSOtJnZS0BZm++bDBwDHE+oMHyrzWtNrY/7D+0lAI8758LvbWYOuA74e2cZ2vN5xWAx8BMz+9yTFJ1z+8xsjpn9gVBJeBgYB2x1zr3dGrQP8BzwAvCn1uGQC1r/cjOAqa0rfqdzrrLN2/+w9S+1h1BpEBER6Qpt93v3AA+a2R5Co9Z/Bu4mNFy/f4Qs38x+RegQ2VfN7Hra7Pecc01mtsDMfk/oF977zOwKQvu9t4E/AInAo865YOs+8L+cc18HMLO7gOdbnz/unLsW+NDMfgQcBSwC7gV+ZWY7gH3OuZ8TGjW4FxjunGs5MENnK0BTIouIiBwiM3veOXf5/kev80SSioGIiIiEaYIjERERCVMxEBERkTAVAxEREQlTMRAREZEwFQMREREJUzEQERGRsP8Pk2wM3szc5mIAAAAASUVORK5CYII=\n",
+ "text/plain": [
+ "<Figure size 504x432 with 1 Axes>"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "import timeit\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt \n",
+ "import pandas\n",
+ "import teqp\n",
+ "\n",
+ "def get_critical_curve(ipure):\n",
+ " \"\"\" Return curve as pandas DataFrame \"\"\"\n",
+ " names = ['Nitrogen', 'Ethane']\n",
+ " model = teqp.build_multifluid_model(names, teqp.get_datapath())\n",
+ " T0 = model.get_Tcvec()[ipure]\n",
+ " rho0 = np.array([1.0/model.get_vcvec()[ipure]]*2)\n",
+ " rho0[1-ipure] = 0\n",
+ " o = teqp.TCABOptions() \n",
+ " o.init_dt = 1.0 # step in the parameter\n",
+ " o.rel_err = 1e-8\n",
+ " o.abs_err = 1e-5\n",
+ " o.integration_order = 5\n",
+ " o.calc_stability = True\n",
+ " o.polish = True\n",
+ " curveJSON = teqp.trace_critical_arclength_binary(model, T0, rho0, '', o)\n",
+ " df = pandas.DataFrame(curveJSON)\n",
+ " rhotot = df['rho0 / mol/m^3']+df['rho1 / mol/m^3']\n",
+ " df['z0 / mole frac.'] = df['rho0 / mol/m^3']/rhotot\n",
+ " return df\n",
+ "\n",
+ "fig, ax = plt.subplots(1,1,figsize=(7, 6))\n",
+ "tic = timeit.default_timer()\n",
+ "for ipure in [1,0]:\n",
+ " df = get_critical_curve(ipure)\n",
+ " first_unstable = np.argmax(~df['locally stable'])\n",
+ " df = df.iloc[0:(first_unstable if first_unstable else len(df))]\n",
+ " line, = plt.plot(df['T / K'], df['p / Pa']/1e6, '-')\n",
+ " plt.plot(df['T / K'].iloc[0], df['p / Pa'].iloc[0]/1e6, 'd', \n",
+ " color=line.get_color())\n",
+ "\n",
+ "elap = timeit.default_timer()-tic\n",
+ "plt.gca().set(xlabel='$T$ / K', ylabel='$p$ / MPa',\n",
+ " xlim=(100, 350), ylim=(1, 1e3))\n",
+ "plt.yscale('log')\n",
+ "plt.tight_layout(pad=0.2)\n",
+ "plt.gcf().text(0,0,f'time: {elap:0.1f} s', ha='left', va='bottom', fontsize=7)\n",
+ "plt.gcf().text(1,0,f'teqp: {teqp.__version__}', ha='right', va='bottom', fontsize=7);"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "6120531f",
+ "metadata": {},
+ "source": [
+ "And now for something a bit more interesting: ethane + alkane critical curves"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "92020440",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import timeit\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt \n",
+ "import pandas\n",
+ "import teqp\n",
+ "\n",
+ "def get_critical_curve(names, ipure):\n",
+ " \"\"\" Return curve as pandas DataFrame \"\"\"\n",
+ " model = teqp.build_multifluid_model(names, teqp.get_datapath())\n",
+ " T0 = model.get_Tcvec()[ipure]\n",
+ " rho0 = np.array([1.0/model.get_vcvec()[ipure]]*2)\n",
+ " rho0[1-ipure] = 0\n",
+ " o = teqp.TCABOptions() \n",
+ "# print(dir(o))\n",
+ " o.init_dt = 1.0 # step in the parameter\n",
+ " o.rel_err = 1e-6 # relative error on the step\n",
+ " o.abs_err = 1e-6 # absolute error on the step\n",
+ " o.max_dt = 100 # cap the size of the allowed step\n",
+ " o.calc_stability = True\n",
+ " o.polish = True\n",
+ " curveJSON = teqp.trace_critical_arclength_binary(model, T0, rho0, '', o)\n",
+ " df = pandas.DataFrame(curveJSON)\n",
+ " rhotot = df['rho0 / mol/m^3']+df['rho1 / mol/m^3']\n",
+ " df['z0 / mole frac.'] = df['rho0 / mol/m^3']/rhotot\n",
+ " return df\n",
+ "\n",
+ "fig, ax = plt.subplots(1,1,figsize=(7, 6))\n",
+ "tic = timeit.default_timer()\n",
+ "name0 = 'ETHANE'\n",
+ "for othername in ['METHANE','PROPANE','BUTANE','PENTANE','HEXANE']:\n",
+ " for ipure in [1]:\n",
+ " df = get_critical_curve([name0, othername], ipure)\n",
+ " line, = plt.plot(df['T / K'], df['p / Pa']/1e6, '-')\n",
+ " plt.plot(df['T / K'].iloc[0], df['p / Pa'].iloc[0]/1e6, 'd', \n",
+ " color=line.get_color())\n",
+ "\n",
+ "elap = timeit.default_timer()-tic\n",
+ "plt.gca().set(xlabel='$T$ / K', ylabel='$p$ / MPa')#,xlim=(100, 350), ylim=(1, 1e3))\n",
+ "plt.tight_layout(pad=0.2)\n",
+ "plt.gcf().text(0,0,f'time: {elap:0.1f} s', ha='left', va='bottom', fontsize=7)\n",
+ "plt.gcf().text(1,0,f'teqp: {teqp.__version__}', ha='right', va='bottom', fontsize=7);"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "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.10.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/doc/source/algorithms/index.rst b/doc/source/algorithms/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..13225b9ff81375465bb8140e44e6e13a58b07862
--- /dev/null
+++ b/doc/source/algorithms/index.rst
@@ -0,0 +1,14 @@
+Algorithms
+==========
+
+.. toctree::
+ :maxdepth: 2
+ :caption: Contents:
+
+ VLE
+ critical_curves
+
+Information
+-----------
+
+The algorithms are written in a very generic way; they take an instance of a thermodynamic model, and the necessary derivatives are calculated from this model with automatic differentiation (or similar). In that way, implementing a model is all that is required to enable its use in the calculation of critical curves or to trace the phase equilibria. Determining the starting values, on the other hand, may require model-specific assistance, for instance with superancillary equations.
\ No newline at end of file
diff --git a/doc/source/conf.py b/doc/source/conf.py
new file mode 100644
index 0000000000000000000000000000000000000000..c0474e69ce2e88279148890d4bd1a60e5e8eef70
--- /dev/null
+++ b/doc/source/conf.py
@@ -0,0 +1,75 @@
+# Configuration file for the Sphinx documentation builder.
+#
+# This file only contains a selection of the most common options. For a full
+# list see the documentation:
+# https://www.sphinx-doc.org/en/master/usage/configuration.html
+
+import os, subprocess
+on_rtd = os.environ.get('READTHEDOCS') == 'True'
+
+# -- Path setup --------------------------------------------------------------
+
+# If extensions (or modules to document with autodoc) are in another directory,
+# add these directories to sys.path here. If the directory is relative to the
+# documentation root, use os.path.abspath to make it absolute, like shown here.
+#
+# import os
+# import sys
+# sys.path.insert(0, os.path.abspath('.'))
+
+
+# -- Project information -----------------------------------------------------
+
+project = 'teqp'
+copyright = '2022, Ian Bell'
+author = 'Ian Bell'
+
+# The full version, including alpha/beta/rc tags
+import teqp
+release = teqp.__version__
+
+# -- Execute all notebooks --------------------------------------------------
+if on_rtd:
+ # subprocess.check_output(f'jupyter nbconvert --version', shell=True)
+ for path, dirs, files in os.walk('.'):
+ for file in files:
+ if file.endswith('.ipynb') and '.ipynb_checkpoints' not in path:
+ subprocess.check_output(f'jupyter nbconvert --to notebook --output {file} --execute {file}', shell=True, cwd=path)
+ # --ExecutePreprocessor.allow_errors=True (this allows you to allow errors globally, but a raises-exception cell tag is better)
+
+### -- Auto-generate API documentation -----------------------------------------
+here = os.path.dirname(__file__)
+subprocess.check_output(f'sphinx-apidoc -f -o api {os.path.dirname(teqp.__file__)}', shell=True, cwd=here)
+
+# -- General configuration ---------------------------------------------------
+
+# Add any Sphinx extension module names here, as strings. They can be
+# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
+# ones.
+extensions = [
+'nbsphinx',
+'sphinx.ext.autodoc'
+]
+
+# Add any paths that contain templates here, relative to this directory.
+templates_path = ['_templates']
+
+# List of patterns, relative to source directory, that match files and
+# directories to ignore when looking for source files.
+# This pattern also affects html_static_path and html_extra_path.
+exclude_patterns = []
+
+# -- Options for HTML output -------------------------------------------------
+
+# The theme to use for HTML and HTML Help pages. See the documentation for
+# a list of builtin themes.
+#
+if on_rtd:
+ html_theme = 'default'
+else:
+ html_theme = 'insipid'
+
+# Add any paths that contain custom static files (such as style sheets) here,
+# relative to this directory. They are copied after the builtin static files,
+# so a file named "default.css" will overwrite the builtin "default.css".
+html_static_path = ['_static']
\ No newline at end of file
diff --git a/doc/source/derivs/derivs.ipynb b/doc/source/derivs/derivs.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..a2015a904564a8e106710bef0bb1e54d56aba7e0
--- /dev/null
+++ b/doc/source/derivs/derivs.ipynb
@@ -0,0 +1,332 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "1afc4517",
+ "metadata": {},
+ "source": [
+ "# Thermodynamic Derivatives\n",
+ "\n",
+ "## Helmholtz energy derivatives\n",
+ "Thermodynamic derivatives are at the very heart of teqp. All models are defined in the form $\\alpha^r(T, \\rho, z)$, where $\\rho$ is the molar density, and z are mole fractions. There are exceptions for models for which the independent variables are in simulation units (Lennard-Jones and its ilk).\n",
+ "\n",
+ "Thereofore, to obtain the residual pressure, it is obtained as a derivative:\n",
+ "\n",
+ "$$ p^r = \\rho R T \\left( \\rho \\left(\\frac{\\partial \\alpha^r}{\\partial \\rho}\\right)_{T}\\right)$$\n",
+ "\n",
+ "and other residual thermodynamic properties are defined likewise.\n",
+ "\n",
+ "We can define the concise derivative\n",
+ "\n",
+ "$$ \\Lambda^r_{xy} = (1/T)^x(\\rho)^y\\left(\\frac{\\partial^{x+y}(\\alpha^r)}{\\partial (1/T)^x\\partial \\rho^y}\\right)$$\n",
+ "\n",
+ "so we can re-write the derivative above as \n",
+ "\n",
+ "$$ p^r = \\rho R T \\Lambda^r_{01}$$\n",
+ "\n",
+ "Similar definitions apply for all the other thermodynamic properties, with the tot superscript indicating it is the sum of the residual and ideal-gas (not included in teqp) contributions:\n",
+ "\n",
+ "$$\n",
+ "\\frac{p}{\\rho R T} = 1 + \\Lambda^r_{01}\n",
+ "$$\n",
+ "Internal energy ($u= a+Ts$):\n",
+ "$$\n",
+ "\t\\frac{u}{RT} = \\Lambda^{\\rm tot}_{10}\n",
+ "$$\n",
+ "Enthalpy ($h= u+p/\\rho$):\n",
+ "$$\n",
+ "\\frac{h}{RT} = 1+\\Lambda^r_{01} + \\Lambda^{\\rm tot}_{10}\n",
+ "$$\n",
+ "Entropy ($s\\equiv -(\\partial a/\\partial T)_v$):\n",
+ "$$\n",
+ "\\frac{s}{R} = \\Lambda^{\\rm tot}_{10}-\\Lambda^{\\rm tot}_{00}\n",
+ "$$\n",
+ "Gibbs energy ($g= h-Ts$):\n",
+ "$$\n",
+ "\t\\frac{g}{RT} = 1+\\Lambda^r_{01}+\\Lambda^{\\rm tot}_{00}\n",
+ "$$\n",
+ "Derivatives of pressure:\n",
+ "$$\n",
+ "\t\\left(\\frac{\\partial p}{\\partial \\rho}\\right)_{T} = RT\\left(1+2\\Lambda^r_{01}+\\Lambda^r_{02}\\right)\n",
+ "$$\n",
+ "$$\n",
+ "\t\\left(\\frac{\\partial p}{\\partial T}\\right)_{\\rho} = R\\rho\\left(1+\\Lambda^r_{01}-\\Lambda^r_{11}\\right)\n",
+ "$$\n",
+ "Isochoric specific heat ($c_v\\equiv (\\partial u/\\partial T)_v$):\n",
+ "$$\n",
+ "\\frac{c_v}{R} = -\\Lambda^{\\rm tot}_{20}\n",
+ "$$\n",
+ "Isobaric specific heat ($c_p\\equiv (\\partial h/\\partial T)_p$; see Eq. 3.56 from Span \\cite{Span-BOOK-2000} for the derivation):\n",
+ "$$\n",
+ "\\frac{c_p}{R} = -\\Lambda^{\\rm tot}_{20}+\\frac{(1+\\Lambda^r_{01}-\\Lambda^r_{11})^2}{1+2\\Lambda^r_{01}+\\Lambda^r_{02}}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "ca35cc29",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2022-07-06T18:39:50.124207Z",
+ "iopub.status.busy": "2022-07-06T18:39:50.124207Z",
+ "iopub.status.idle": "2022-07-06T18:39:50.261696Z",
+ "shell.execute_reply": "2022-07-06T18:39:50.261181Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "import teqp\n",
+ "import numpy as np"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "8d254c40",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2022-07-06T18:39:50.264736Z",
+ "iopub.status.busy": "2022-07-06T18:39:50.264736Z",
+ "iopub.status.idle": "2022-07-06T18:39:50.277541Z",
+ "shell.execute_reply": "2022-07-06T18:39:50.276780Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "Tc_K = [300]\n",
+ "pc_Pa = [4e6]\n",
+ "acentric = [0.01]\n",
+ "model = teqp.canonical_PR(Tc_K, pc_Pa, acentric)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "49ee3453",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2022-07-06T18:39:50.280537Z",
+ "iopub.status.busy": "2022-07-06T18:39:50.280537Z",
+ "iopub.status.idle": "2022-07-06T18:39:50.292555Z",
+ "shell.execute_reply": "2022-07-06T18:39:50.292555Z"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "-0.06836660379313926"
+ ]
+ },
+ "execution_count": 3,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "z = np.array([1.0])\n",
+ "model.get_Ar01(300,300,z)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "498a2350",
+ "metadata": {},
+ "source": [
+ "And there are additional methods to obtain all the derivatives up to a given order:\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "3ae755e1",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "[-0.06966138343515413,\n",
+ " -0.06836660379313926,\n",
+ " 0.0025357822532378147,\n",
+ " -0.00015701162203571184,\n",
+ " 1.6818628788290574e-05,\n",
+ " -2.2305940927885907e-06,\n",
+ " 3.8259258513417917e-07]"
+ ]
+ },
+ "execution_count": 4,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "model.get_Ar06n(300,300,z) # derivatives 00, 01, 02, ... 06"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "cd644a3e",
+ "metadata": {},
+ "source": [
+ "But more derivatives are slower than fewer:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "4e0609ae",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "1.9 µs ± 38.5 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)\n",
+ "2.07 µs ± 24.8 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)\n"
+ ]
+ }
+ ],
+ "source": [
+ "%timeit model.get_Ar01(300,300,z)\n",
+ "%timeit model.get_Ar04n(300,300,z)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d05e9070",
+ "metadata": {},
+ "source": [
+ "Note: calling overhead is usually on the order of 1 microsecond"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "502e9830",
+ "metadata": {},
+ "source": [
+ "## Virial coefficients\n",
+ "\n",
+ "Virial coefficients represent the thermodynamics of the interaction of two-, three-, ... bodies interacting with each other. They can be obtained rigorously if the potential energy surface of interaction is fully known. In general, such a surface can only be constructed for small rigid molecules. Many simple thermodynamic models do a poor job of predicting the thermodynamics captured by the virial coefficients. \n",
+ "\n",
+ "The i-th virial coefficient is defined by \n",
+ "$$\n",
+ "\tB_i = \\frac{(\\alpha^r)^{(i-1)}}{(i-2)!}\n",
+ "$$\n",
+ "with the concise derivative term\n",
+ "$$\n",
+ "\t(\\alpha^r)^{(i)} = \\lim_{\\rho \\to 0} \\left(\\frac{\\partial^{i}\\alpha^r}{\\partial \\rho^{i}}\\right)_{T,\\vec{x}}\n",
+ "$$\n",
+ "\n",
+ "teqp supports the virial coefficient directly, there is the ``get_B2vir`` method for the second virial coefficient:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "id": "720e120c",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "-0.00023661263734465424"
+ ]
+ },
+ "execution_count": 6,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "model.get_B2vir(300, z)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "2fba4369",
+ "metadata": {},
+ "source": [
+ "And the ``get_Bnvir`` method that allows for the calculation of higher virial coefficients:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "92d01992",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "{2: -0.0002366126373446542,\n",
+ " 3: 3.001768410777936e-08,\n",
+ " 4: -3.2409760373816355e-12,\n",
+ " 5: 3.9617816466337214e-16,\n",
+ " 6: -4.552923983836698e-20,\n",
+ " 7: 5.3759278511184914e-24}"
+ ]
+ },
+ "execution_count": 7,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "model.get_Bnvir(7, 300, z)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d50998b1",
+ "metadata": {},
+ "source": [
+ "The ``get_Bnvir`` method was implemented because when doing automatic differentiation, all the intermediate derivatives are also retained.\n",
+ "\n",
+ "There is also a method to calculate temperature derivatives of a given virial coefficient"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "id": "fa4ff386",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "1.0095625628421257e-10"
+ ]
+ },
+ "execution_count": 8,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "model.get_dmBnvirdTm(2, 3, 300, z) # third temperature derivative of the second virial coefficient"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "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.10.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/doc/source/derivs/index.rst b/doc/source/derivs/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..8dccdb9bef2cb5289f2817e4eb7020b1b9d4fa0d
--- /dev/null
+++ b/doc/source/derivs/index.rst
@@ -0,0 +1,8 @@
+Derivatives
+===========
+
+.. toctree::
+ :maxdepth: 2
+ :caption: Contents:
+
+ derivs
\ No newline at end of file
diff --git a/doc/source/getting_started/index.rst b/doc/source/getting_started/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..cfba90910fb921350b56007e0aa558091097fb04
--- /dev/null
+++ b/doc/source/getting_started/index.rst
@@ -0,0 +1,39 @@
+Getting Started
+===============
+
+Introduction
+------------
+
+teqp (phonetically: tɛk pi) is a C++-based library with wrappers. It was written because implementing EOS (particularly working out the derivatives) is a painful, error-prone, and slow process. The advent of open-source automatic differentiation libraries makes the implementation of EOS as fast as hand-written derivatives, and much easier to implement without errors.
+
+There is a paper about teqp: https://doi.org/10.1021/acs.iecr.2c00237
+
+The documentation is based on the Python wrapper because it can be readily integrated with the documentation tools (sphinx in this case) and can be auto-generated at documentation build time.
+
+Installation
+------------
+
+Python
+^^^^^^
+
+The library can be installed with:
+
+.. code::
+
+ pip install teqp
+
+because the binary wheels for all major platforms are provided on pypi.
+
+If you desire to build teqp yourself, it is recommended to pull from github and build a binary wheel, and then subsequently install that wheel:
+
+.. code::
+
+ git clone --recursive https://github.com/usnistgov/teqp
+ cd teqp
+ python setup.py bdist_wheel
+ pip install dist/*.whl # or replace with the appropriate binary wheel
+
+C++
+^^^
+
+The build is cmake based. There are targets available for an interface library, etc. Please see ``CMakeLists.txt``
\ No newline at end of file
diff --git a/doc/source/index.rst b/doc/source/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..c537fabc96e31273233d4b77b45e43a965384da3
--- /dev/null
+++ b/doc/source/index.rst
@@ -0,0 +1,19 @@
+Welcome to teqp's documentation!
+================================
+
+.. toctree::
+ :maxdepth: 2
+ :caption: Contents:
+
+ getting_started/index
+ models/index
+ derivs/index
+ algorithms/index
+ api/modules
+
+Indices and tables
+==================
+
+* :ref:`genindex`
+* :ref:`modindex`
+* :ref:`search`
diff --git a/doc/source/models/CPA.ipynb b/doc/source/models/CPA.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..47e2d9361da8a696fa3208cc44bb1d82fc31a753
--- /dev/null
+++ b/doc/source/models/CPA.ipynb
@@ -0,0 +1,35 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "df447d16",
+ "metadata": {},
+ "source": [
+ "# Cubic Plus Association (CPA)\n",
+ "\n",
+ "The combination of a cubic EOS with association"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "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.10.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/doc/source/models/PCSAFT.ipynb b/doc/source/models/PCSAFT.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..973f43bd867d0fefdfb2d68d33553a5f96c9a32a
--- /dev/null
+++ b/doc/source/models/PCSAFT.ipynb
@@ -0,0 +1,155 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "f406bbb5",
+ "metadata": {},
+ "source": [
+ "# PC-SAFT\n",
+ "\n",
+ "The PC-SAFT implementation in teqp is based on the implementation of Gross and Sadowski (https://doi.org/10.1021/ie0003887), with the typo from their paper fixed. It does NOT include the association contribution, only the dispersive contributions.\n",
+ "\n",
+ "The model in teqp requires the user to specify the values of ``sigma``, ``epsilon/kB``, and ``m`` for each substance. A very few substances are hardcoded in teqp, for testing purposes. "
+ ]
+ },
+ {
+ "cell_type": "raw",
+ "id": "d9efd027",
+ "metadata": {
+ "raw_mimetype": "text/restructuredtext"
+ },
+ "source": [
+ "The Python class is here: :py:class:`PCSAFTEOS <teqp.teqp.PCSAFTEOS>`"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "984925ce",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import teqp\n",
+ "import numpy as np\n",
+ "teqp.__version__"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "7bbd7129",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "TeXkey = 'Gross-IECR-2001'\n",
+ "ms = [1.0, 1.6069, 2.0020]\n",
+ "eoverk = [150.03, 191.42, 208.11]\n",
+ "sigmas = [3.7039, 3.5206, 3.6184]\n",
+ "\n",
+ "coeffs = []\n",
+ "for i in range(len(ms)):\n",
+ " c = teqp.SAFTCoeffs()\n",
+ " c.m = ms[i]\n",
+ " c.epsilon_over_k = eoverk[i]\n",
+ " c.sigma_Angstrom = sigmas[i]\n",
+ " coeffs.append(c)\n",
+ " \n",
+ "model = teqp.PCSAFTEOS(coeffs)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "578630c8",
+ "metadata": {},
+ "source": [
+ "The model parameters can be queried:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d4e47e54",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "model.get_m(), model.get_epsilon_over_k_K(), model.get_sigma_Angstrom()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "5cbb382d",
+ "metadata": {},
+ "source": [
+ "## Adjusting k_ij\n",
+ "\n",
+ "Fine-tuned values of $k_{ij}$ can be provided when instantiating the model. A complete matrix of all the $k_{ij}$ values must be provided. This allows for asymmetric mixing models in which $k_{ij}\\neq k_{ji}$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a32c41b5",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "k_01 = 0.01; k_10 = k_01\n",
+ "kmat = [[0,k_01,0],[k_10,0,0],[0,0,0]]\n",
+ "teqp.PCSAFTEOS(coeffs, kmat)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ca52e844",
+ "metadata": {},
+ "source": [
+ "## Superancillary\n",
+ "\n",
+ "The superancillary equation for PC-SAFT has been developed, and is much more involved than that of the cubic EOS. As a consequence, the superancillary equation has been provided as a separate package rather than integrating it into to teqp to minimize the binary size of teqp. It can be installed from PYPI with: ``pip install PCSAFTsuperanc``"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "0bdf568f",
+ "metadata": {},
+ "source": [
+ "## Maximum density\n",
+ "\n",
+ "The maximum number density allowed by the EOS is defined based on the packing fraction. To get a molar density, divide by Avogadro's number. The function is conveniently exposed in Python:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "3c8491a9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "max_rhoN = teqp.PCSAFTEOS(coeffs).max_rhoN(130.0, np.array([0.3, 0.3, 0.4]))\n",
+ "display(max_rhoN)\n",
+ "max_rhoN/6.022e23 # the maximum molar density in mol/m^3"
+ ]
+ }
+ ],
+ "metadata": {
+ "celltoolbar": "Raw Cell Format",
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "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.10.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/doc/source/models/cubics.ipynb b/doc/source/models/cubics.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..d9c6ee29c6a4cf06207b95ee529de2ae3460f449
--- /dev/null
+++ b/doc/source/models/cubics.ipynb
@@ -0,0 +1,150 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "df447d16",
+ "metadata": {},
+ "source": [
+ "# General cubics\n",
+ "\n",
+ "The reduced residual Helmholtz energy for the main cubic EOS (van der Waals, Peng-Robinson, and Soave-Redlich-Kwong) can be written in a common form (see https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7365965/)\n",
+ "\n",
+ "$$ \\alpha^r = \\psi^{(-)} - \\dfrac{\\tau a_m}{RT_r } \\psi^{(+)} $$\n",
+ "\n",
+ "$$ \\psi^{(-)} =-\\ln(1-b_m\\rho ) $$\n",
+ "\n",
+ "$$ \\psi^{(+)} = \\dfrac{\\ln\\left(\\dfrac{\\Delta_1 b_m\\rho+1}{\\Delta_2b_m\\rho+1}\\right)}{b_m(\\Delta_1-\\Delta_2)} $$\n",
+ "\n",
+ "with the constants given by:\n",
+ "\n",
+ "* vdW: $\\Delta_1=0$, $\\Delta_2=0$\n",
+ "* SRK: $\\Delta_1=1$, $\\Delta_2=0$\n",
+ "* PR: $\\Delta_1=1+\\sqrt{2}$, $\\Delta_2=1-\\sqrt{2}$\n",
+ "\n",
+ "The quantities $a_m$ and $b_m$ are described (with exact solutions for the numerical coefficients) for each of these EOS in https://pubs.acs.org/doi/abs/10.1021/acs.iecr.1c00847.\n",
+ "\n",
+ "The models in teqp are instantiated based on knowledge of the critical temperature, pressure, and acentric factor. Thereafter all quantities are obtained from derivatives of $\\alpha^r$."
+ ]
+ },
+ {
+ "cell_type": "raw",
+ "id": "a65e7316",
+ "metadata": {
+ "raw_mimetype": "text/restructuredtext"
+ },
+ "source": [
+ "The Python class is here: :py:class:`GeneralizedCubic <teqp.teqp.GeneralizedCubic>`"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9a884dab",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import teqp\n",
+ "\n",
+ "# Values taken from http://dx.doi.org/10.6028/jres.121.011\n",
+ "Tc_K = [ 190.564, 154.581, 150.687 ]\n",
+ "pc_Pa = [ 4599200, 5042800, 4863000 ]\n",
+ "acentric = [ 0.011, 0.022, -0.002 ]\n",
+ "\n",
+ "# Instantiate Peng-Robinson model\n",
+ "modelPR = teqp.canonical_PR(Tc_K, pc_Pa, acentric)\n",
+ "\n",
+ "# Instantiate Soave-Redlich-Kwong model\n",
+ "modelSRK = teqp.canonical_SRK(Tc_K, pc_Pa, acentric)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "615c7830",
+ "metadata": {},
+ "source": [
+ "## Adjusting k_ij\n",
+ "\n",
+ "Fine-tuned values of $k_{ij}$ can be provided when instantiating the model, for Peng-Robinson and SRK. A complete matrix of all the $k_{ij}$ values must be provided. This allows for asymmetric mixing models in which $k_{ij}\\neq k_{ji}$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a8c87890",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "k_12 = 0.01\n",
+ "kmat = [[0,k_12,0],[k_12,0,0],[0,0,0]]\n",
+ "teqp.canonical_PR(Tc_K, pc_Pa, acentric, kmat)\n",
+ "teqp.canonical_SRK(Tc_K, pc_Pa, acentric, kmat)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c3a62795",
+ "metadata": {},
+ "source": [
+ "## Superancillary\n",
+ "\n",
+ "The superancillary equation gives the co-existing liquid and vapor (orthobaric) densities as a function of temperature. The set of Chebyshev expansions was developed in https://pubs.acs.org/doi/abs/10.1021/acs.iecr.1c00847 . These superancillary equations are more accurate than iterative calculations in double precision arithmetic and also at least 10 times faster to calculate, and cannot fail in iterative routines, even extremely close to the critical point. \n",
+ "\n",
+ "The superancillary equation is only exposed for pure fluids to remove ambiguity when considering mixtures. The returned tuple is the liquid and vapor densities"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1b02ef72",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "teqp.canonical_PR([Tc_K[0]], [pc_Pa[0]], [acentric[0]]).superanc_rhoLV(100)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "da00ff20",
+ "metadata": {},
+ "source": [
+ "## a and b\n",
+ "\n",
+ "For the cubic EOS, it can be useful to obtain the a and b parameters directly. The b parameter is particularly useful because 1/b is the maximum allowed density in the EOS"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b5f6aaea",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "z = np.array([0.3, 0.4, 0.3])\n",
+ "modelPR.get_a(140, z), modelPR.get_b(140, z)"
+ ]
+ }
+ ],
+ "metadata": {
+ "celltoolbar": "Raw Cell Format",
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "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.10.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/doc/source/models/index.rst b/doc/source/models/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..305635865bbe29b0377261f7ec1ac6f97a564d80
--- /dev/null
+++ b/doc/source/models/index.rst
@@ -0,0 +1,11 @@
+Models
+======
+
+.. toctree::
+ :maxdepth: 2
+ :caption: Contents:
+
+ cubics
+ CPA
+ PCSAFT
+ multifluid
\ No newline at end of file
diff --git a/doc/source/models/multifluid.ipynb b/doc/source/models/multifluid.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..1aa94510216f3a22d57f02997b6e97de6a39b985
--- /dev/null
+++ b/doc/source/models/multifluid.ipynb
@@ -0,0 +1,347 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Multi-fluid EOS\n",
+ "\n",
+ "Peering into the innards of teqp"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2022-07-06T18:40:06.606089Z",
+ "iopub.status.busy": "2022-07-06T18:40:06.605465Z",
+ "iopub.status.idle": "2022-07-06T18:40:07.147629Z",
+ "shell.execute_reply": "2022-07-06T18:40:07.146669Z"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "'0.9.2'"
+ ]
+ },
+ "execution_count": 1,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "import timeit, json\n",
+ "import pandas\n",
+ "import numpy as np\n",
+ "import teqp\n",
+ "teqp.__version__"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Pure fluid loading"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2022-07-06T18:40:07.163785Z",
+ "iopub.status.busy": "2022-07-06T18:40:07.163785Z",
+ "iopub.status.idle": "2022-07-06T18:40:17.280030Z",
+ "shell.execute_reply": "2022-07-06T18:40:17.278333Z"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "12.9 ms ± 4.32 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
+ ]
+ }
+ ],
+ "source": [
+ "# By default teqp looks for fluids relative to the set of fluids in ROOT/dev/fluids\n",
+ "# The name (case-sensitive) should match the .json file, without the json extension.\n",
+ "%timeit model = teqp.build_multifluid_model([\"Methane\", \"Ethane\"], teqp.get_datapath())"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2022-07-06T18:40:17.286031Z",
+ "iopub.status.busy": "2022-07-06T18:40:17.285350Z",
+ "iopub.status.idle": "2022-07-06T18:40:27.300085Z",
+ "shell.execute_reply": "2022-07-06T18:40:27.299050Z"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "120 ms ± 19.9 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
+ ]
+ }
+ ],
+ "source": [
+ "# And if you provide valid aliases, alias lookup will be used to resolve the name\n",
+ "# But beware, this is rather a lot slower than the above because all fluid files need to be read\n",
+ "# in to build the alias map\n",
+ "%timeit model = teqp.build_multifluid_model([\"n-C1H4\", \"n-C3H8\"], teqp.get_datapath())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "So, how to make it faster? Only do it once and cache"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2022-07-06T18:40:27.307115Z",
+ "iopub.status.busy": "2022-07-06T18:40:27.307115Z",
+ "iopub.status.idle": "2022-07-06T18:40:34.876155Z",
+ "shell.execute_reply": "2022-07-06T18:40:34.875146Z"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "85.4 ms ± 12.3 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Here is the set of possible aliases to absolute paths of files\n",
+ "# Building this map takes a little while (somewhat faster in C++) due to all the file reads\n",
+ "# If you know your files will not change, good idea to build this alias map yourself.\n",
+ "%timeit aliasmap = teqp.build_alias_map(teqp.get_datapath())\n",
+ "aliasmap = teqp.build_alias_map(teqp.get_datapath())"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2022-07-06T18:40:34.879157Z",
+ "iopub.status.busy": "2022-07-06T18:40:34.879157Z",
+ "iopub.status.idle": "2022-07-06T18:40:42.930387Z",
+ "shell.execute_reply": "2022-07-06T18:40:42.930079Z"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "9.9 ms ± 164 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Then load the absolute paths from the alias map, \n",
+ "# which will guarantee that you hit exactly what you were looking for,\n",
+ "# resolving aliases as needed\n",
+ "identifiers = [aliasmap[n] for n in [\"Neon\", \"Hydrogen\"]]\n",
+ "%timeit model = teqp.build_multifluid_model(identifiers, teqp.get_datapath())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "At some point soon teqp will support in-memory loading of JSON data for the pure components, without requiring reads from the operating system"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2022-07-06T18:40:42.932720Z",
+ "iopub.status.busy": "2022-07-06T18:40:42.932720Z",
+ "iopub.status.idle": "2022-07-06T18:40:42.946939Z",
+ "shell.execute_reply": "2022-07-06T18:40:42.945927Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "# And you can also load the JSON that teqp is loading for the pure fluids\n",
+ "pureJSON = teqp.collect_component_json(['Neon','Hydrogen'], teqp.get_datapath())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Mixture model loading"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2022-07-06T18:40:42.949975Z",
+ "iopub.status.busy": "2022-07-06T18:40:42.948972Z",
+ "iopub.status.idle": "2022-07-06T18:40:42.962308Z",
+ "shell.execute_reply": "2022-07-06T18:40:42.961547Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "# Load the default JSON for the binary interaction parameters\n",
+ "BIP = json.load(open(teqp.get_datapath()+'/dev/mixtures/mixture_binary_pairs.json'))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2022-07-06T18:40:42.965595Z",
+ "iopub.status.busy": "2022-07-06T18:40:42.964944Z",
+ "iopub.status.idle": "2022-07-06T18:40:42.978409Z",
+ "shell.execute_reply": "2022-07-06T18:40:42.977232Z"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "{'BibTeX': 'Kunz-JCED-2012',\n",
+ " 'CAS1': '74-82-8',\n",
+ " 'CAS2': '74-84-0',\n",
+ " 'F': 1.0,\n",
+ " 'Name1': 'Methane',\n",
+ " 'Name2': 'Ethane',\n",
+ " 'betaT': 0.996336508,\n",
+ " 'betaV': 0.997547866,\n",
+ " 'function': 'Methane-Ethane',\n",
+ " 'gammaT': 1.049707697,\n",
+ " 'gammaV': 1.006617867}"
+ ]
+ },
+ "execution_count": 8,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# You can obtain interaction parameters either by pairs of names, where name is the name that teqp uses, the [\"INFO\"][\"NAME\"] field\n",
+ "params, swap_needed = teqp.get_BIPdep(BIP, ['Methane','Ethane'])\n",
+ "params"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2022-07-06T18:40:42.981288Z",
+ "iopub.status.busy": "2022-07-06T18:40:42.981288Z",
+ "iopub.status.idle": "2022-07-06T18:40:42.993879Z",
+ "shell.execute_reply": "2022-07-06T18:40:42.993281Z"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "{'BibTeX': 'Kunz-JCED-2012',\n",
+ " 'CAS1': '74-82-8',\n",
+ " 'CAS2': '74-84-0',\n",
+ " 'F': 1.0,\n",
+ " 'Name1': 'Methane',\n",
+ " 'Name2': 'Ethane',\n",
+ " 'betaT': 0.996336508,\n",
+ " 'betaV': 0.997547866,\n",
+ " 'function': 'Methane-Ethane',\n",
+ " 'gammaT': 1.049707697,\n",
+ " 'gammaV': 1.006617867}"
+ ]
+ },
+ "execution_count": 9,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# Or also by CAS#\n",
+ "params, swap_needed = teqp.get_BIPdep(BIP, ['74-82-8','74-84-0'])\n",
+ "params"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2022-07-06T18:40:42.996338Z",
+ "iopub.status.busy": "2022-07-06T18:40:42.996338Z",
+ "iopub.status.idle": "2022-07-06T18:40:43.423368Z",
+ "shell.execute_reply": "2022-07-06T18:40:43.422895Z"
+ },
+ "tags": [
+ "raises-exception"
+ ]
+ },
+ "outputs": [
+ {
+ "ename": "ValueError",
+ "evalue": "Can't match the binary pair for: 74-82-8/Ethane",
+ "output_type": "error",
+ "traceback": [
+ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
+ "\u001b[1;31mValueError\u001b[0m Traceback (most recent call last)",
+ "Input \u001b[1;32mIn [10]\u001b[0m, in \u001b[0;36m<cell line: 2>\u001b[1;34m()\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[38;5;66;03m# But mixing is not allowed\u001b[39;00m\n\u001b[1;32m----> 2\u001b[0m params, swap_needed \u001b[38;5;241m=\u001b[39m \u001b[43mteqp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget_BIPdep\u001b[49m\u001b[43m(\u001b[49m\u001b[43mBIP\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43m74-82-8\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mEthane\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 3\u001b[0m params\n",
+ "\u001b[1;31mValueError\u001b[0m: Can't match the binary pair for: 74-82-8/Ethane"
+ ]
+ }
+ ],
+ "source": [
+ "# But mixing is not allowed\n",
+ "params, swap_needed = teqp.get_BIPdep(BIP, ['74-82-8','Ethane'])\n",
+ "params"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "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.10.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}