{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a24b4c77",
   "metadata": {},
   "source": [
    "# Bulk Properties\n",
    "\n",
    "`yaeos` is a library based on Equation-of-State calculations. An Equation of\n",
    "State (EOS) is a mathematical model that describes the relationship between\n",
    "thermodynamic properties of a system, typically pressure (P), volume (V), and\n",
    "temperature (T). These equations are fundamental in physics, chemistry, and\n",
    "engineering, especially in thermodynamics and fluid mechanics.\n",
    "\n",
    "EOS models are used to predict phase behavior, thermodynamic properties, and\n",
    "equilibrium conditions of gases, liquids, solids, and mixtures. They are widely\n",
    "applied in areas such as chemical engineering, petroleum engineering, and\n",
    "material science. From the ideal gas law, to cubic equations of state and SAFT,\n",
    "to more complex models, EOS are essential tools for the study and design of\n",
    "processes.\n",
    "\n",
    "The classic way of defining an EOS is though a mathematical equation that\n",
    "express the pressure as a function of the volume and temperature. For example,\n",
    "the Van der Waals EOS is defined as:\n",
    "\n",
    "$$\n",
    "    P = \\frac{RT}{v - b} - \\frac{a}{v^2}\n",
    "$$\n",
    "\n",
    "Where $P$ is the pressure, $v$ is the molar volume, $T$ is the\n",
    "temperature, $R$ is the gas constant, and $a$ and $b$ are parameters that\n",
    "depend on the substance. For mixtures, the $a$ and $b$ parameters are\n",
    "calculated from the pure components parameters using a mixing rule.\n",
    "\n",
    "The modern way of defining an EOS is through a mathematical model that express\n",
    "the residual Helmholtz free energy as a function of the mole number vector,\n",
    "volume and temperature. This approach is more general and allows for a more\n",
    "flexible definition of the EOS in the context of implementing it in code. All\n",
    "thermodynamic properties can be calculated from the residual Helmholtz free\n",
    "energy and its derivatives. `yaeos` is based on this approach.\n",
    "\n",
    "$$\n",
    "    A^{r}(\\vec{n}, V, T) = -\\int_{\\infty}^{V} \\left(P(\\vec{n},V,T) - \\frac{nRT}{V} \\right) dV\n",
    "$$\n",
    "\n",
    "If you want to know how each bulk property is calculated from the residual\n",
    "Helmholtz free energy, please refer to the [Fortran User Documentation](https://ipqa-research.github.io/yaeos/page/index.html)\n",
    "\n",
    "## Example model\n",
    "First, we are going the import the `yaeos` library and instantiate a\n",
    "Peng-Robinson (1976) model with two components. Later in this tutorial you will\n",
    "learn got to instantiate different models, but for now we will use this one as\n",
    "an example.\n",
    "\n",
    "**The important thing to understand is that all the EoS can evaluate\n",
    "the thermodynamic properties explained in this section.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ff68069c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import yaeos\n",
    "\n",
    "\n",
    "# Model parameters\n",
    "Tc = [320, 375]      # critical temperatures [K]\n",
    "Pc = [30, 45]        # critical pressures [bar]\n",
    "w = [0.0123, 0.045]  # acentric factors [-]\n",
    "\n",
    "model = yaeos.PengRobinson76(Tc, Pc, w)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca9c3a11",
   "metadata": {},
   "source": [
    "## Properties function of Volume and Temperature $\\left(V, T \\right)$\n",
    "We can obtain different bulk properties from a $A^r(\\vec{n}, V, T)$ model\n",
    "specifying the number of moles, volume and temperature. But also, derivatives\n",
    "could be obtained if asked for.\n",
    "\n",
    "All not asked derivatives will be set as `None`. The derivatives are returned\n",
    "as a Python dictionary with the keys `dv`, `dt`, and `dn`, etc. respectively.\n",
    "The compositional derivatives are always returned as `Numpy` arrays.\n",
    "\n",
    "\n",
    "The following properties can be calculated:\n",
    "\n",
    "### Pressure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "86c6989c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "56.000067835576885"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = [5.0, 5.0]  # number of moles [mol]\n",
    "V = 1.0         # volume [L]\n",
    "T = 300.0       # temperature [K]\n",
    "\n",
    "model.pressure(n, V, T) # Calculate pressure [bar]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba929595",
   "metadata": {},
   "source": [
    "You can obtain the next derivatives:\n",
    "\n",
    "$$\n",
    "\\left(\\frac{\\partial P}{\\partial V}\\right)_{T, \\vec{n}} \\quad\n",
    "\\left(\\frac{\\partial P}{\\partial T}\\right)_{V, \\vec{n}} \\quad\n",
    "\\left(\\frac{\\partial P}{\\partial n_i}\\right)_{V, T, n_{j \\neq i}}\n",
    "$$\n",
    "\n",
    "Ask for derivatives as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6905a67e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pressure:  56.000067835576885\n",
      "Derivatives:  {'dv': -647.5628198115064, 'dt': 2.900736382506016, 'dn': array([78.64160081, 50.87096316])}\n"
     ]
    }
   ],
   "source": [
    "# Asking for all derivatives\n",
    "pressure, derivatives = model.pressure(n, V, T, dv=True, dt=True, dn=True)\n",
    "\n",
    "# Displaying results\n",
    "print(\"Pressure: \", pressure)\n",
    "print(\"Derivatives: \", derivatives)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b21e6aa9",
   "metadata": {},
   "source": [
    "### Logarithm of fugacity coefficients $ln \\, \\phi(\\vec{n},V,T)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6a1beaa6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.12084851, -1.61086   ])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = [5.0, 5.0]  # number of moles [mol]\n",
    "V = 1.0         # volume [L]\n",
    "T = 300.0       # temperature [K]\n",
    "\n",
    "model.lnphi_vt(n, V, T)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c27c00c",
   "metadata": {},
   "source": [
    "You can obtain the next derivatives:\n",
    "\n",
    "$$\n",
    "\\left(\\frac{\\partial ln \\, \\hat{\\phi}_i}{\\partial T}\\right)_{P, \\vec{n}} \\quad\n",
    "\\left(\\frac{\\partial ln \\, \\hat{\\phi}_i}{\\partial P}\\right)_{T, \\vec{n}} \\quad\n",
    "\\left(\\frac{\\partial ln \\, \\hat{\\phi}_i}{\\partial n_j}\\right)_{P, T, n_{k \\neq j}}\n",
    "$$\n",
    "\n",
    "\n",
    "Ask for derivatives as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "bf18174c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Logarithm of Fugacity Coefficients:  [-1.12084851 -1.61086   ]\n",
      "Derivatives:  {'dt': array([0.01380036, 0.01926589]), 'dp': array([-0.0129884 , -0.01470769]), 'dn': array([[-0.00904545,  0.00904545],\n",
      "       [ 0.00904545, -0.00904545]])}\n"
     ]
    }
   ],
   "source": [
    "# Asking for all the derivatives\n",
    "lnphi, derivatives = model.lnphi_vt(n, V, T, dt=True, dp=True, dn=True)\n",
    "\n",
    "# Displaying results\n",
    "print(\"Logarithm of Fugacity Coefficients: \", lnphi)\n",
    "print(\"Derivatives: \", derivatives)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6399604f",
   "metadata": {},
   "source": [
    "You can evaluate fugacity by simply:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f38aea95",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fugacities:  [9.12809679 5.59204872]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "p = model.pressure(n, V, T)  # Calculate pressure [bar]\n",
    "\n",
    "z = n / np.sum(n)            # Calculate mole fractions\n",
    "\n",
    "fugacities = np.exp(lnphi) * z * p\n",
    "\n",
    "print(\"Fugacities: \", fugacities)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66dafc71",
   "metadata": {},
   "source": [
    "### Residual Helmholtz free energy $A^r(\\vec{n},V,T)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "06994b8f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-519.8710586930237"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = [5.0, 5.0]  # number of moles [mol]\n",
    "V = 1.0         # volume [L]\n",
    "T = 300.0       # temperature [K]\n",
    "\n",
    "model.helmholtz_residual_vt(n, V, T) # Calculate residual Helmholtz energy [bar L]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c11d453c",
   "metadata": {},
   "source": [
    "You can obtain the next derivatives:\n",
    "\n",
    "$$\n",
    "\\left(\\frac{\\partial A^{r}_{(\\vec{n}, V, T)}}{\\partial V}\\right)_{T, \\vec{n}} \\quad\n",
    "\\left(\\frac{\\partial A^{r}_{(\\vec{n}, V, T)}}{\\partial T}\\right)_{V, \\vec{n}} \\quad\n",
    "\\left(\\frac{\\partial A^{r}_{(\\vec{n}, V, T)}}{\\partial n_i}\\right)_{V, T, n_{j \\neq i}} \\quad\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\left(\\frac{\\partial^2 A^{r}_{(\\vec{n}, V, T)}}{\\partial V \\partial T}\\right)_{ \\vec{n}} \\quad\n",
    "\\left(\\frac{\\partial^2 A^{r}_{(\\vec{n}, V, T)}}{\\partial V^2}\\right)_{T, \\vec{n}} \\quad\n",
    "\\left(\\frac{\\partial^2 A^{r}_{(\\vec{n}, V, T)}}{\\partial T^2}\\right)_{V, \\vec{n}} \\quad\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\left(\\frac{\\partial^2 A^{r}_{(\\vec{n}, V, T)}}{\\partial n_i \\partial n_j}\\right)_{V, T, n_{k \\neq i,j}} \\quad\n",
    "\\left(\\frac{\\partial^2 A^{r}_{(\\vec{n}, V, T)}}{\\partial V \\partial n_i}\\right)_{T, n_{j \\neq i}} \\quad\n",
    "\\left(\\frac{\\partial^2 A^{r}_{(\\vec{n}, V, T)}}{\\partial T \\partial n_i}\\right)_{V, n_{j \\neq i}} \\quad\n",
    "$$\n",
    "\n",
    "Ask for derivatives as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7d79134e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Residual Helmholtz free energy:  -519.8710586930237\n",
      "Derivatives:  {'dt': 1.746237812013294, 'dv': 193.43381070442314, 'dn': array([-65.21921351, -77.44176037]), 'dtv': -2.069290120706016, 'dv2': 398.1289412715064, 'dt2': -0.0021897892913883494, 'dvn': array([-53.69821295, -25.9275753 ]), 'dtn': array([0.39595811, 0.36714748]), 'dn2': array([[6.8304639 , 3.90917869],\n",
      "       [3.90917869, 1.27633637]])}\n"
     ]
    }
   ],
   "source": [
    "# Asking for all the derivatives\n",
    "Ar , derivatives = model.helmholtz_residual_vt(\n",
    "    n, V, T, dv=True, dt=True, dn=True, dtv=True, dvn=True, dtn=True, dn2=True,\n",
    "    dv2=True, dt2=True\n",
    ")\n",
    "\n",
    "# Displaying results\n",
    "print(\"Residual Helmholtz free energy: \", Ar)\n",
    "print(\"Derivatives: \", derivatives)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45d01fc4",
   "metadata": {},
   "source": [
    "### Residual Gibbs free energy $G^r(\\vec{n}, V, T)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "de21a88b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-713.3048693974467"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = [5.0, 5.0]  # number of moles [mol]\n",
    "V = 1.0         # volume [L]\n",
    "T = 300.0       # temperature [K]\n",
    "\n",
    "model.gibbs_residual_vt(n, V, T) # Calculate residual Gibbs energy [bar L]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42bea797",
   "metadata": {},
   "source": [
    "You can obtain the next derivatives:\n",
    "\n",
    "$$\n",
    "\\left(\\frac{\\partial G^{r}_{(\\vec{n}, V, T)}}{\\partial T}\\right)_{V, \\vec{n}} \\quad\n",
    "\\left(\\frac{\\partial G^{r}_{(\\vec{n}, V, T)}}{\\partial V}\\right)_{T, \\vec{n}} \\quad\n",
    "\\left(\\frac{\\partial G^{r}_{(\\vec{n}, V, T)}}{\\partial n_j}\\right)_{V, T, n_{i \\neq j}}\n",
    "$$\n",
    "\n",
    "Ask for derivatives as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "1376ab42",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Residual Gibbs Energy:  -713.3048693974467\n",
      "Derivatives:  {'dt': 3.81552793271931, 'dv': -398.1289412715064, 'dn': array([-11.52100055, -51.51418507])}\n"
     ]
    }
   ],
   "source": [
    "# Asking for all the derivatives\n",
    "Gr, derivatives = model.gibbs_residual_vt(n, V, T, dt=True, dv=True, dn=True)\n",
    "\n",
    "# Displaying results\n",
    "print(\"Residual Gibbs Energy: \", Gr)\n",
    "print(\"Derivatives: \", derivatives)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "003e87ae",
   "metadata": {},
   "source": [
    "### Residual enthalpy $H^r(\\vec{n}, V, T)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "0e3b5062",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1237.176213001435"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = [5.0, 5.0]  # number of moles [mol]\n",
    "V = 1.0         # volume [L]\n",
    "T = 300.0       # temperature [K]\n",
    "\n",
    "model.enthalpy_residual_vt(n, V, T) # Calculate residual enthalpy [bar L]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afb48efb",
   "metadata": {},
   "source": [
    "You can obtain the next derivatives:\n",
    "\n",
    "$$\n",
    "\\left(\\frac{\\partial H^{r}_{(\\vec{n}, V, T)}}{\\partial T}\\right)_{V, \\vec{n}} \\quad\n",
    "\\left(\\frac{\\partial H^{r}_{(\\vec{n}, V, T)}}{\\partial V}\\right)_{T, \\vec{n}} \\quad\n",
    "\\left(\\frac{\\partial H^{r}_{(\\vec{n}, V, T)}}{\\partial n_j}\\right)_{V, T, n_{i \\neq j}}\n",
    "$$\n",
    "\n",
    "Ask for derivatives as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f15b0810",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Residual Enthalpy:  -1237.176213001435\n",
      "Derivatives:  {'dt': 2.7262269081225208, 'dv': 222.65809494029838, 'dn': array([-130.30843347, -161.65842812])}\n"
     ]
    }
   ],
   "source": [
    "# Asking for all the derivatives\n",
    "Hr, derivatives = model.enthalpy_residual_vt(n, V, T, dt=True, dv=True, dn=True)\n",
    "\n",
    "# Displaying results\n",
    "print(\"Residual Enthalpy: \", Hr)\n",
    "print(\"Derivatives: \", derivatives)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b643141",
   "metadata": {},
   "source": [
    "### Residual internal energy $U^r(\\vec{n}, V, T)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "33656746",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1043.7424022970117"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = [5.0, 5.0]  # number of moles [mol]\n",
    "V = 1.0         # volume [L]\n",
    "T = 300.0       # temperature [K]\n",
    "\n",
    "model.internal_energy_residual_vt(n, V, T) # Calculate residual internal energy [bar L]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e736d9aa",
   "metadata": {},
   "source": [
    "You can obtain the next derivatives:\n",
    "\n",
    "$$\n",
    "\\left(\\frac{\\partial U^{r}_{(\\vec{n}, V, T)}}{\\partial V}\\right)_{T, \\vec{n}} \\quad\n",
    "\\left(\\frac{\\partial U^{r}_{(\\vec{n}, V, T)}}{\\partial T}\\right)_{V, \\vec{n}} \\quad\n",
    "\\left(\\frac{\\partial U^{r}_{(\\vec{n}, V, T)}}{\\partial n_i}\\right)_{V, T, n_{j \\neq i}}\n",
    "$$\n",
    "\n",
    "Ask for derivatives as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "7a7832c9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Residual Internal Energy:  -1043.7424022970117\n",
      "Derivatives:  {'dt': 0.6569367874165049, 'dv': 814.2208469162279, 'dn': array([-184.00664642, -187.58600342])}\n"
     ]
    }
   ],
   "source": [
    "# Asking for all the derivatives\n",
    "Ur, derivatives = model.internal_energy_residual_vt(n, V, T, dt=True, dv=True, dn=True)\n",
    "\n",
    "# Displaying results\n",
    "print(\"Residual Internal Energy: \", Ur)\n",
    "print(\"Derivatives: \", derivatives)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c67a3d0d",
   "metadata": {},
   "source": [
    "### Residual entropy $S^r(\\vec{n}, V, T)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "5fa1d37f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.746237812013294"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = [5.0, 5.0]  # number of moles [mol]\n",
    "V = 1.0         # volume [L]\n",
    "T = 300.0       # temperature [K]\n",
    "\n",
    "model.entropy_residual_vt(n, V, T) # Calculate residual entropy [bar L / K]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30224622",
   "metadata": {},
   "source": [
    "You can obtain the next derivatives:\n",
    "\n",
    "$$\n",
    "\\left(\\frac{\\partial S^{r}_{(\\vec{n}, V, T)}}{\\partial T}\\right)_{V, \\vec{n}} \\quad\n",
    "\\left(\\frac{\\partial S^{r}_{(\\vec{n}, V, T)}}{\\partial V}\\right)_{T, \\vec{n}} \\quad\n",
    "\\left(\\frac{\\partial S^{r}_{(\\vec{n}, V, T)}}{\\partial n_j}\\right)_{V, T, n_{i \\neq j}}\n",
    "$$\n",
    "\n",
    "Ask for derivatives as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "d41409c5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Residual Entropy:  -1.746237812013294\n",
      "Derivatives:  {'dt': 0.0021897892913883494, 'dv': 2.069290120706016, 'dn': array([-0.39595811, -0.36714748])}\n"
     ]
    }
   ],
   "source": [
    "# Asking for all the derivatives\n",
    "Sr, derivatives = model.entropy_residual_vt(n, V, T, dt=True, dv=True, dn=True)\n",
    "\n",
    "# Displaying results\n",
    "print(\"Residual Entropy: \", Sr)\n",
    "print(\"Derivatives: \", derivatives)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "280a119b",
   "metadata": {},
   "source": [
    "### Residual heat capacity at constant volume $C_V^r(\\vec{n}, V, T)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "e7cc6b25",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6569367874165049"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = [5.0, 5.0]  # number of moles [mol]\n",
    "V = 1.0         # volume [L]\n",
    "T = 300.0       # temperature [K]\n",
    "\n",
    "model.cv_residual_vt(n, V, T) # Calculate residual heat capacity [bar L / K]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebd4fa3b",
   "metadata": {},
   "source": [
    "### Residual heat capacity at constant pressure $C_P^r(\\vec{n}, V, T)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "15d3ba36",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.72361653132665"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = [5.0, 5.0]  # number of moles [mol]\n",
    "V = 1.0         # volume [L]\n",
    "T = 300.0       # temperature [K]\n",
    "\n",
    "model.cp_residual_vt(n, V, T) # Calculate residual heat capacity [bar L / K]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce69f064",
   "metadata": {},
   "source": [
    "## Properties function of Pressure and Temperature $\\left(P, T \\right)$\n",
    "We can obtain different bulk properties from a $A^r$ model specifying the\n",
    "number of moles, pressure and temperature. But also, derivatives could be \n",
    "obtained if asked for.\n",
    "\n",
    "All not asked derivatives will be set as `None`. The derivatives are returned\n",
    "as a Python dictionary with the keys `dp`, `dt`, and `dn`, etc. respectively.\n",
    "The compositional derivatives are always returned as `Numpy` arrays.\n",
    "\n",
    "When we are calculating properties with pressure as specified variable, we\n",
    "have to ask for a specific volume root. The possible roots are:\n",
    "\n",
    "- \"liquid\"\n",
    "- \"vapor\"\n",
    "- \"stable\" (default if not specified)\n",
    "\n",
    "The \"stable\" root refers to the stable phase of the system. If this options is\n",
    "set, both liquid and vapor roots will be calculated and the phase with the\n",
    "lower residual Gibbs free energy $(G^{r}_{(\\vec{n},P,T)})$ will be returned as\n",
    "the stable phase.\n",
    "\n",
    "The following properties can be calculated:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18d3967e",
   "metadata": {},
   "source": [
    "### Volume"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "d74a70ff",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Volume (liquid) in liters:  1.1546008756262833\n",
      "Volume (vapor) in liters:  245.62683294458913\n",
      "Volume (stable phase) in liters:  245.62683294458913\n"
     ]
    }
   ],
   "source": [
    "n = [5.0, 5.0]  # number of moles [mol]\n",
    "P = 1.0         # pressure [bar]\n",
    "T = 300.0       # temperature [K]\n",
    "\n",
    "volume_liquid = model.volume(n, P, T, root=\"liquid\")\n",
    "volume_vapor = model.volume(n, P, T, root=\"vapor\")\n",
    "volume_stable = model.volume(n, P, T, root=\"stable\")\n",
    "\n",
    "# Displaying results\n",
    "print(\"Volume (liquid) in liters: \", volume_liquid)\n",
    "print(\"Volume (vapor) in liters: \", volume_vapor)\n",
    "print(\"Volume (stable phase) in liters: \", volume_stable)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b802cff",
   "metadata": {},
   "source": [
    "### Logarithm of fugacity coefficients $ln \\, \\phi(\\vec{n},P,T)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "a8922cde",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Logarithm of fugacity coefficients (liquid):  [2.6074475  2.24385923]\n",
      "Logarithm of fugacity coefficients (vapor):  [-0.01507129 -0.01531375]\n",
      "Logarithm of fugacity coefficients (stable phase):  [-0.01507129 -0.01531375]\n"
     ]
    }
   ],
   "source": [
    "n = [5.0, 5.0]  # number of moles [mol]\n",
    "P = 1.0         # pressure [bar]\n",
    "T = 300.0       # temperature [K]\n",
    "\n",
    "ln_phi_l = model.lnphi_pt(n, P, T, root=\"liquid\")\n",
    "ln_phi_v = model.lnphi_pt(n, P, T, root=\"vapor\")\n",
    "ln_phi_s = model.lnphi_pt(n, P, T, root=\"stable\")\n",
    "\n",
    "# Displaying results\n",
    "print(\"Logarithm of fugacity coefficients (liquid): \", ln_phi_l)\n",
    "print(\"Logarithm of fugacity coefficients (vapor): \", ln_phi_v)\n",
    "print(\"Logarithm of fugacity coefficients (stable phase): \", ln_phi_s)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2a4e93c",
   "metadata": {},
   "source": [
    "You can obtain the next derivatives **(on the specified volume root)**:\n",
    "\n",
    "$$\n",
    "\\left(\\frac{\\partial ln \\, \\hat{\\phi}_i}{\\partial T}\\right)_{P, \\vec{n}} \\quad\n",
    "\\left(\\frac{\\partial ln \\, \\hat{\\phi}_i}{\\partial P}\\right)_{T, \\vec{n}} \\quad\n",
    "\\left(\\frac{\\partial ln \\, \\hat{\\phi}_i}{\\partial n_j}\\right)_{P, T, n_{k \\neq j}}\n",
    "$$\n",
    "\n",
    "\n",
    "Ask for derivatives as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "78841079",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Logarithm of fugacity coefficients (stable phase):  [-0.01507129 -0.01531375]\n",
      "Derivatives:  {'dt': array([0.00013302, 0.00013235]), 'dp': array([-0.01513287, -0.01539262]), 'dn': array([[-8.47067747e-08,  8.47067748e-08],\n",
      "       [ 8.47067748e-08, -8.47067747e-08]])}\n"
     ]
    }
   ],
   "source": [
    "# Asking for all the derivatives\n",
    "ln_phi_s, derivatives = model.lnphi_pt(n, P, T, root=\"stable\", dt=True, dp=True, dn=True)\n",
    "\n",
    "# Displaying results\n",
    "print(\"Logarithm of fugacity coefficients (stable phase): \", ln_phi_s)\n",
    "print(\"Derivatives: \", derivatives)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "yaeos",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}