{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Benchmarking models for fugacity coefficients\n",
    "\n",
    "This notebook benchmarks the models for fugacity coefficients in VolFe where possible.\n",
    "\n",
    "Functions from Sulfur_X (Ding et al., 2022) and EVo (Liggins et al., 2022) have been copied into \"function_from_other_tools.py\" to allow comparison.\n",
    "\n",
    "## Python set-up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import VolFe as vf\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import function_from_other_tools as of"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models for fugacity coefficient of O2\n",
    "\n",
    "options = y_O2, function = y_O2\n",
    "\n",
    "### 'ideal' Treat as ideal gas, y = 1 at all P."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_models = [[\"y_O2\", \"ideal\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "\n",
    "PT = {\"P\":500.} # bar\n",
    "PT[\"T\"]=500. # 'C\n",
    "vf.y_O2(PT,models=my_models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 'Shi92' Shi & Saxena (1992) AmMin 77(9-10):1038-1049\n",
    "\n",
    "Comparison with results from EVo (Liggins et al., 2022) - values match perfectly.\n",
    "\n",
    "Calculated with EVo:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Conditions\n",
    "PT = {\"P\":1000.} # bar\n",
    "PT[\"T\"]=1200. # C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 1.209857826330598, 0, 0, 0, 0, 0, 0, 0, 0, 0)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# EVo\n",
    "species_list = [\"O2\"]\n",
    "of.find_Y(PT['P'], PT[\"T\"]+273.15, species_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.209857826330598"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# VolFe\n",
    "my_models = [[\"y_O2\", \"Shi92\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "vf.y_O2(PT,models=my_models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models for fugacity coefficient of H2\n",
    "\n",
    "options = y_H2, function = y_H2\n",
    "\n",
    "### 'Shaw64' Eq. (4) from Shaw & Wones (1964)\n",
    "\n",
    "Note: used a value of -0.011901 instead of -0.11901 as reported for calculation of C3 to match data in Table 2.\n",
    "\n",
    "Table 3 in Shaw & Wone's contains values for the fugacity coefficient for various P and T.\n",
    "\n",
    "We spot check ten values here (using values for bars not atmospheres) - they all match to three decimal places.\n",
    "\n",
    "y = 1.133 (P = 500 bar, T = 500)\n",
    "\n",
    "y = 1.085 (P = 500 bar, T = 850)\n",
    "\n",
    "y = 1.073 (P = 500 bar, T = 1000)\n",
    "\n",
    "y = 1.481 (P = 1600 bar, T = 500)\n",
    "\n",
    "y = 1.296 (P = 1600 bar, T = 850)\n",
    "\n",
    "y = 1.252 (P = 1600 bar, T = 1000)\n",
    "\n",
    "y = 1.749 (P = 2300 bar, T = 500)\n",
    "\n",
    "y = 1.378 (P = 2300 bar, T = 1000)\n",
    "\n",
    "y = 1.615 (P = 3000 bar, T = 850)\n",
    "\n",
    "y = 1.515 (P = 3000 bar, T = 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y =  1.1327128811340401 (P = 500.0 bars, T = 500.0 )\n",
      "y =  1.0853513685035343 (P = 500.0 bars, T = 850.0 )\n",
      "y =  1.0733936356394882 (P = 500.0 bars, T = 1000.0 )\n",
      "y =  1.4812384375778227 (P = 1600.0 bars, T = 500.0 )\n",
      "y =  1.295939045430643 (P = 1600.0 bars, T = 850.0 )\n",
      "y =  1.2516574282929442 (P = 1600.0 bars, T = 1000.0 )\n",
      "y =  1.7496048740312717 (P = 2300.0 bars, T = 500.0 )\n",
      "y =  1.3780736491913488 (P = 2300.0 bars, T = 1000.0 )\n",
      "y =  1.6148638426007402 (P = 3000.0 bars, T = 850.0 )\n",
      "y =  1.5154237931859629 (P = 3000.0 bars, T = 1000.0 )\n"
     ]
    }
   ],
   "source": [
    "my_models = [[\"y_H2\", \"Shaw64\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "\n",
    "PT = {\"P\":500.} # bar\n",
    "PT[\"T\"]=500. # 'C\n",
    "y = vf.y_H2(PT,models=my_models)\n",
    "print('y = ', y, '(P =',PT['P'],'bars, T =', PT['T'],')')\n",
    "PT[\"T\"]=850. # 'C\n",
    "y = vf.y_H2(PT,models=my_models)\n",
    "print('y = ', y, '(P =',PT['P'],'bars, T =', PT['T'],')')\n",
    "PT[\"T\"]=1000. # 'C\n",
    "y = vf.y_H2(PT,models=my_models)\n",
    "print('y = ', y, '(P =',PT['P'],'bars, T =', PT['T'],')')\n",
    "PT = {\"P\":1600.} # bar\n",
    "PT[\"T\"]=500. # 'C\n",
    "y = vf.y_H2(PT,models=my_models)\n",
    "print('y = ', y, '(P =',PT['P'],'bars, T =', PT['T'],')')\n",
    "PT[\"T\"]=850. # 'C\n",
    "y = vf.y_H2(PT,models=my_models)\n",
    "print('y = ', y, '(P =',PT['P'],'bars, T =', PT['T'],')')\n",
    "PT[\"T\"]=1000. # 'C\n",
    "y = vf.y_H2(PT,models=my_models)\n",
    "print('y = ', y, '(P =',PT['P'],'bars, T =', PT['T'],')')\n",
    "PT = {\"P\":2300.} # bar\n",
    "PT[\"T\"]=500. # 'C\n",
    "y = vf.y_H2(PT,models=my_models)\n",
    "print('y = ', y, '(P =',PT['P'],'bars, T =', PT['T'],')')\n",
    "PT[\"T\"]=1000. # 'C\n",
    "y = vf.y_H2(PT,models=my_models)\n",
    "print('y = ', y, '(P =',PT['P'],'bars, T =', PT['T'],')')\n",
    "PT = {\"P\":3000.} # bar\n",
    "PT[\"T\"]=850. # 'C\n",
    "y = vf.y_H2(PT,models=my_models)\n",
    "print('y = ', y, '(P =',PT['P'],'bars, T =', PT['T'],')')\n",
    "PT[\"T\"]=1000. # 'C\n",
    "y = vf.y_H2(PT,models=my_models)\n",
    "print('y = ', y, '(P =',PT['P'],'bars, T =', PT['T'],')')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 'ideal' Treat as ideal gas, y = 1 at all P."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_models = [[\"y_H2\", \"ideal\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "\n",
    "PT = {\"P\":500.} # bar\n",
    "PT[\"T\"]=500. # 'C\n",
    "vf.y_H2(PT,models=my_models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models for fugacity coefficient of S2\n",
    "\n",
    "options = y_S2, function = y_S2\n",
    "\n",
    "### 'ideal' Treat as ideal gas, y = 1 at all P."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_models = [[\"y_S2\", \"ideal\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "\n",
    "PT = {\"P\":500.} # bar\n",
    "PT[\"T\"]=500. # 'C\n",
    "vf.y_S2(PT,models=my_models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 'Shi92' Shi & Saxena (1992) AmMin 77(9-10):1038-1049\n",
    "\n",
    "Comparison with results from EVo (Liggins et al., 2022) - values match perfectly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "PT = {\"P\":1000.} # bar\n",
    "PT[\"T\"]=1200. # C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 0, 0, 0, 0, 0, 1.1913268980917349, 0, 0, 0, 0)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# EVo\n",
    "species_list = [\"S2\"]\n",
    "of.find_Y(PT['P'], PT[\"T\"]+273.15, species_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.1913268980917349"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# VolFe\n",
    "my_models = [[\"y_S2\", \"Shi92\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "vf.y_S2(PT,models=my_models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models for fugacity coefficient of CO\n",
    "\n",
    "options = y_CO, function = y_CO\n",
    "\n",
    "### 'ideal' Treat as ideal gas, y = 1 at all P."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_models = [[\"y_CO\", \"ideal\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "\n",
    "PT = {\"P\":500.} # bar\n",
    "PT[\"T\"]=500. # 'C\n",
    "vf.y_CO(PT,models=my_models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 'Shi92' Shi & Saxena (1992) AmMin 77(9-10):1038-1049\n",
    "\n",
    "Comparison with results from EVo (Liggins et al., 2022) - values match perfectly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "PT = {\"P\":1000.} # bar\n",
    "PT[\"T\"]=1200. # C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 0, 0, 1.2675079568495566, 0, 0, 0, 0, 0, 0, 0)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# EVo\n",
    "species_list = [\"CO\"]\n",
    "of.find_Y(PT['P'], PT[\"T\"]+273.15, species_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.2675079568495566"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# VolFe\n",
    "my_models = [[\"y_CO\", \"Shi92\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "vf.y_CO(PT,models=my_models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models for fugacity coefficient of H2O\n",
    "\n",
    "options = y_H2O, function = y_H2O\n",
    "\n",
    "### 'ideal' Treat as ideal gas, y = 1 at all P."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_models = [[\"y_H2O\", \"ideal\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "\n",
    "PT = {\"P\":500.} # bar\n",
    "PT[\"T\"]=500. # 'C\n",
    "vf.y_H2O(PT,models=my_models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 'Holland91' Holland & Powell (1991) CMP 109:265-273\n",
    "\n",
    "Comparison with results from EVo (Liggins et al., 2022) - values match perfectly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Conditions\n",
    "PT = {\"P\":1000.} # bar\n",
    "PT[\"T\"]=1200. # C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.9988859662849531, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# EVo\n",
    "species_list = [\"H2O\"]\n",
    "of.find_Y(PT['P'], PT[\"T\"]+273.15, species_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9988859662849531"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# VolFe\n",
    "my_models = [[\"y_H2O\", \"Holland91\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "vf.y_H2O(PT,models=my_models)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.9499443446696176"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Conditions\n",
    "PT = {\"P\":5000.} # bar\n",
    "PT[\"T\"]=1000. # C\n",
    "my_models = [[\"y_H2O\", \"Holland91\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "vf.gas_molar_volume('H2O', PT, my_models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models for fugacity coefficient of CO2\n",
    "\n",
    "options = y_CO2, function = y_CO2\n",
    "\n",
    "### 'Shi92' Shi & Saxena (1992)\n",
    "\n",
    "Nothing available in original paper to benchmark too.\n",
    "\n",
    "### 'Holland91' Holland & Powell (1991) CMP 109:265-273 10.1007/BF00306484\n",
    "\n",
    "Nothing available for Eq. (8,9) and Table 2 or Eq. (4,A1-3) and Table 2.\n",
    "\n",
    "Comparison with results from EVo (Liggins et al., 2022) - values match perfectly for Eq. (8) and Table 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Conditions\n",
    "PT = {\"P\":1000.} # bar\n",
    "PT[\"T\"]=1200. # C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 0, 0, 0, 1.301912410138897, 0, 0, 0, 0, 0, 0)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# EVo\n",
    "species_list = [\"CO2\"]\n",
    "of.find_Y(PT['P'], PT[\"T\"]+273.15, species_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.301912410138897"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# VolFe\n",
    "my_models = [[\"y_CO2\", \"Holland91_eq8_tab1\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "vf.y_CO2(PT,models=my_models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 'ideal' Treat as ideal gas, y = 1 at all P."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_models = [[\"y_CO2\", \"ideal\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "\n",
    "PT = {\"P\":500.} # bar\n",
    "PT[\"T\"]=500. # 'C\n",
    "vf.y_CO2(PT,models=my_models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models for fugacity coefficient of SO2\n",
    "\n",
    "options = y_SO2, function = y_SO2\n",
    "\n",
    "### 'ideal' Treat as ideal gas, y = 1 at all P."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_models = [[\"y_SO2\", \"ideal\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "\n",
    "PT = {\"P\":500.} # bar\n",
    "PT[\"T\"]=500. # 'C\n",
    "vf.y_SO2(PT,models=my_models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 'Shi92' Shi & Saxena (1992) AmMin 77(9-10):1038-1049\n",
    "\n",
    "Comparison with results from Sulfur_X (Ding et al., 2023) - values match perfectly.\n",
    "\n",
    "Calculated using Sulfur_X function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.2064562376797612"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PT = {\"P\":1000.} # bar\n",
    "PT[\"T\"]=1200.+273.15 # K\n",
    "of.phiso2(PT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculated using VolFe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.206456237679761"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_models = [[\"y_SO2\", \"Shi92\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "\n",
    "PT = {\"P\":1000.} # bar\n",
    "PT[\"T\"]=1200. # 'C\n",
    "vf.y_SO2(PT,models=my_models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 'Shi92_Hughes23' Fig.S1 modified from Shi & Saxena (1992) from Hughes et al. (2023) JGSL 180(3)\n",
    "\n",
    "This is the first time this model has been coded up so VolFe represents the benchmark for future use."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models for fugacity coefficient of H2S\n",
    "\n",
    "options = y_H2S, function = y_H2S\n",
    "\n",
    "### 'ideal' Treat as ideal gas, y = 1 at all P."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_models = [[\"y_H2S\", \"ideal\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "\n",
    "PT = {\"P\":500.} # bar\n",
    "PT[\"T\"]=500. # 'C\n",
    "vf.y_H2S(PT,models=my_models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 'Shi92' Shi & Saxena (1992) AmMin 77(9-10):1038-1049\n",
    "\n",
    "Comparison with results from EVo (Liggins et al., 2022) - values match perfectly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "PT = {\"P\":1000.} # bar\n",
    "PT[\"T\"]=1200. # C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 0, 0, 0, 0, 0, 0, 0, 1.489941993252676, 0, 0)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# EVo\n",
    "species_list = [\"H2S\"]\n",
    "of.find_Y(PT['P'], PT[\"T\"]+273.15, species_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.489941993252676"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# VolFe\n",
    "my_models = [[\"y_H2S\", \"Shi92\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "vf.y_H2S(PT,models=my_models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 'Shi92_Hughes23' Fig.S1 modified from Shi & Saxena (1992) Hughes et al. (2024) AmMin 109(3):422-438 doi:10.2138/am-2023-8739\n",
    "\n",
    "This is the first time this model has been coded up so VolFe represents the benchmark for future use."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models for fugacity coefficient of CH4\n",
    "\n",
    "options = y_CH4, function = y_CH4\n",
    "\n",
    "### 'ideal' Treat as ideal gas, y = 1 at all P."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_models = [[\"y_CH4\", \"ideal\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "\n",
    "PT = {\"P\":500.} # bar\n",
    "PT[\"T\"]=500. # 'C\n",
    "vf.y_CH4(PT,models=my_models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 'Shi92' Shi & Saxena (1992) AmMin 77(9-10):1038-1049\n",
    "\n",
    "Comparison with results from EVo (Liggins et al., 2022) - values match perfectly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "# conditions\n",
    "PT = {\"P\":1000.} # bar\n",
    "PT[\"T\"]=1200. # C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 0, 0, 0, 0, 1.2858217235364475, 0, 0, 0, 0, 0)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# EVo\n",
    "species_list = [\"CH4\"]\n",
    "of.find_Y(PT['P'], PT[\"T\"]+273.15, species_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.2858217235364475"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_models = [[\"y_CH4\", \"Shi92\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "vf.y_CH4(PT,models=my_models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models for fugacity coefficient of OCS\n",
    "\n",
    "options = y_OCS, function = y_OCS\n",
    "\n",
    "### 'ideal' Treat as ideal gas, y = 1 at all P."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_models = [[\"y_OCS\", \"ideal\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "\n",
    "PT = {\"P\":500.} # bar\n",
    "PT[\"T\"]=500. # 'C\n",
    "vf.y_OCS(PT,models=my_models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 'Shi92' Shi & Saxena (1992) AmMin 77(9-10):1038-1049\n",
    "\n",
    "Comparison with results from EVo (Liggins et al., 2022) - values match perfectly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "# conditions\n",
    "PT = {\"P\":1000.} # bar\n",
    "PT[\"T\"]=1200 # K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.3560000299582256)"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# EVo\n",
    "species_list = [\"OCS\"]\n",
    "of.find_Y(PT['P'], PT[\"T\"]+273.15, species_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.3560000299582256"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# VolFe\n",
    "my_models = [[\"y_OCS\", \"Shi92\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "vf.y_OCS(PT,models=my_models)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "volfe-dev",
   "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.13.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}