Can't get properties at critical point

[grrhendra]

I think I'm probably missing something obvious, but I can't get the properties at the critical point:

import iapws

fluid = iapws.IAPWS95()

criticalState = fluid(rho=fluid.rhoc,T=fluid.Tc)

Running this raises ZeroDivisionError: 0.0 cannot be raised to a negative power on line 1917 of iapws95.py. Is there another function I should be using to get the critical point properties?

[jjgomera]

Thanks for report, it was a bug, fixed in last commit, try using it now.

[grrhendra]

Thanks, I don't get an error any more! It looks like there might be a related issue in _saturation, though. Minimum working example:

import iapws

fluid = iapws.IAPWS95()

rhol, rhov, Ps = fluid._saturation(fluid.Tc)

The output has no errors, but there are 3 warnings beginning with RuntimeWarning: divide by zero encountered in scalar power DeltaBd = b*Delta**(b-1)*Deltad. If this function worked all the way to the critical point, then one could calculate the saturation temperature given the saturation pressure via Brent's method:

import iapws
from scipy.optimize import brentq

fluid = iapws.IAPWS95()

def Tsat(P): # P in kPa to T in K
	def PsatError(T):
		return fluid._saturation(T)[2] - P

	return brentq(PsatError, 0.01+273.15, fluid.Tc) # Find a T for which 273.16 <= T <= fluid.Tc and PsatError(T) is within floating point precision of zero

With the current version of iapws, the Tsat(P) method gives some warnings regardless of the given P because it tries to calculate the saturation properties at the upper bound fluid.Tc as part of the first iteration. If the final argument in the code above is changed from fluid.Tc to fluid.Tc-1e-10, then Tsat(P) works without warnings, but it might be for users of iapws to not have to make these kinds of corrections.

[jjgomera]

Thanks, applying original fix too in saturation routine to avoid these warnings

[grrhendra]

This fix works perfectly - thanks!