From 16b4806d00aabd5c031452222387902f2b02d1db Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Juan=20Jos=C3=A9=20G=C3=B3mez=20Romera?=
 <jjgomera@gmail.com>
Date: Mon, 9 Dec 2024 21:26:22 +0100
Subject: [PATCH] Apply bug #77 fixing too in saturation routine

---
 iapws/iapws95.py | 14 +++++++++++---
 1 file changed, 11 insertions(+), 3 deletions(-)

diff --git a/iapws/iapws95.py b/iapws/iapws95.py
index d5ed26d..3b120e4 100644
--- a/iapws/iapws95.py
+++ b/iapws/iapws95.py
@@ -180,7 +180,11 @@ def _phird(tau, delta, coef):
         Delta = Tita**2+B*((delta-1)**2)**a
         Deltad = (delta-1)*(A*Tita*2/bt*((delta-1)**2)**(0.5/bt-1)
                             + 2*B*a*((delta-1)**2)**(a-1))
-        DeltaBd = b*Delta**(b-1)*Deltad
+
+        if Delta == 0:
+            DeltaBd = 0
+        else:
+            DeltaBd = b*Delta**(b-1)*Deltad
 
         fird += n*(Delta**b*(F+delta*Fd)+DeltaBd*delta*F)
 
@@ -257,7 +261,12 @@ def _phirt(tau, delta, coef):
         F = exp(-C*(delta-1)**2-D*(tau-1)**2)
         Ft = -2*D*F*(tau-1)
         Delta = Tita**2+B*((delta-1)**2)**a
-        DeltaBt = -2*Tita*b*Delta**(b-1)
+
+        if Delta == 0:
+            DeltaBt = 0
+        else:
+            DeltaBt = -2*Tita*b*Delta**(b-1)
+
         firt += n*delta*(DeltaBt*F+Delta**b*Ft)
 
     return firt
@@ -441,7 +450,6 @@ class MEoS(_fase):
     cp0_cv = None
     gamma0 = None
 
-
     @classmethod
     def from_list(cls, p1name, p1val, p2name, p2val):
         """Speed up method using multiprocessing for multiple point calculation