{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Implementing an equation of state in python\n",
"\n",
"> In `FeOs`, you can implement your equation of state in python, register it to the Rust backend, and compute properties and phase equilbria as if you implemented it in Rust.\n",
"> In this tutorial, we will implement the Peng-Robinson equation of state."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Table of contents \n",
"\n",
"- [Implementation](#Implementation)\n",
"- [Computing properties](#Computing-properties)\n",
"- [Critical point](#Critical-point)\n",
"- [Phase equilibria and phase diagrams](#Phase-equilibria-and-phase-diagrams)\n",
"- [Mixtures](#Mixtures)\n",
"- [Comparison to Rust implementation](#Comparison-to-Rust-implementation)"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from feos.eos import *\n",
"from feos.si import *\n",
"import numpy as np\n",
"\n",
"optional = True\n",
"\n",
"import warnings\n",
"warnings.filterwarnings('ignore')\n",
"\n",
"if optional:\n",
" import matplotlib.pyplot as plt\n",
" import pandas as pd\n",
" import seaborn as sns\n",
"\n",
" sns.set_style(\"ticks\")\n",
" sns.set_palette(\"Dark2\")\n",
" sns.set_context(\"talk\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Implementation \n",
"[↑ Back to top](#toc)\n",
"\n",
"To implement an equation of state in python, we have to define a `class` which has to have the following methods:\n",
"\n",
"```python\n",
"class MyEquationOfState:\n",
" def helmholtz_energy(self, state: StateHD) -> D\n",
" \n",
" def components(self) -> int\n",
" \n",
" def subset(self, indices: List[int]) -> Self\n",
" \n",
" def molar_weight(self) -> SIArray1\n",
" \n",
" def max_density(self, moles: SIArray1) -> f64\n",
"``` \n",
"\n",
"- `components(self) -> int`: Returns the number of components (usually inferred from the shape of the input parameters).\n",
"- `molar_weight(self) -> SIArray1`: Returns an `SIArray1` with size equal to the number of components containing the molar mass of each component.\n",
"- `max_density(self, moles: np.ndarray[float]) -> float`: Returns the maximum allowed number density in units of `molecules/Angstrom³`.\n",
"- `subset(self, indices: List[int]) -> self`: Returns a new equation of state with parameters defined in `indices`.\n",
"- `helmholtz_energy(self, state: StateHD) -> Dual`: Returns the helmholtz energy as (hyper)-dual number given a `StateHD`."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"SQRT2 = np.sqrt(2)\n",
"\n",
"class PyPengRobinson: \n",
" def __init__(\n",
" self, critical_temperature, critical_pressure, \n",
" acentric_factor, molar_weight, delta_ij=None\n",
" ):\n",
" \"\"\"Peng-Robinson Equation of State\n",
" \n",
" Parameters\n",
" ----------\n",
" critical_temperature : SIArray1\n",
" critical temperature of each component.\n",
" critical_pressure : SIArray1\n",
" critical pressure of each component.\n",
" acentric_factor : np.array[float] \n",
" acentric factor of each component (dimensionless).\n",
" molar_weight: SIArray1\n",
" molar weight of each component.\n",
" delta_ij : np.array[[float]], optional\n",
" binary parameters. Shape=[n, n], n = number of components.\n",
" defaults to zero for all binary interactions.\n",
" \n",
" Raises\n",
" ------\n",
" ValueError: if the input values have incompatible sizes.\n",
" \"\"\"\n",
" self.n = len(critical_temperature)\n",
" if len(set((\n",
" len(critical_temperature), \n",
" len(critical_pressure), \n",
" len(acentric_factor)\n",
" ))) != 1:\n",
" raise ValueError(\"Input parameters must all have the same lenght.\")\n",
" \n",
" self.tc = critical_temperature / KELVIN\n",
" self.pc = critical_pressure / PASCAL\n",
" self.omega = acentric_factor\n",
" self.mw = molar_weight / GRAM * MOL\n",
"\n",
" self.a_r = (0.45724 * critical_temperature**2 * RGAS \n",
" / critical_pressure / ANGSTROM**3 / NAV / KELVIN)\n",
" self.b = (0.07780 * critical_temperature * RGAS \n",
" / critical_pressure / ANGSTROM**3 / NAV)\n",
" self.kappa = (0.37464 \n",
" + (1.54226 - 0.26992 * acentric_factor) * acentric_factor)\n",
" self.delta_ij = (np.zeros((self.n, self.n)) \n",
" if delta_ij is None else delta_ij)\n",
" \n",
" def helmholtz_energy(self, state):\n",
" \"\"\"Return helmholtz energy.\n",
" \n",
" Parameters\n",
" ----------\n",
" state : StateHD\n",
" The thermodynamic state.\n",
" \n",
" Returns\n",
" -------\n",
" helmholtz_energy: float | any dual number\n",
" The return type depends on the input types.\n",
" \"\"\" \n",
" n = np.sum(state.moles)\n",
" x = state.molefracs\n",
" tr = 1.0 / self.tc * state.temperature\n",
" ak = ((1.0 - np.sqrt(tr)) * self.kappa + 1.0)**2 * self.a_r\n",
" ak_mix = 0.0\n",
" if self.n > 1:\n",
" for i in range(self.n):\n",
" for j in range(self.n):\n",
" ak_mix += (np.sqrt(ak[i] * ak[j]) \n",
" * (x[i] * x[j] * (1.0 - self.delta_ij[i, j])))\n",
" else:\n",
" ak_mix = ak[0]\n",
" b = np.sum(x * self.b)\n",
" v = state.volume\n",
" rho = np.sum(state.partial_density)\n",
" a = n * (-np.log(1.0 - b * rho) \n",
" - ak_mix / (b * SQRT2 * 2.0 * state.temperature) \n",
" * np.log((1.0 + (1.0 + SQRT2) * b * rho) / (1.0 + (1.0 - SQRT2) * b * rho)))\n",
" return a\n",
" \n",
" def components(self) -> int: \n",
" \"\"\"Number of components.\"\"\"\n",
" return self.n\n",
" \n",
" def subset(self, i: [int]):\n",
" \"\"\"Return new equation of state containing a subset of all components.\"\"\"\n",
" if self.n > 1:\n",
" tc = self.tc[i] \n",
" pc = self.pc[i]\n",
" mw = self.mw[i]\n",
" omega = self.omega[i]\n",
" return PyPengRobinson(\n",
" SIArray1(tc * KELVIN), \n",
" SIArray1(pc * PASCAL), \n",
" omega, \n",
" SIArray1(mw * GRAM / MOL)\n",
" )\n",
" else:\n",
" return self\n",
" \n",
" def molar_weight(self) -> SIArray1:\n",
" return SIArray1(self.mw * GRAM / MOL)\n",
" \n",
" def max_density(self, moles:[float]) -> float:\n",
" b = np.sum(moles * self.b) / np.sum(moles);\n",
" return 0.9 / b "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Computing properties \n",
"[↑ Back to top](#toc)\n",
"\n",
"Let's compute some properties. First, we have to instanciate the class and register it to Rust.\n",
"This is done using the `EquationOfState.python` method."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"tc = SIArray1(369.96 * KELVIN)\n",
"pc = SIArray1(4250000.0 * PASCAL)\n",
"omega = np.array([0.153])\n",
"molar_weight = SIArray1(44.0962 * GRAM / MOL)\n",
"\n",
"# create an instance of our python class and hand it over to rust\n",
"pr = PyPengRobinson(tc, pc, omega, molar_weight)\n",
"eos = EquationOfState.python_residual(pr)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Thermodynamic state: the `State` object\n",
"\n",
"Before we can compute a property, we create a `State` object. This can be done in several ways depending on what control variables we need.\n",
"If no total amount of substance is defined, it is set to $n = \\frac{1}{N_{AV}}$.\n",
"For possible input combinations, you can inspect the signature of the constructor using `State?`."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"\u001b[0;31mInit signature:\u001b[0m \u001b[0mState\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m/\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mDocstring:\u001b[0m \n",
"A thermodynamic state at given conditions.\n",
"\n",
"Parameters\n",
"----------\n",
"eos : Eos\n",
" The equation of state to use.\n",
"temperature : SINumber, optional\n",
" Temperature.\n",
"volume : SINumber, optional\n",
" Volume.\n",
"density : SINumber, optional\n",
" Molar density.\n",
"partial_density : SIArray1, optional\n",
" Partial molar densities.\n",
"total_moles : SINumber, optional\n",
" Total amount of substance (of a mixture).\n",
"moles : SIArray1, optional\n",
" Amount of substance for each component.\n",
"molefracs : numpy.ndarray[float]\n",
" Molar fraction of each component.\n",
"pressure : SINumber, optional\n",
" Pressure.\n",
"molar_enthalpy : SINumber, optional\n",
" Molar enthalpy.\n",
"molar_entropy : SINumber, optional\n",
" Molar entropy.\n",
"molar_internal_energy: SINumber, optional\n",
" Molar internal energy\n",
"density_initialization : {'vapor', 'liquid', SINumber, None}, optional\n",
" Method used to initialize density for density iteration.\n",
" 'vapor' and 'liquid' are inferred from the maximum density of the equation of state.\n",
" If no density or keyword is provided, the vapor and liquid phase is tested and, if\n",
" different, the result with the lower free energy is returned.\n",
"initial_temperature : SINumber, optional\n",
" Initial temperature for temperature iteration. Can improve convergence\n",
" when the state is specified with pressure and molar entropy or enthalpy.\n",
"\n",
"Returns\n",
"-------\n",
"State : state at given conditions\n",
"\n",
"Raises\n",
"------\n",
"Error\n",
" When the state cannot be created using the combination of input.\n",
"\u001b[0;31mType:\u001b[0m type\n",
"\u001b[0;31mSubclasses:\u001b[0m \n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"State?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we use input variables other than $\\mathbf{N}, V, T$ (the natural variables of the Helmholtz energy), creating a state is an iterative procedure.\n",
"For example, we can create a state for a give $T, p$, which will result in a iteration of the volume (density)."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$1.6605\\times10^{-24}\\,\\mathrm{ mol}$"
],
"text/plain": [
"1.6605390671738466e-24 mol"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# If no amount of substance is given, it is set to 1/NAV.\n",
"s = State(eos, temperature=300*KELVIN, pressure=1*BAR)\n",
"s.total_moles"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$1\\,\\mathrm{ mol}$"
],
"text/plain": [
"1 mol"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s_pt = State(\n",
" eos, \n",
" temperature=300*KELVIN, \n",
" pressure=1*BAR, \n",
" total_moles=1*MOL\n",
")\n",
"s_pt.total_moles"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can use other variables as well. For example, we can create a state at given $h, p$ (using the enthalpy from the prior computation as input):"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"thread '' panicked at 'No ideal gas model initialized!', src/eos.rs:"
]
},
{
"ename": "PanicException",
"evalue": "No ideal gas model initialized!",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mPanicException\u001b[0m Traceback (most recent call last)",
"Input \u001b[0;32mIn [9]\u001b[0m, in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0m s_ph \u001b[38;5;241m=\u001b[39m \u001b[43mState\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m 2\u001b[0m \u001b[43m \u001b[49m\u001b[43meos\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\n\u001b[1;32m 3\u001b[0m \u001b[43m \u001b[49m\u001b[43mpressure\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mBAR\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\n\u001b[1;32m 4\u001b[0m \u001b[43m \u001b[49m\u001b[43mmolar_enthalpy\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43ms_pt\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmolar_enthalpy\u001b[49m\u001b[43m(\u001b[49m\u001b[43mContributions\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mResidual\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 5\u001b[0m \u001b[43m)\u001b[49m\n",
"\u001b[0;31mPanicException\u001b[0m: No ideal gas model initialized!"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"49:10\n"
]
}
],
"source": [
"s_ph = State(\n",
" eos, \n",
" pressure=1*BAR, \n",
" molar_enthalpy=s_pt.molar_enthalpy(Contributions.Residual)\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"rel. dev.\n",
"entropy : 1.2028878239185388e-16\n",
"density : -1.8532974044002563e-15\n",
"temperature: 1.5158245029548805e-15\n"
]
}
],
"source": [
"# check if states are equal\n",
"print(\"rel. dev.\")\n",
"print(\n",
" \"entropy : \", \n",
" (s_ph.molar_entropy() - s_pt.molar_entropy()) / s_pt.molar_entropy()\n",
")\n",
"print(\n",
" \"density : \", \n",
" (s_ph.mass_density() - s_pt.mass_density()) / s_pt.mass_density()\n",
")\n",
"print(\n",
" \"temperature: \",\n",
" (s_ph.temperature - s_pt.temperature) / s_pt.temperature\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Critical point \n",
"[↑ Back to top](#toc)\n",
"\n",
"To generate a state at critical conditions, we can use the `critical_point` constructor."
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Critical point\n",
"temperature: 369.9506174234607 K\n",
"density : 198.1862458057177 kg/m³\n",
"pressure : 4.249677749116942 MPa\n"
]
}
],
"source": [
"s_cp = State.critical_point(eos)\n",
"print(\"Critical point\")\n",
"print(\"temperature: \", s_cp.temperature)\n",
"print(\"density : \", s_cp.mass_density())\n",
"print(\"pressure : \", s_cp.pressure())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Phase equilibria and phase diagrams\n",
"[↑ Back to top](#toc)\n",
"\n",
"We can also create an object, `PhaseEquilibrium`, that contains states that are in equilibrium."
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/markdown": [
"||temperature|density|\n",
"|-|-|-|\n",
"|phase 1|300.00000 K|488.99014 mol/m³|\n",
"|phase 2|300.00000 K|11.53399 kmol/m³|\n"
],
"text/plain": [
"phase 0: T = 300.00000 K, ρ = 488.99014 mol/m³\n",
"phase 1: T = 300.00000 K, ρ = 11.53399 kmol/m³"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"vle = PhaseEquilibrium.pure(eos, 300.0*KELVIN)\n",
"vle"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Each phase is a `State` object. We can simply access these states and compute properties, just like before."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/markdown": [
"|temperature|density|\n",
"|-|-|\n",
"|300.00000 K|11.53399 kmol/m³|"
],
"text/plain": [
"T = 300.00000 K, ρ = 11.53399 kmol/m³"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"vle.liquid # the high density phase `State`"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/markdown": [
"|temperature|density|\n",
"|-|-|\n",
"|300.00000 K|488.99014 mol/m³|"
],
"text/plain": [
"T = 300.00000 K, ρ = 488.99014 mol/m³"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"vle.vapor # the low density phase `State`"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Heat of vaporization: 14.782343503305114 kJ/mol\n",
"for T = 300 K\n",
"and p = 9.95 bar\n"
]
}
],
"source": [
"# we can now easily compute any property:\n",
"print(\"Heat of vaporization: \", vle.vapor.molar_enthalpy() - vle.liquid.molar_enthalpy())\n",
"print(\"for T = {}\".format(vle.liquid.temperature))\n",
"print(\"and p = {:.2f} bar\".format(vle.liquid.pressure() / BAR))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can also easily compute **vapor pressures** and **boiling temperatures**:"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"vapor pressure (T = 300 K): 994.7761635610093 kPa\n",
"boiling temperature (p = 3 bar): 247.84035574956746 K\n"
]
}
],
"source": [
"# This also works for mixtures, in which case the pure component properties are computed.\n",
"# Hence, the result is a list - that is why we use an index [0] here.\n",
"print(\"vapor pressure (T = 300 K):\", PhaseEquilibrium.vapor_pressure(eos, 300*KELVIN)[0])\n",
"print(\"boiling temperature (p = 3 bar):\", PhaseEquilibrium.boiling_temperature(eos, 2*BAR)[0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Phase Diagram\n",
"\n",
"We could repeatedly compute `PhaseEquilibrium` states for different temperatures / pressures to generate a phase diagram.\n",
"Because this a common task, there is a object for that as well.\n",
"\n",
"The `PhaseDiagram` object creates multiple `PhaseEquilibrium` objects (`npoints` of those) between a given lower temperature and the critical point."
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [],
"source": [
"dia = PhaseDiagram.pure(eos, 230.0 * KELVIN, 500)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can access each `PhaseEquilbrium` and then conveniently comput any property we like:"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [],
"source": [
"enthalpy_of_vaporization = [\n",
" (vle.vapor.molar_enthalpy() - vle.liquid.molar_enthalpy()) / (KILO * JOULE) * MOL \n",
" for vle in dia.states\n",
"]"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAd0AAAEZCAYAAADSRP/MAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6e0lEQVR4nO3deVxVdf4/8Ne97PsuyCIiCAgquICKBini7rjk16xJrdScJvu5pNVMy8yUbTOlpblmWmaZpbkwlgriLqK4IyCrbAqCsu9wz+8P5I4I4mW599x7eT0fDx/pueec+/50PLz8nOXzkQiCIICIiIiUTip2AURERF0FQ5eIiEhFGLpEREQqwtAlIiJSEYYuERGRiuiKXYC68fHxgUwmg6mpqdilEBGRGigrK4NUKkV8fHyH98We7iNkMhn4FhURETUSBAEymaxT9sWe7iMae7ixsbEiV0JEROpg8ODBnbYv9nSJiIhUhKFLRESkIgxdIiIiFWHoEhERqQhDl4iISEUYukRERCrCV4aUoKquFqfvpMBc3whelvawMDASuyQiIlIDDF0l2HEzBv88/1/5n7sbW8DLyh7eVg7wfvBfD4tuMNTVE7FKIiJSNYauEgzu5ope5rZIKykAANypKMadimIcz0mSryOVSOBmbgsvy4YQ9rKyRx8rB7ia2UBHyqv+RETaSCJwzMMmGkce6YwRqUprqpBUlIeEwlzcLMxDYmEuEgtzUVhd8dhtDHR04Wlp3ySMva0c4GBsDolE0uGaiIiobTozF9jTVSIzfUMM6uaKQd1c5csEQUB+ZRluFuU+COGGME4qykNlXS2q6+tw/V4Ort/LabIvC30j+aVpH2tH+Fh3h7eVPYx09VXdLCIiaieGropJJBJ0MzZDN2MzPOXYW75cJsiQWVqImw96wzeLGsI4tbgA9YIMxTWViMm7hZi8W/JtpBIJepnbwse6uzyIfa27o5uRGXvFRERqiKGrJqQSKXqa26CnuQ3GuvrKl1fX1yG1OF9+eTr+/h3EF95BXkUJZIKAlOJ8pBTn40D6Nfk2NoYm8LV2RB/r7vIgdrewg55UR4ymERHRAwxdNWego/ugJ9u9yfJ7VWUNAXz/Dm48+G9K0V3UCTLcqyrHydvJOHk7Wb6+vlQHXlYO8LH+3+VpH6vufJ2JiEiFGLoaysbQFE859m5yibq6vg7JRXnyEL5x/zYS7t9BcU0VamT1D90rvijfxtXMBv1sHNHPxunBL0dYGZqI0CIiIu3H0NUiBjq66GvjhL42TvJlgiDgdnkxbty/LQ/i+Pu5yCi9BwDIKL2HjNJ7+O+t6/JtnE0tHwphJ/S3dYKNoanK20NEpG0YulpOIpHAydQSTqaWGNPDR768tKYK8ffvyHu/1+/lIKU4HzJBQHZZEbLLivBHxg35+t2NLdDf1gl9bRzR38YZfW0cYW9sLkaTiIg0FkO3izLTN8QQBzcMcXCTL6uorWkSxNfu5SC56C7qBVnDAB+ZxTicGS9f397IDH1tnNDP1gn+ts7ws3WGnZGZGM0hItIIDF2SM9bTx2B7Vwy2/997xZV1tUgovIO4e7dxrSAbcfduI7EwF3WCDHmVpcjLTsTR7ET5+s6mlvC3dYGfrTP87VzQ38YJJnoGYjSHiEjtMHSpVUa6ehho1wMD7XrIl1XX1+FmYS6u3cvB9YKGXnFCYS5qZfXyS9ON94ilEgl6W3SDv50LBti6wN/OGV5WDnx9iYi6JIYutZmBji762zqjv60z4NWwrKquFvGFd3AlPwtXCrJxJT8LaSUFkAkCbhbl4WZRHnYlx8q372fTeEm6IYh7mtlwQA8i0noMXeoUhi30iIurK3HtXjau5GfjSkEWLudn4W5lKarr6xB7NwOxdzPk61oaGGOArQsGdeuBwd1c4W/nAlNeliYiLcPQJaWxMDBq8i6xIAi4U1HyoDechasF2bhakI2y2moUVVfgWM5NHMu5CaDhsrSXpT0Gd3PFoG49MNDOFW7m7A0TkWZj6JLKSCQSOJpYwNHEAhN69gXQMOZ0anEBLudn4uLdTFzMz8TNwjzIBAEJhblIKMzFDzdjAADWBiYY1K3Hg1+u8Ld15oQPRKRROLXfIzpzCidqn5KaKlzJz0Ls3QxczM/E5fxMlNRUNVtPRyKFr3V3DOzWA4PsXBFg7wpnUysRKiYibdaZucDQfQRDV/3IBBmSi/JxMT8DF+9m4OLdTKQU57e4rqOJBQLt3TDEvicC7d3Q29IOUolUxRUTkTZh6CoRQ1czFFZX4NLdTFzKz0Ts3Qxcyc9CeV1Ns/WsDIwRaN/zwS839LVx5OtKRNQmnMSeujwrA2OEungj1MUbAFAnq0fC/VzE5KUjJi8d5/Nu4V5VOQqrK3A4M14+kpaxrj4G2vXAEIeGEB5o58L7wkSkMuzpPoI9Xe0gCALSSgoaAjj3Fs7n3UJm2f1m6+lJddDPxgmB9j0xzKEXAu17wkzfUISKiUhd8fKyEjF0tdft8mJcyLvV0BvOTcfNorxm6+hIpOhv64QgB3cEde+FgG49YazHnjBRV8bQVSKGbtdRWF2B2LxbiHkQxNcKclAvyJqsoyfVgb+tM4K6uyPIoRcGdnOFka6eSBUTkRgYukrE0O26ymqrEZObjrO5aYi+k4q4+7che+T0MNDRxQA7Fwzv7o4gB3cMsHOBvg4fjSDSZgxdJWLoUqPi6krE5KXj7J1UnM1NQ/z9O83WMdTRQ4C9K4Z390CIY2/42nTnK0pEWoahq0QMXXqc+1XliM5Nw9k7aYjOTUVS0d1m61gbmGCEozuecuyNYMfecDK1VH2hRNSpGLpKxNAlReVXliL6ThrO3EnFqdspLT4d7W5hh6ccG3rBw7q7cxIHIg3E0FUihi61162Sezh5Oxknc5JxNje12dCVuhIpBnbrgeAHk0D42TpBlwN1EKk9hq4SMXSpM9TJ6nG1IAcnbyfh1O0UXLqbibpHnow21zfEiO4eGOXshaedveBgbC5StUTUGoauEjF0SRlKa6oQnZuGk7eTcep2ClJbGDva17o7Rjl7Y5SzFwbYubAXTKQmGLpKxNAlVcgpK8KJ20k4np2Ek7eTUVZb3eRzC30jhDj1bugFO3nB1shUpEqJiKGrRAxdUrVaWT1i72YgKusmorITm42UJYEE/W2dMMrZCyOdveBn4wwdKV9LIlIVhq4SMXRJbDllRTiWcxPHsm/i1O0UVDwye5K1gQlCXbwwpocPgh17w4RPRBMpFUNXiRi6pE6q6+twPi8dx7JvIir7ZrN5hA10dDG8uzvG9PDBaJc+fBiLSAm0PnRzc3OxZcsW3LhxA4mJiaioqMD27dsxZMiQJuvNnj0b58+fb7b9hAkTsHr16nZ9N0OX1FlG6T0czUpERFYCou+kNXsi2s/WGWNc+iCshw/6WDlAIpGIVCmR9tD6+XQzMjJw8OBB+Pj4YOjQoYiKinrsuj179sRnn33WZJmVlZWySyQShauZDV72GY6XfYajpKYKx7Jv4khWPI5l30RJTRWuFmTjakE2/nM5Ai6mVghz6YMxPXwwxMENenwamkh0ahm6AQEBiI6OBgBERka2GrqGhobw9/dXUWVE6sNc3xBTevlhSi8/1MrqEZObjiOZ8TiSFY/ssiJklRVia8JZbE04Cwt9Q4x26YMJrn0R7OTJmZKIRKKWoSvlk5lEbaIn1cEIRw+McPTAv4ZMRkJhLiIy43EkKwFXC7JRXFOFPamXsSf1Mkx09RHq4o0Jrn0x0tmLD2IRqZBCobtv37527Xzq1Knt2q4t0tPTERAQgPLycjg7O2Pq1KlYsGAB9PRa/pd847X5xyktLYWZmZkySiVSCYlEAh/r7vCx7o7F/qG4U16MI5nx+D0jDtG5aSivq8GB9Gs4kH4NBjq6eNrJExN69kOYSx+Y6xuKXT6RVlModN9++21IJBK05ZkriUSi9NAdNGgQJkyYgF69eqGiogKRkZFYs2YNbty4gXXr1in1u4k0RXcTC8ztMwxz+wzD/apyHM6Mx++34nD6Tgqq6+twODMehzPjoSfVwVOOHpjg2hdje/jAytBE7NKJtI5Cobt9+3Zl19EuS5YsafLnkSNHwtbWFhs3bkRsbGyLvdonPX32pJ4wkSazNjTBc54BeM4zAMXVlYjISsAfGXE4npOE6vo6RD14Nemts3vxlKMHprj5YayrL3vARJ1EodANDAxUdh2dZurUqdi4cSOuXLnCACVqhYWBEWZ4DMQMj4Eor61GVPZN/H4rDkezE1FRV4PjOUk4npME/bO/YaSzF6a4+WG0Sx8Y6+mLXTqRxlLLB6k6QiZreG+RD2MRKc5EzwCT3fpjslt/VNbV4lj2TRxIv4qIrIQml6CNdfUR1qMPprj5IcTJEwY6WvcjhEip2n3GVFRUYMuWLYiIiEB2djYAwNnZGWPGjMG8efNgbGzcaUW2xf79+wEAfn5+onw/kaYz0tXDhJ59MaFnX5TVViMiMwEH0q/ieE4SKupqsD/tKvanXYW5viHGu/riT25+GN7dnbMiESmgXaFbVFSEP//5z0hNTYWVlRX69OkDALh16xbWrVuHQ4cO4ccff4SlpWW7Czt06BAA4Pr16wCACxcuoLCwEEZGRggJCUFsbCw2b96MMWPGwMnJCRUVFTh69Ch+++03jBs3DoMGDWr3dxNRA1M9A0xz98c0d38UVlfgUMYNHEi/ijN3UlFSU4VdyRexK/ki7IxMMbWXP2a4D4SPdXeOhEX0GO0aBvKDDz7Azp078e6772LWrFnQ0Wn4F259fT127dqFlStX4vnnn8e7777b7sK8vLxaXO7k5ISoqChkZGTgo48+QmJiIgoLCyGVSuHm5oapU6di9uzZ8praisNAEj1ZfmUpDt6Kw4H0qzifd6vJZ16W9pjhMRBTe/mju4mFOAUSdSLRx15++umnERwcjA8++KDFz9977z2cOnUKx48f72h9KsfQJWqb7LJC7E29gt2pl5D60IQMEkjwlKMHprsPwHhXXw7CQRpL9LGXCwoK5JeUW+Lj44O9e/e2uygi0hzOplZ43W8kFvV/GlcLsrE79RL2p11FYXUFTt5Oxsnbyfh7tD7Gu/pihsdABDm4cz5g6rLaFbq2trZISEh47OcJCQmwtbVtd1FEpHkkEgn87Vzgb+eC9wMm4nhOEvakXkZEZjwq6mrkw1A6GJtjpscgPOs5GK5mNmKXTaRS7QrdkSNHYteuXfDx8cHMmTPlr+fIZDL8+uuv2LNnD5599tlOLZSINIe+ji7G9PDBmB4+KKquwH9vXceelEu4cDcDuRUlWHPtGNZcO4bh3d0xq3cAxrv6wpCTMFAX0K57uoWFhZg1axYyMzNhbW0NNzc3AA3jIN+/fx89evTAzz//rJFT7PGeLpHy3Cq5h1+SY/FLykXkVpTIl1voG2Gauz+e6x0AXxtHESskak70B6kAoKysDN988w0iIyPl7+m6uLggNDQUCxYsgKmpaYeLEwNDl0j56mUyHM9Jws6kC4jMSkCdIJN/1s/GCbN6D8bUXv6wMDASsUqiBmoRutqKoUukWvmVpdiTchk/J19AykNPPxvo6GJSz36Y4z0MA+1c+O4viYahq0QMXSJxCIKAi3czsTP5PA6kX0NlXa38M1/r7pjtPRTTevnz1SNSOYauEjF0icRXVluNfWlXsCMxBnH3b8uXm+oZ4Bn3gZjjPRReVvYiVkhdiVqEbnh4OH788UdkZGSgqKio+Y4lEsTHx3e0PpVj6BKpD0EQcLkgC9sTziH81jVU19fJPxti3xOzvYdivGtfTrxASiX64Bjr16/H2rVrYWNjgwEDBsDCgkO9EVHnk0gkGGjXAwPteuAfgRPxS8pF/JAYg1ul9xCTdwsxebdgY2iC5z0DMcd7KIedJLXXrp7uiBEj4O7uji1btkBPT7verWNPl0i9yQQZTt9OxfbEc4jISkD9gyefdSVSTHTrh/k+IzDAzkXkKkmbiN7TLS8vx/jx47UucIlI/UklUgQ79UawU2/cLi/GjzdjsONmDO5VlcunHRxo1wPzfYZjfM++0OOUg6RG2hW6ffr0wZ07dzq7FiKiNnE0scCKgWPwev+ROJB+FVvizyD+/h1cys/EX09kovsFC8ztMwx/9gyAlaGJ2OUSte/y8vnz5/H6669j69at8PX1VUZdouHlZSLNJQgCzuWlY8uN0ziSmQABDT/eDHX08Iz7AMz3HYHelt1ErpI0jVo8vRwZGYklS5bA398fTk5O8vGX5TuWSPDxxx93uEBVY+gSaYeM0nv4LiEaPyddQGlttXz5GJc+eLVfCALse4pXHGkU0UP36tWrmD9/PkpLSx+/Y4mk1ZmI1BVDl0i7lNVW45fkWHwbfxYZpffkywd3c8Vf+4VgtIs3pBJONUiPJ3rozpw5E1lZWfjoo48wePBgmJubd7gQdcHQJdJO9TIZ/si8gQ3XT+BqQbZ8uYeFHf7SNxjT3AfwfV9qUWfmQrv+eXfz5k28/PLLGDVqlFYFLhFpLx2pFJN69sN/J72GX8YtwNNOngCAlOJ8LD+zB0G/fob110+gpKZK5EpJm7UrdG1sbPi6EBFpJIlEgqDu7tgx5mUcmbIY090HQEciRV5lKT6O/QNDf/0U/7l0BIXVFWKXSlqoXaE7ffp0HDhwAHV1dU9emYhITflYd8ea4GdxZsYKzPcZDmNdfZTUVOGrq1EY+sun+Dj2DxRUloldJmmRdt3TjY6OxhdffAGZTIbnn38ezs7O0NFp/gJ6QEBApxSpSrynS9R1FVaVY0v8GWyNPyN/4tlQRw9zvIdgYd9g2BvzdlpXJPqDVN7e3k138sg8l4Ig8OllItJYxdWV2JpwBltunEFxTSWAhvl9n/MMwF/7hsDR1FLcAkmlRA/dvXv3KrTetGnT2lyQ2Bi6RNSotKYK3yeew+a4U7hfXQ4A0JPq4HnPALzuNwoO7Pl2CaKHrjZj6BLRoypqa7DjZgw2xp3E3cqG8QkMdHTxYp8gvNYvBNYcYlKrif7KEBFRV2Ksp49X+j6FMzPexD8CJ8LG0ATV9XXYFHcSw379DP+5dISvGpFCGLpERAoy0tXDAt+ncHbGm3hr4FhY6BuivK4GX12NwrBfP8PX146horZG7DJJjTF0iYjayETPAK/7jcTZGW9hsd8oGOvqo7imEp9ePIyg3f/Gt/FnUFPPVyqpOYYuEVE7WRgYYcXAMYj+vzfxiu8IGOjooqCqDP+ICcfIvasQnn4NfGyGHsbQJSLqIBtDU7wfOAmnn1mB2V5DoCORIqP0Pl49/hMm/3c9YnLTxS6R1ARDl4iok3Q3scAnQdNwdOoSjO3hAwC4UpCFZ/7YhHlHtyOl6K7IFZLYFA7dkJAQfPjhh4iOjkZ9fb0yayIi0mgelt3wbegc7Bm/EAPsXAAAhzPjEbrvS/zt7F7kVz5+WlTSbgqH7qhRoxAZGYmXXnoJw4YNw4oVKxAREYHKykpl1kdEpLGGOLjhwMS/YsPTz8PVzBr1ggw/3IzBU3s+x8brJ/mwVRfU5sExrl27hoiICERGRiI9PR0GBgYICgpCWFgYRo4cCSsrK2XVqhIcHIOIlKGmvg4/3IzB6itHUfRgBqNe5rb455DJGOXsJXJ11Bq1GZEqNTUVkZGRiIyMRFxcHKRSKQYOHIiwsDCEhobCycmpwwWqGkOXiJSpsLoCqy5HYnviOdQLMgBAqLM3/hE4Cb0sbEWujlqiNqH7sLy8PHkPODY2FvX19fD29sbSpUsRHBzcGV+hEgxdIlKFxMJc/CMmHGfupAJoGNN5vs8ILPYfBVM9A5Gro4epZeg+rKSkBFFRUYiMjMSAAQMwb968zv4KpWHoEpGqCIKA3zPi8OGFg8guKwIAdDMywzsBEzC9l3+zGdxIHGofupqMoUtEqlZZV4tNcSfx9bXjqKqvBQAM7+6OT4ZNRS8LO5GrI054QESkRYx09bDEPxQnpr+BCa59AQBn7qRi9L4v8cXlCFTV1YpcIXUWhi4RkZpwMrXE5lEv4PvRL8LF1Ao1snqsvnIUYfu/wunbKWKXR51ALUM3NzcXK1euxHPPPYcBAwbAy8sLMTExLa575swZzJw5E/3798ewYcPw/vvvo6SkRMUVExF1nlAXb0RNW4q/9guBrkSK9JICzDq8Ba+f+JkDa2g4tQzdjIwMHDx4EMbGxhg6dOhj14uJicErr7wCBwcHbNy4EW+99RaioqLwyiuvQCaTqbBiIqLOZaSrj78PHo9DU/4fArq5AgD2pl3ByL2rsSf1MidS0FC6YhfQkoCAAERHRwMAIiMjERUV1eJ6//nPf9C7d298+eWXkEob/v1gZ2eHl19+GYcOHcKECRNUVjMRkTJ4Wzlgz4SF2JV8ESsv/I6i6gosPrkLB9Ku4pOgaXA0sRC7RGoDtezpNgZoa/Ly8nD9+nVMmTKlyfrDhw+Hvb09Dh8+rMwSiYhURiqR4jnPABybtgzjXX0BAEezExG6dxV+SjrPXq8GUbinu23btjbtWCqVwtzcHJ6envD19W1zYU+SlJQEAOjdu3ezzzw9PZGcnNzp30lEJKZuxmbYPPIFHLx1He+c2497VeV488xvCE+/hn8HTYeLmbXYJdITKBy6n332Wbu+QCKRwMfHBxs2bEC3bt3atY+WFBUVAQAsLJpfWrGwsEB8fHyL2zW+b/U4paWlMDMz63B9RETKIJFIMMmtP4K6u+O9mAPYn3YVp26nIHTfl/hH4CQ87xnAQTXUmMKhu3379jbtWBAElJWV4cqVK9i2bRs+++wzfPHFF20u8Eke95eLf+mISJtZG5pgXchz+FPP/vh79D7kVZbirbO/4WhWAv49/BnYGpmKXSK1QOHQDQwMbNOOy8rKYGpqitDQUNTW1mLv3r1tLq41lpaWAP7X431YcXFxiz1g4MkjijypJ0xEpE7Guvoi0MEN70Tvw4H0aziSlYCL+1bj8+HPIKyHj9jl0SPa9SDVypUrW/28rKwM8+fPl/+5b9++nX6jv/Febkv3bpOSklq810tEpI2sDIyxLuQ5rA2eBXN9Q9yrKsdLR7fjzTO/oby2Wuzy6CHtCt0dO3Zgy5YtLX5WUVGBBQsWICEhQb5s0qRJOH/+fPsqfAwHBwf07dsX4eHhTd7JjY6ORl5eHsaMGdOp30dEpM4kEgmmufsjcsoSBDn0AgD8lHQeY/avwaX8TJGro0btCt3Fixfjiy++wIEDB5osr6qqwsKFCxEXF4c1a9Z0qLBDhw7h0KFDuHz5MgDgwoULOHToEE6cOCFfZ/ny5UhMTMSyZcsQHR2Nffv2YcWKFfDz88O4ceM69P1ERJrI0dQSP4+bj38EToSBji4ySu9h+sGN2HD9BGQCBw0SW7tnGfrnP/+JPXv2YNOmTQgKCkJ1dTUWLlyI2NhYrFmzBqNGjepQYV5eXi0ud3JyajJYxsmTJ7F27VokJibCxMQEo0ePxooVKx57T/dJOMsQEWmLxMJcvHZ8J24W5QEARjp54cvg/4ONIR+yagu1mNpPEAQsWrQI586dwzfffIN169YhJiYGX3zxBcaOHdvhwsTC0CUibVJZV4N/xPwXPyU13OKzNzbH2uBnEdTdXeTKNIdahC4AVFdX48UXX8TVq1chlUrx73//W+OHXmToEpE22p92FW+d/Q1ltdWQSiRY7DcKS/xCoaPACIBdncpD98KFC4/9rKioCG+//TYmTpyIyZMnN/ksICCgwwWqGkOXiLRVekkB/np8J67fywEABDv2xrqQWbAyNBG5MvWm8tD19vZudbCJxl00riMIAiQSSZMnmDUFQ5eItFl1fR1WXvgd2xLOAgCcTS3xzcjZ6GfrJHJl6qszc0GhwTE++eSTDn8RERGJz0BHFx8O/RMGdeuBFWf2ILusCFN/34BPhk3FzN4cHEjZFArdadOmKbsOIiJSoam9/OFl6YD5UT8go/Qelp3ejSsF2fhn4CTo66jlrK9agXfQiYi6qD7WDvh98iKMdvEGAGxPPIf/+2Mz7laUilyZ9mLoEhF1YRYGRtgaOgdvDBgNCSS4mJ+Jyf9dh/j7d8QuTSsxdImIujipRIql/qOxNXQ2jHX1kVNehGkHNyAis+UpUqn9GLpERAQACOvhg30T/wJHEwuU19Xg5aM/YHPcqU6fsKYrY+gSEZGcj7Uj/jtpEQbYuUCAgA8uHMSbZ39DTX2d2KVpBYYuERE10c3YDL+MewVTevkBAHYmXcCciO9QWlMlcmWaT+HQDQkJwYcffojo6GjU19crsyYiIhKZka4evg6ehTcGjAYAnL6Tghl/bOKTzR2kcOiOGjUKkZGReOmllzBs2DCsWLECERERqKysVGZ9REQkEolEgqX+o7FqxAzoSKS4cf8Oph7cgLTifLFL01htnvDg2rVriIiIQGRkJNLT02FgYICgoCCEhYVh5MiRsLKyUlatKsFhIImImovKvomFx3agsq4W1gYm2B72IvztXMQuSyXUZpah1NRUREZGIjIyEnFxcZBKpRg4cCDCwsIQGhoKJyfNG8uToUtE1LJL+ZmYG/EdCqsrYKSrh29GzcbTTp5il6V0ahO6D8vLy5P3gGNjY1FfXw9vb28sXboUwcHBnfEVKsHQJSJ6vLTifPz5yFZklRVCT6qDDU8/j3GuvmKXpVSdmQud9vSyvb09XnjhBXz33Xc4e/YsPvnkEzg5OSE5ObmzvoKIiETWy8IO+ya+Ci9Le9TK6rHw2I84kHZV7LI0Rqf1dLUFe7pERE92v6oczx3eghv370AqkeCL4TPwf70HiV2WUqhlT5eIiLoOa0MT7Bq3AP62LpAJApae/hU7EmPELkvtMXSJiKhdLA2MsXPsPATa9wQAvB29F1vjz4hblJpj6BIRUbuZ6RtiR9jLGN7dHQDwfkw4dtxkj/dxGLpERNQhxnr6+G70i/Lg/dvZfdidclHkqtQTQ5eIiDrMSFcP20LnIqCbKwQIWHZ6N8LTr4ldltppd+jW19cjMzMTeXl5nVkPERFpKGM9fWwPewl+ts6QCQJeP/EzDmfcELsstaLbno0uXryIxYsX4969ewAAExMT9OnTB76+vvDx8YGvry969eoFiUTSqcUSEZF6M9M3xI4xL+PZQ98g/v4dvHr8J/w4dh6GOfQSuzS10K73dKdPn4709HTMnTsX+vr6SE5Oxo0bN5CVlQVBECCRSGBkZIRLly4po2al4nu6REQdd6+qDNN/34TU4nyY6xtiz/i/oI+1g9hltUtn5kK7erppaWmYP38+Fi1a1GR5aWkpbty4gRs3biA+Pr7DxRERkWayMTTFjrCXMOXgBtytLMXsiK04MPGvcDS1FLs0UbXrnq6trS2sra2bLTczM8PQoUMxb948fPHFFx0ujoiINJeLmTV+CHsJpnoGyK0owQsRW1FYXSF2WaJqV+hOmDAB58+f7+xaiIhIy/jaOOLbUbOhJ9VBUtFdzDu6HVV1tWKXJRqFQnfy5Ml488038e233+L06dPye7rh4eHKro+IiDTccEcPfPnUTADA+bxbWHFmD7rqsP8K3dPV09PDoUOHcODAAfkTyfr6+njzzTdx5MgRjB8/Hn5+fho5fy4RESnflF5+yCkvwsexf2Bv2hV4WTlgUf+nxS5L5RR+erm+vh4pKSlISEiQ/0pMTERJSUnDjiQSmJubw8fHR/7a0IQJE5RavDLw6WUiIuUQBAFLT/2K3amXIIEEW0a9gLEaMBevWk1in52djcTERMTHx8vDODc3FxKJBAkJCR0uUNUYukREylNdX4eZf2zGxfxMGOvqY9/EV+Fj3V3sslqlstAtKCiAra1tm3daVFSEhIQEDBs2rEPFiYGhS0SkXPmVpZgY/jVulxejh6k1fv/TIlgaGItd1mOpbD7dsWPHYtu2baivr2/TTi0tLTUycImISPnsjMywLXQuDHR0kVl2H0tP/QqZIBO7LJVoNXRXrVqFX375BZMmTcKZM5wjkYiIOoevjSM+GTYVABCRlYCNcafELUhFWg3dkJAQhIeHY+bMmViyZAkWLVqEnJwcVdVGRERabGbvwZjVu+HS7WcXDyM6N03kipTvie/p6urq4qWXXsLhw4dhZWWFSZMmYc2aNaiurlZFfUREpMU+HDoFPtbdUS/I8NrxnbhbUSp2SUql8IhU1tbW+PDDD/HTTz/hwoULGD9+PA4dOqTM2oiISMsZ6eph08g/w0zPAHcrS7Hs9K9aPXBGm4eBNDExwdy5c2FlZYWlS5cqoyYiIupC3Mxt8e/hzwAAjuck4buEaJErUp5WR6Q6deoUkpOTkZycjKSkJKSlpaGmpgbOzs7w8vLCyJEjVVVni2JiYjBnzpwWP/v999/h7u6u4oqIiKg9Jrv1R2RWAvakXsbK2N8x3NEdnpb2YpfV6VoN3eXLl8PT0xNeXl6YNWsWvLy84OnpCUNDQ1XVp5Dly5cjICCgyTJnZ2eRqiEiovZYOXQKzufdQlZZIRad+Bnhk16DgU67ZqBVW622JiYmRlV1dIibmxv8/f3FLoOIiDrATN8QXwU/ixl/bEL8/Tv4+toxvDEgTOyyOlW7pvYDGsbQTEvT/se7iYhIdQLte2KhbzAAYO3VY0i4nytyRZ2r3aFbW1uLiRMndmYt7fb+++/Dx8cHgwYNwsKFCxEXFyd2SURE1E7LBoyGm7kt6gQZlp/ZjTpZ20ZFVGcdulgu9mPdZmZmmDt3LgIDA2FpaYnU1FRs3rwZzz33HHbs2AE/P79m2zSOofk4paWlMDMzU1bJRET0BEa6evjP8OmY8cdmXC3IxrfxZ7Cwb7DYZXWKVnu6MTExKC19/IvKjXPrisXHxwd///vfMXr0aAwePBjPPvssfv75ZxgbG2P16tWi1kZERO031KEX5noPBQD851IEskrvi1xR52i1p7tgwQLU1tbC2dlZPkdu4y9jY/WcEcLOzg4jRoxAVFRUi58/aZaIJ/WEiYhINf42eDwOZ8Yjt6IEH1w4iG9GzRa7pA5rNXQvXbqElJQU3LhxA/Hx8YiKisKGDRtQVVUFe3v1fX9KJusas1UQEWkzUz0DvBswAYtO/Iw/Mm7g1O1kPOXYW+yyOqTV0NXV1YW3tze8vb3xzDMNo4XIZDKkpaUhLi6u1Unqa2pqoK+v37nVKiA/Px9nz57lK0RERFpgipsffkiMQUxeOt4/F44jUxdDT6ojdlnt1uYHqaRSKTw8PODh4YGpU6c2+zwuLg67d+/GH3/8ofT3fN944w24uLjA19cX5ubmSEtLwzfffIOqqiosW7ZMqd9NRETKJ5FI8OHQyRh3YC2Si+/iu4SzWOD7lNhltVunDPVRVFSEAwcOYPfu3UhOToYgCCp5yMrLywsHDx7Ejh07UFlZCUtLSwQGBuLVV1+Fp6en0r+fiIiUz8faEbO9huD7xHNYdTkSM9wHwsrQROyy2kUidOC9n1OnTmHPnj2IiopCTU0NTExMMHbsWNTU1ODgwYOtXn5WV40PUj3pgSsiIlKdwuoKDN/9b5TUVGFh32C8FzBBZd/dmbnQ5p5uVlYWfvvtN+zbtw+5ubmQSqUYNmwYpk6dirCwMBgYGGDz5s04ePBgh4sjIiICACsDY/y1Xwg+vXgY3yWcxbw+QXA0tRS7rDZTOHTDw8Oxe/duXLhwATKZDJ6enpg9ezYmT54MOzs7ZdZIRESEl/sMx9b4s7hbWYpVVyLx+YgZYpfUZgqH7ooVK2BgYIA5c+ZgypQp6NOnjzLrIiIiasJYTx9L/UPxt+h9+CXlIv7a72n0srAVu6w2UXjsZT09PVRXVyMqKgpRUVHIyclRZl1ERETNzPIMgIupFWSCgA1xJ8Qup80UDt0zZ87g3XffhYmJCdauXYvRo0dj9uzZ2L17N8rKypRZIxEREQBAT6qDV/uFAAB2p1zC7bIicQtqI4VD19zcHC+88AL27t2LPXv2YNasWUhKSsK7776LESNG4I033sDJkyc5GhQRESnVTI9BsDMyRa2sHhvjTopdTpt06JWhmpoaHD58GLt378b58+cBADY2NujWrRsSEhL4yhARESnFhusn8FHsHzDU0UPMzLdgY2iqtO/qzFxo93y6AKCvr4/Jkyfj+++/R0REBP7yl79AT08P8fHxHS6MiIjocWZ7D4WFviGq6mvx080LYpejsA6F7sOcnZ2xePFiREVFYfPmzRgzZkxn7ZqIiKgJUz0DPNu7oQf6w81zGjPRfaeFbiOJRILg4GB89dVXnb1rIiIiudneQyGBBLfLixGRpRm3Mzs9dImIiFTBzdwWI50bxtn/PuGcyNUohqFLREQa68U+QQCA03dSkFx0V+RqnoyhS0REGutpp95wNbMGAOxKVv+3Thi6RESksaQSKWZ4DAQA7E29jHo1HyuCoUtERBptuvsAAEBeZSlO30kRuZrWMXSJiEijuZrZINC+JwBgT+plcYt5AoYuERFpvOm9Gnq7RzLjUV1fJ3I1j8fQJSIijTemhw8kkKCsthpn76SKXc5jMXSJiEjjdTM2w6BuPQAAhzPVdyhihi4REWmFcT18ATRcYpYJ6vkUM0OXiIi0wjjXhtC9W1mKqwU5IlfTMoYuERFphZ7mNuhpZgMAOKOm93UZukREpDWGd3cHALV9mIqhS0REWiPoQeiez7uFGjV8dYihS0REWiOoey8AQFV9LS7nZ4lcTXMMXSIi0hp2RmbwtOwGQD3v6zJ0iYhIqwxzaOjtXszPFLmS5hi6RESkVfxsnQEA1wqyIQiCyNU0xdAlIiKt0s+mIXQLqyuQU14kbjGPYOgSEZFW6W1pB0MdPQDANTUbJIOhS0REWkVXqiN/mCql+K7I1TTF0CUiIq3Ty8IOAJBSnC9yJU0xdImISOt4PAjdVIYuERGRcrk/FLrq9AQzQ5eIiLRO48QHZbXVKKqpFLma/2HoEhGR1rEyNJb/vriaoUtERKQ05vpG8t+XsKdLRESkPKZ6+pBKJACAYoYuERGR8kglUpjpGQIAimuqRK7mfxi6RESklSwNGi4x8/JyJyovL8fKlSsxYsQI9O/fH9OnT8fRo0fFLouIiETWeF+XD1J1okWLFiE8PByLFy/Gpk2b4OHhgUWLFuHEiRNil0ZERCKy0G/s6arP5WVdsQvoiBMnTuDs2bP4+uuvERYWBgAYOnQosrKy8OmnnyIkJETkComISCzm+o33dNnT7RQREREwMzNDaGiofJlEIsG0adOQlpaGlJQUEasjIiIxWTy4p6tOoavRPd3k5GR4eHhAKm36bwcvLy8AQFJSEjw8PJp8Nnjw4Fb3WVpaCjMzs84tlIiIVI73dDtZUVERLCwsmi1vXFZUVKTiioiISF04m1gCAPSlOuIW8hCN7ukCDZeT2/JZbGxsq/t7Uk+YiIg0w8zeg1AvyDDK2VvsUuQ0OnQtLS1b7M0WFxcDQIu9YCIi6hpM9Aww33eE2GU0odGXlz08PJCamgqZTNZkeVJSEgDA09NTjLKIiIhapNGhGxYWhpKSEkRFRTVZvm/fPri5uTV7iIqIiEhMGn15OSQkBEOGDME777yDoqIiODs7Y9++fbh48SLWr18vdnlERERNaHToSiQSrF+/HqtWrcLq1atRUlICDw8PfP311xg1apTY5RERETUhEQRBELsIddL49PKTnnImIqKuoTNzQaPv6RIREWkSjb68rAxlZWUQBIHv6xIREYCGkQpbGxOiLdjTfYRUKm32P7e0tBSlpaUiVSSurtx2oGu3vyu3Heja7Wfbm7ZdIpE0G264vXhPVwFd+T5vV2470LXb35XbDnTt9rPtyms7e7pEREQqwtAlIiJSEYYuERGRijB0iYiIVIShS0REpCIMXSIiIhVh6BIREakI39MlIiJSEfZ0iYiIVIShS0REpCIMXSIiIhXpUrMMRUdHY//+/bh8+TJyc3NhYWGB/v374/XXX4eXlxeAhlmGtm/fjrNnzyItLQ2VlZVwcXHB9OnT8fzzz0NfX1++v5iYGMyZM6fF7/r999/h7u6uknYpSpH2A8Ds2bNx/vz5ZttPmDABq1evbrKsvLwcq1evxqFDh1BSUgIPDw+89tprCA0NVXp72kKRtmdnZ7da97PPPosPPvgAgOYd+0uXLmHdunVISkpCUVERTExM4OnpiXnz5iEkJKTJumfOnMFXX32FxMREmJiYICwsDMuXL4e5uXmT9TTl2CvSdm0+7xU99tp43ivSdlWf910qdHfu3ImioiK8+OKLcHd3R0FBAbZs2YIZM2bghx9+gL+/P27fvo3t27djypQpeOmll2BsbIxz587h888/x/nz57F+/fpm+12+fDkCAgKaLHN2dlZVsxSmSPsb9ezZE5999lmT7a2srJrtc9GiRYiPj8fy5cvh7OyMvXv3YtGiRdi4cWOzH+ZiUqTt3bp1w65du5ptu3fvXvz8888YPXp0s8805diXlJTAzc0N06dPh62tLUpKSrBr1y688sorWLVqFSZOnAig4YfKK6+8gtDQUCxZsgR3797F559/jqSkJPz0009NZlrRlGOvSNu1+bxX9NgD2nfeK9J2lZ/3QhdSUFDQbFlxcbEwePBgYdGiRYIgCEJ5eblQXl7ebL21a9cKnp6eQmJionzZuXPnBE9PTyEiIkJ5RXciRdovCILwwgsvCH/605+euL/jx48Lnp6ewpEjR+TLZDKZMGvWLGHcuHGdU3QnUbTtLZk8ebIQHBws1NfXy5dp2rFvSW1trRAcHCzMnj1bvuyZZ54RpkyZ0qStp0+fFjw9PYWDBw/Kl2nSsW/Jo23X5vO+JS0de20871vSUttboqzzvkvd07WxsWm2zNzcHK6ursjNzQUAGBsbw9jYuNl6/fr1AwD5eppIkfa3RUREBMzMzJpcmpFIJJg2bRrS0tKQkpLSoXo7U3vbfv36ddy8eRPTpk3rtPk01YWuri7MzMygp6cHAMjLy8P169cxZcqUJm0dPnw47O3tcfjwYfkyTTr2LXm07dp83rfk0fa3hbYd+5Yo87zXrp8i7XD//n0kJyejd+/era537tw5SCQSeHh4NPvs/fffh4+PDwYNGoSFCxciLi5OWeV2use1Pz09HQEBAfDx8cGYMWOwfv161NbWNlknOTkZHh4ezf5SNt4jTUpKUm7xHaTIsd+zZw8kEgmeeeaZFj/XtGMvk8lQV1eHvLw8rFmzBrdu3cLcuXMB/O94tfT/w9PTE8nJyfI/a+Kxb63tj6NN570i7dfW876tx16Z532Xuqf7KEEQ8N5770Emk2HevHmPXe/atWv44YcfMGXKFDg5OcmXm5mZYe7cuQgMDISlpSVSU1OxefNmPPfcc9ixYwf8/PxU0Yx2e1z7Bw0ahAkTJqBXr16oqKhAZGQk1qxZgxs3bmDdunXy9YqKitCzZ89m+7WwsJB/rq4UOfbV1dU4ePAgAgMD4eLi0uQzTT32S5YskfdYTU1N8eWXXyI4OBjA/45X4/F7mIWFBeLj4+V/1sRj31rbW6Jt5/2T2q/N531bjr3Sz/t2X5jWAp9++qng6ekp7Nmz57Hr3Lp1SwgKChImT54slJaWPnGfd+/eFQIDA4W5c+d2YqXKoUj7G61atUrw9PQULly4IF82ZswYYeHChc3WTU9PFzw9PYWffvqpU+vtTIq0/cCBA4Knp6ewb98+hfapCcc+MzNTuHr1qnD06FFh6dKlgq+vrxAeHi4Iwv/ae/369WbbLVu2TAgKCpL/WROPfWttf5Q2nvdtaX8jbTnv29J2ZZ/3Xfby8urVq7F161a88847mD59eovrZGVlYc6cOTA3N8e2bdtgamr6xP3a2dlhxIgRuHr1ameX3KkUaf/Dpk6dCgC4cuWKfJmlpWWL/6otLi4G0HKPSR0o2vY9e/bAzMwMY8eOVWi/mnDsXVxc0L9/f4waNQqrVq3CiBEj8MEHH0Amk8HS0hJAyz2V4uLiJsdTE499a21/mLae94q2/2Hact63pe3KPu+7ZOh+9dVX2LhxI1asWPHYd64aTzwDAwN89913LT6I8zit/SVWB4q0/1GNbXr4Po6HhwdSU1Obtbfxno6np2cnVdx5FG17Tk4Ozp07h4kTJ8LQ0FDh/av7sX9Uv379UFxcjPv378vv5T5877ZRUlJSk3u9mnjsH/Vw2xtp83n/qJba/yhtOe8f9bi2q+K873Kh+/XXX2P9+vVYvHgx5s+f3+I6OTk5mDt3LqRSKb7//nvY29srvP/8/HycPXu2yTuv6kSR9rdk//79ANDknkVYWBhKSkoQFRXVZN19+/bBzc2txYdPxNSWtv/2228QBOGxD1K0RN2P/aMEQcD58+dhbm4OS0tLODg4oG/fvggPD2/yQyQ6Ohp5eXkYM2aMfJmmHftHPdp2QLvP+0e11P6WaMN5/6jW2q6K875LPUi1detWrF27FiNHjkRQUFCTSyb6+vrw8fHBvXv3MHfuXNy7dw8ff/wx8vLykJeXJ1+vR48esLa2BgC88cYbcHFxga+vL8zNzZGWloZvvvkGVVVVWLZsmaqb90SKtD82NhabN2/GmDFj4OTkhIqKChw9ehS//fYbxo0bh0GDBsm3CQkJwZAhQ/DOO++gqKgIzs7O2LdvHy5evNjiYAJiUqTtjQRBwN69e+Hp6Yn+/fu3uD9NO/ZvvPEGnJyc4OvrCysrK+Tn52Pv3r04d+4c3nvvPejqNvwoWL58OebNm4dly5bh2WefRV5eHj7//HP4+flh3Lhx8v1p0rFXpO3afN4r0n5tPe8V/XsPqO6871JT+z1umDMAcHJyQlRUVKvDfAHAJ598Ir8PuHnzZhw8eBA5OTmorKyEpaUlAgMD8eqrr6rlJRZF2p+RkYGPPvoIiYmJKCwshFQqhZubG6ZOnYrZs2dDR0enyXZlZWVYtWoVDh8+3GQ4uJZGcRGTIm1vFB0djRdffBF/+9vf8OKLL7a4jaYd+x07diA8PBy3bt1CaWkpzMzM0LdvX/z5z3/GqFGjmqx78uRJrF27Vj4M5OjRo7FixYpm9+o05dgr0nZtPu8Vab+2nvdt+XuvqvO+S4UuERGRmLrcPV0iIiKxMHSJiIhUhKFLRESkIgxdIiIiFWHoEhERqQhDl4iISEUYukRERCrC0CUiIlIRhi6RFikvL0efPn3g5eWl0K/GmWEe51//+heeeuoptDaGTkxMDLy8vPDtt982++z8+fMYNGgQRowYgcTExA63j0jTdamxl4m0XX19PT799NMmy3bu3InLly/jrbfeajJrjr6+fqvTsAmCgKNHjyI0NBQSiaTNtRw7dgyLFy+Gra0ttm3bBldX1zbvg0jbMHSJtIi5uTmmTJnSZNn3338PAwMDzJkzp8kA709y/fp15OXltWs83fDwcLz99tvo2bMntm7d2qYZe4i0GS8vE2mx2tpaJCcnw8vLq02BCwAREREwMzPDkCFD2rTdTz/9hDfffBM+Pj7YsWMHA5foIezpEmmxlJQU1NTUoE+fPm3eNiIiAiEhIdDT01N4m02bNmHVqlUYOnQo1q9fDxMTkzZ/L5E2Y+gSabH4+HgAgK+vb5u2S01NRXp6OpYsWaLwNjt37kRWVhZGjx6N1atXQ19fv03fSdQV8PIykRZrDN229nQjIyOhr6+Pp556SuFt8vPzAQAuLi4MXKLHYE+XSIvFx8dDV1cXXl5ebdouIiICw4cPb9Pl4QULFuDChQvYtm0bAODtt99u03cSdQXs6RJpKZlMhsTERPTq1QsGBgYKb5ebm4u4uDiEhoa26fuMjIywadMmDBs2DNu2bcPHH3/c1pKJtB5Dl0hL3bp1CxUVFW2+nxsZGQmJRNLm0AUAQ0NDbNy4EUFBQfj+++/x0UcftXkfRNqMoUukpdp7PzciIgKDBg2CtbV1u77X0NAQGzZswPDhw7F9+3asXLmyXfsh0kYMXSIt1Z4nl4uLixEbG9uuATEe1hi8I0aMwA8//IAPPvigQ/sj0hYMXSItlZCQAIlEAm9vb4W3OXbsGOrq6tp1aflRBgYGWL9+PUaMGIEff/wR//rXv1odw5moK5AIPAuI6IHXXnsN2dnZ2L9/v9ilEGklvjJERHL+/v6YOXOm2GUQaS32dImIiFSE93SJiIhUhKFLRESkIgxdIiIiFWHoEhERqQhDl4iISEUYukRERCrC0CUiIlKR/w/W4zWtG9dFTQAAAABJRU5ErkJggg==",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(figsize=(7, 4))\n",
"sns.lineplot(x=dia.vapor.temperature / KELVIN, y=enthalpy_of_vaporization, ax=ax);\n",
"ax.set_ylabel(r\"$\\Delta^{LV}h$ / kJ / mol\")\n",
"ax.set_xlabel(r\"$T$ / K\");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A more convenient way is to create a dictionary. The dictionary can be used to build pandas `DataFrame` objects.\n",
"This is a bit less flexible, because the units of the properties are rigid. You can inspect the method signature to check what units are used."
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"\u001b[0;31mSignature:\u001b[0m \u001b[0mdia\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_dict\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mDocstring:\u001b[0m\n",
"Returns the phase diagram as dictionary.\n",
"\n",
"Units\n",
"-----\n",
"temperature : K\n",
"pressure : Pa\n",
"densities : mol / m³\n",
"molar enthalpies : kJ / mol\n",
"molar entropies : kJ / mol / K\n",
"\n",
"Returns\n",
"-------\n",
"dict[str, list[float]]\n",
" Keys: property names. Values: property for each state.\n",
"\n",
"Notes\n",
"-----\n",
"xi: liquid molefraction of component i\n",
"yi: vapor molefraction of component i\n",
"i: component index according to order in parameters.\n",
"\u001b[0;31mType:\u001b[0m builtin_function_or_method\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"dia.to_dict?"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
" \n",
" \n",
" \n",
" | \n",
" density liquid | \n",
" density vapor | \n",
" molar enthalpy vapor | \n",
" temperature | \n",
" molar enthalpy liquid | \n",
" pressure | \n",
" molar entropy liquid | \n",
" molar entropy vapor | \n",
" \n",
" \n",
" \n",
" \n",
" | 0 | \n",
" 14125.988947 | \n",
" 52.208491 | \n",
" 22.140400 | \n",
" 230.000000 | \n",
" 3.376293 | \n",
" 96625.278174 | \n",
" 0.039106 | \n",
" 0.120689 | \n",
" \n",
" \n",
" | 1 | \n",
" 14118.006852 | \n",
" 52.811929 | \n",
" 22.135738 | \n",
" 230.280462 | \n",
" 3.383021 | \n",
" 97830.133956 | \n",
" 0.039135 | \n",
" 0.120569 | \n",
" \n",
" \n",
" | 2 | \n",
" 14110.010220 | \n",
" 53.420767 | \n",
" 22.131064 | \n",
" 230.560924 | \n",
" 3.389761 | \n",
" 99046.729400 | \n",
" 0.039164 | \n",
" 0.120449 | \n",
" \n",
" \n",
" | 3 | \n",
" 14101.999011 | \n",
" 54.035036 | \n",
" 22.126380 | \n",
" 230.841386 | \n",
" 3.396514 | \n",
" 100275.143120 | \n",
" 0.039193 | \n",
" 0.120330 | \n",
" \n",
" \n",
" | 4 | \n",
" 14093.973182 | \n",
" 54.654773 | \n",
" 22.121683 | \n",
" 231.121849 | \n",
" 3.403278 | \n",
" 101515.453964 | \n",
" 0.039221 | \n",
" 0.120211 | \n",
" \n",
" \n",
" \n",
" |