From a9ad75e974b7d1aed642a23f5b58ff7bd27a2d8d Mon Sep 17 00:00:00 2001
From: Jeff Santner <jeffreysantner@gmail.com>
Date: Fri, 18 May 2018 10:06:51 -0700
Subject: [PATCH] [Python] Add get_equivalence_ratio function

---
 interfaces/cython/cantera/test/test_thermo.py | 22 +++++++++
 interfaces/cython/cantera/thermo.pyx          | 47 +++++++++++++++++++
 2 files changed, 69 insertions(+)

diff --git a/interfaces/cython/cantera/test/test_thermo.py b/interfaces/cython/cantera/test/test_thermo.py
index 1dc919b15..01bbc1e6a 100644
--- a/interfaces/cython/cantera/test/test_thermo.py
+++ b/interfaces/cython/cantera/test/test_thermo.py
@@ -262,6 +262,28 @@ class TestThermoPhase(utilities.CanteraTest):
         self.assertNear(gas['CH'].X[0],  31.0/424.0)
         self.assertNear(gas['OH'].X[0],  11.0/212.0)
 
+    def test_get_equivalence_ratio(self):
+        gas = ct.Solution('gri30.xml')
+        for phi in np.linspace(0.5, 2.0, 5):
+            gas.set_equivalence_ratio(phi, 'CH4:0.8, CH3OH:0.2', 'O2:1.0, N2:3.76')
+            self.assertNear(phi, gas.get_equivalence_ratio())  
+        # Check sulfur species
+        sulfur_species = [k for k in ct.Species.listFromFile('nasa_gas.xml') if k.name in ("SO", "SO2")]
+        gas = ct.Solution(thermo='IdealGas', kinetics='GasKinetics',
+                          species=ct.Species.listFromFile('gri30.xml') + sulfur_species)
+        for phi in np.linspace(0.5, 2.0, 5):
+            gas.set_equivalence_ratio(phi, 'CH3:0.5, SO:0.25, OH:0.125, N2:0.125', 'O2:0.5, SO2:0.25, CO2:0.125')
+            self.assertNear(phi, gas.get_equivalence_ratio())
+        gas.X = 'CH4:1, N2:1, CO2:1, H2O:1'
+        self.assertEqual(gas.get_equivalence_ratio(), np.inf)
+        # Check behavior with oxidizers besides O2, and check optional oxidizer arguments
+        gas.set_equivalence_ratio(0.5, 'CH4:0.8, CH3OH:0.2', 'O2:1.0, N2:3.76, NO:0.1')
+        self.assertNear(0.5, gas.get_equivalence_ratio())
+        gas.X = 'CH4:1, O2:2, NO:0.1'
+        self.assertNear(1.0, gas.get_equivalence_ratio(ignore=['NO']))
+        self.assertNear(0.975, gas.get_equivalence_ratio(oxidizers=['O2']))
+        self.assertNear(gas.get_equivalence_ratio(), gas.get_equivalence_ratio(oxidizers=['O2', 'NO']))
+        
     def test_full_report(self):
         report = self.phase.report(threshold=0.0)
         self.assertIn(self.phase.name, report)
diff --git a/interfaces/cython/cantera/thermo.pyx b/interfaces/cython/cantera/thermo.pyx
index 3bbe34674..b80be1737 100644
--- a/interfaces/cython/cantera/thermo.pyx
+++ b/interfaces/cython/cantera/thermo.pyx
@@ -633,6 +633,53 @@ cdef class ThermoPhase(_SolutionBase):
         Xr = phi * Xf + stoichAirFuelRatio * Xo
         self.TPX = None, None, Xr
 
+    def get_equivalence_ratio(self, oxidizers=[], ignore=[]):
+        """
+        Get the composition of a fuel/oxidizer mixture. This gives the
+        equivalence ratio of an unburned mixture. This is not a quantity that is
+        conserved after oxidation. Considers the oxidation of C to CO2, H to H2O
+        and S to SO2. Other elements are assumed not to participate in oxidation
+        (i.e. N ends up as N2).
+
+        :param oxidizers:
+            List of oxidizer species names as strings. Default: with
+            ``oxidizers=[]``, every species that contains O but does not contain
+            H, C, or S is considered to be an oxidizer.
+        :param ignore:
+            List of species names as strings to ignore.
+
+        >>> gas.set_equivalence_ratio(0.5, 'CH3:0.5, CH3OH:.5, N2:0.125', 'O2:0.21, N2:0.79, NO:0.01')
+        >>> gas.get_equivalence_ratio()
+        0.50000000000000011
+        >>> gas.get_equivalence_ratio(['O2'])  # Only consider O2 as the oxidizer instead of O2 and NO
+        0.48809523809523814
+        >>> gas.X = 'CH4:1, O2:2, NO:0.1'
+        >>> gas.get_equivalence_ratio(ignore=['NO'])
+        1.0
+        """
+        if not oxidizers:  # Default behavior, find all possible oxidizers
+            oxidizers = [s.name for s in self.species() if
+                         all(y not in s.composition for y in ['C', 'H', 'S'])]
+        alpha = 0
+        mol_O = 0
+        for k, s in enumerate(self.species()):
+            if s.name in ignore:
+                continue
+            elif s.name in oxidizers:
+                mol_O += s.composition.get('O', 0) * self.X[k]
+            else:
+                nC = s.composition.get('C', 0)
+                nH = s.composition.get('H', 0)
+                nO = s.composition.get('O', 0)
+                nS = s.composition.get('S', 0)
+
+                alpha += (2 * nC + nH / 2 + 2 * nS - nO) * self.X[k]
+
+        if mol_O == 0:
+            return float('inf')
+        else:
+            return alpha / mol_O
+
     def elemental_mass_fraction(self, m):
         r"""
         Get the elemental mass fraction :math:`Z_{\mathrm{mass},m}` of element