From 6145f3adcc7aca4d4a56ed05d150a561a855969b Mon Sep 17 00:00:00 2001
From: Yoel <yoelcortes@gmail.com>
Date: Mon, 5 Sep 2022 12:50:54 -0500
Subject: [PATCH] enhancements to setting thermo

---
 thermosteam/_settings.py                      | 20 +++++++++++++++----
 thermosteam/_thermo.py                        |  9 ++++++---
 .../poyinting_correction_factors.py           | 12 +++++++----
 3 files changed, 30 insertions(+), 11 deletions(-)

diff --git a/thermosteam/_settings.py b/thermosteam/_settings.py
index e2ea3760..ff975014 100644
--- a/thermosteam/_settings.py
+++ b/thermosteam/_settings.py
@@ -155,9 +155,14 @@ def electricity_price(self, electricity_price):
         bst.PowerUtility.price = electricity_price
     
     def set_thermo(self, thermo: tmo.Thermo|Iterable[str|tmo.Chemical], 
-                   cache: Optional[bool]=None,
-                   skip_checks: Optional[bool]=False, 
-                   ideal: Optional[bool]=False):
+            mixture: Optional=None,
+            Gamma: Optional=None,
+            Phi: Optional=None,
+            PCF: Optional=None,
+            cache: Optional[bool]=None,
+            skip_checks: Optional[bool]=False, 
+            ideal: Optional[bool]=False,
+        ):
         """
         Set the default :class:`~thermosteam.Thermo` object. If `thermo` is 
         not a :class:`~thermosteam.Thermo` object, an attempt is made to 
@@ -174,10 +179,17 @@ def set_thermo(self, thermo: tmo.Thermo|Iterable[str|tmo.Chemical],
         ideal :
             Whether to use ideal phase equilibrium and mixture property 
             algorithms.
+        Gamma : :class:`~thermosteam.equilibrium.activity_coefficients.ActivityCoefficients` subclass, optional
+            Class for computing activity coefficients.
+        Phi : :class:`~thermosteam.equilibrium.fugacity_coefficients.FugacityCoefficients` subclass, optional
+            Class for computing fugacity coefficients.
+        PCF : :class:`~thermosteam.equilibrium.poyinting_correction_factors.PoyintingCorrectionFactors` subclass, optional
+            Class for computing poynting correction factors.
             
         """
         if not isinstance(thermo, (tmo.Thermo, tmo.IdealThermo)):
-            thermo = tmo.Thermo(thermo, cache=cache, skip_checks=skip_checks)
+            thermo = tmo.Thermo(thermo, mixture=mixture, cache=cache, skip_checks=skip_checks,
+                                Gamma=Gamma, Phi=Phi, PCF=PCF)
         if ideal: thermo = thermo.ideal()
         self._thermo = thermo
     
diff --git a/thermosteam/_thermo.py b/thermosteam/_thermo.py
index e3aed54d..e15efb32 100644
--- a/thermosteam/_thermo.py
+++ b/thermosteam/_thermo.py
@@ -128,11 +128,14 @@ class Thermo:
     __slots__ = ('chemicals', 'mixture', 'Gamma', 'Phi', 'PCF', '_ideal', '_original_thermo') 
     
     def __init__(self, chemicals, mixture=None,
-                 Gamma=eq.DortmundActivityCoefficients,
-                 Phi=eq.IdealFugacityCoefficients,
-                 PCF=eq.MockPoyintingCorrectionFactors,
+                 Gamma=None,
+                 Phi=None,
+                 PCF=None,
                  cache=None,
                  skip_checks=False):
+        if Gamma is None: Gamma = eq.DortmundActivityCoefficients
+        if Phi is None: Phi = eq.IdealFugacityCoefficients
+        if PCF is None: PCF = eq.MockPoyintingCorrectionFactors
         if not isinstance(chemicals, Chemicals): chemicals = Chemicals(chemicals, cache)
         if not mixture:
             mixture = Mixture.from_chemicals(chemicals)
diff --git a/thermosteam/equilibrium/poyinting_correction_factors.py b/thermosteam/equilibrium/poyinting_correction_factors.py
index 53afc2e3..81446194 100644
--- a/thermosteam/equilibrium/poyinting_correction_factors.py
+++ b/thermosteam/equilibrium/poyinting_correction_factors.py
@@ -60,13 +60,15 @@ class MockPoyintingCorrectionFactors(PoyintingCorrectionFactors):
     """
     __slots__ = ('_chemicals',)
     
-    def __call__(self, T, P, Psats):
+    def __call__(self, T, P, Psats=None):
         return 1.
 
 
 @njit(cache=True)
 def ideal_gas_poyinting_correction_factors(T, P, vls, Psats):
-    return np.exp(vls / (R * T) * (P - Psats))
+    dP = P - Psats
+    dP[dP < 0.] = 0.
+    return np.exp(vls / (R * T) * dP)
 
 class IdealGasPoyintingCorrectionFactors(PoyintingCorrectionFactors):
     """Create an IdealGasPoyintingCorrectionFactors object that estimates 
@@ -81,6 +83,8 @@ class IdealGasPoyintingCorrectionFactors(PoyintingCorrectionFactors):
     """
     __slots__ = ('_chemicals',)
 
-    def __call__(self, T, P, Psats):
-        vls = np.array([i.V('l', T, P) for i in self._chemicals])
+    def __call__(self, T, P, Psats=None):
+        vls = np.array([i.V('l', T, P) for i in self._chemicals], dtype=float)
+        if Psats is None:
+            Psats = np.array([i.Psat(T) for i in self._chemicals], dtype=float)
         return ideal_gas_poyinting_correction_factors(T, P, vls, Psats)
\ No newline at end of file