/**
 * @file WaterPropsIAPWS.cpp
 * Definitions for a class for calculating the equation of state of water
 * from the IAPWS 1995 Formulation based on the steam tables thermodynamic
 * basis (See class \link Cantera::WaterPropsIAPWS WaterPropsIAPWS\endlink).
 */
/*
 * Copyright (2006) Sandia Corporation. Under the terms of
 * Contract DE-AC04-94AL85000 with Sandia Corporation, the
 * U.S. Government retains certain rights in this software.
 */
#include "cantera/thermo/WaterPropsIAPWS.h"
#include "cantera/base/ctexceptions.h"
#include "cantera/base/stringUtils.h"
#include <cmath>
#include <cstdio>
#include <cstdlib>

namespace Cantera
{
/*
 * Critical Point values of water in mks units
 */
//! Critical Temperature value (kelvin)
const doublereal T_c = 647.096;
//! Critical Pressure (Pascals)
static const doublereal P_c = 22.064E6;
//! Value of the Density at the critical point (kg  m-3)
const doublereal Rho_c = 322.;
//! Molecular Weight of water that is consistent with the paper (kg kmol-1)
static const doublereal M_water = 18.015268;

//! Gas constant that is quoted in the paper
/*
 * Note, this is the Rgas value quoted in the paper. For consistency
 * we have to use that value and not the updated value
 *
 * The Ratio of R/M = 0.46151805 kJ kg-1 K-1 , which is Eqn. (6.3) in the paper.
 */
static const doublereal Rgas = 8.314371E3;   // Joules kmol-1 K-1

// Base constructor
WaterPropsIAPWS:: WaterPropsIAPWS() :
    m_phi(0),
    tau(-1.0),
    delta(-1.0),
    iState(-30000)
{
    m_phi = new WaterPropsIAPWSphi();
}

// Copy constructor
/*
 * @param b Object to be copied
 */
WaterPropsIAPWS::WaterPropsIAPWS(const WaterPropsIAPWS& b) :
    m_phi(0),
    tau(b.tau),
    delta(b.delta),
    iState(b.iState)
{
    m_phi = new WaterPropsIAPWSphi();
    m_phi->tdpolycalc(tau, delta);
}

// assignment constructor
/*
 * @param right Object to be copied
 */
WaterPropsIAPWS& WaterPropsIAPWS::operator=(const WaterPropsIAPWS& b)
{
    if (this == &b) {
        return *this;
    }
    tau = b.tau;
    delta = b.delta;
    iState = b.iState;
    m_phi->tdpolycalc(tau, delta);
    return *this;
}

// destructor
WaterPropsIAPWS::~WaterPropsIAPWS()
{
    delete(m_phi);
    m_phi = 0;
}

/*
 * Calculate the dimensionless temp and rho and store internally.
 *
 * @param temperature   input temperature (kelvin)
 *  @param rho          density in kg m-3
 *
 *  this is a private function
 */
void WaterPropsIAPWS::calcDim(doublereal temperature, doublereal rho)
{
    tau = T_c / temperature;
    delta = rho / Rho_c;
    /*
     * Determine the internal state
     */
    if (temperature > T_c) {
        iState = WATER_SUPERCRIT;
    } else {
        if (delta < 1.0) {
            iState = WATER_GAS;
        } else {
            iState = WATER_LIQUID;
        }
    }
}

// Calculate the Helmholtz free energy in mks units of J kmol-1 K-1,
// using the last temperature and density
doublereal  WaterPropsIAPWS::helmholtzFE() const
{
    doublereal retn = m_phi->phi(tau, delta);
    doublereal temperature = T_c/tau;
    doublereal RT = Rgas * temperature;
    return (retn * RT);
}

/*
 * Calculate the pressure (Pascals), using the
 * current internally stored temperature and density
 *  Temperature: kelvin
 *  rho: density in kg m-3
 */
doublereal  WaterPropsIAPWS::pressure() const
{
    doublereal retn = m_phi->pressureM_rhoRT(tau, delta);
    doublereal rho = delta * Rho_c;
    doublereal temperature = T_c / tau;
    return (retn * rho * Rgas * temperature/M_water);
}

/*
 * Calculates the density given the temperature and the pressure,
 * and a guess at the density. Note, below T_c, this is a
 * multivalued function.
 *
 * parameters:
 *    temperature: Kelvin
 *    pressure   : Pressure in Pascals (Newton/m**2)
 *    phase      : guessed phase of water
 *               : -1: no guessed phase
 *    rhoguess   : guessed density of the water
 *               : -1.0 no guessed density
 *
 * If a problem is encountered, a negative 1 is returned.
 */
doublereal WaterPropsIAPWS::density(doublereal temperature, doublereal pressure,
                                    int phase, doublereal rhoguess)
{

    doublereal deltaGuess = 0.0;
    if (rhoguess == -1.0) {
        if (phase != -1) {
            if (temperature > T_c) {
                rhoguess = pressure * M_water / (Rgas * temperature);
            } else {
                if (phase == WATER_GAS || phase == WATER_SUPERCRIT) {
                    rhoguess = pressure * M_water / (Rgas * temperature);
                } else if (phase == WATER_LIQUID) {
                    /*
                     * Provide a guess about the liquid density that is
                     * relatively high -> convergnce from above seems robust.
                     */
                    rhoguess = 1000.;
                } else if (phase == WATER_UNSTABLELIQUID || phase == WATER_UNSTABLEGAS) {
                    throw Cantera::CanteraError("WaterPropsIAPWS::density",
                                                "Unstable Branch finder is untested");
                } else {
                    throw Cantera::CanteraError("WaterPropsIAPWS::density",
                                                "unknown state: " + Cantera::int2str(phase));
                }
            }
        } else {
            /*
             * Assume the Gas phase initial guess, if nothing is
             * specified to the routine
             */
            rhoguess = pressure * M_water / (Rgas * temperature);
        }

    }
    doublereal p_red = pressure * M_water / (Rgas * temperature * Rho_c);
    deltaGuess = rhoguess / Rho_c;
    setState_TR(temperature, rhoguess);
    doublereal delta_retn = m_phi->dfind(p_red, tau, deltaGuess);
    doublereal density_retn;
    if (delta_retn >0.0) {
        delta = delta_retn;

        /*
         * Dimensionalize the density before returning
         */
        density_retn = delta_retn * Rho_c;
        /*
         * Set the internal state -> this may be
         * a duplication. However, let's just be sure.
         */
        setState_TR(temperature, density_retn);


    } else {
        density_retn = -1.0;
    }
    return density_retn;
}

// Calculates the density given the temperature and the pressure,
// and a guess at the density, while not changing the internal state
/*
 *  Note, below T_c, this is a multivalued function.
 *
 * The #density() function calculates the density that is consistent with
 * a particular value of the temperature and pressure. It may therefore be
 * multivalued or potentially there may be no answer from this function. It therefore
 * takes a phase guess and a density guess as optional parameters. If no guesses are
 *
 * supplied to density(), a gas phase guess is assumed. This may or may not be what
 * is wanted. Therefore, density() should usually at least be supplied with a phase
 * guess so that it may manufacture an appropriate density guess.
 * #density() manufactures the initial density guess, nondimensionalizes everything,
 * and then calls #WaterPropsIAPWSphi::dfind(), which does the iterative calculation
 * to find the density condition that matches the desired input pressure.
 *
 *  @param  pressure   : Pressure in Pascals (Newton/m**2)
 *  @param  phase      : guessed phase of water
 *                     : -1: no guessed phase
 *  @param rhoguess    : guessed density of the water
 *                     : -1.0 no guessed density
 *  @return
 *     Returns the density. If an error is encountered in the calculation
 *     the value of -1.0 is returned.
 */
doublereal WaterPropsIAPWS::density_const(doublereal pressure,
        int phase, doublereal rhoguess) const
{
    doublereal temperature = T_c / tau;
    doublereal deltaGuess = 0.0;
    doublereal deltaSave = delta;
    if (rhoguess == -1.0) {
        if (phase != -1) {
            if (temperature > T_c) {
                rhoguess = pressure * M_water / (Rgas * temperature);
            } else {
                if (phase == WATER_GAS || phase == WATER_SUPERCRIT) {
                    rhoguess = pressure * M_water / (Rgas * temperature);
                } else if (phase == WATER_LIQUID) {
                    /*
                     * Provide a guess about the liquid density that is
                     * relatively high -> convergnce from above seems robust.
                     */
                    rhoguess = 1000.;
                } else if (phase == WATER_UNSTABLELIQUID || phase == WATER_UNSTABLEGAS) {
                    throw Cantera::CanteraError("WaterPropsIAPWS::density",
                                                "Unstable Branch finder is untested");
                } else {
                    throw Cantera::CanteraError("WaterPropsIAPWS::density",
                                                "unknown state: " + Cantera::int2str(phase));
                }
            }
        } else {
            /*
             * Assume the Gas phase initial guess, if nothing is
             * specified to the routine
             */
            rhoguess = pressure * M_water / (Rgas * temperature);
        }

    }
    doublereal p_red = pressure * M_water / (Rgas * temperature * Rho_c);
    deltaGuess = rhoguess / Rho_c;

    delta = deltaGuess;
    m_phi->tdpolycalc(tau, delta);
    //  setState_TR(temperature, rhoguess);

    doublereal delta_retn = m_phi->dfind(p_red, tau, deltaGuess);
    doublereal density_retn;
    if (delta_retn > 0.0) {
        delta = delta_retn;

        /*
         * Dimensionalize the density before returning
         */
        density_retn = delta_retn * Rho_c;

    } else {
        density_retn = -1.0;
    }

    delta = deltaSave;
    m_phi->tdpolycalc(tau, delta);
    return density_retn;
}

// Returns the density (kg m-3)
/*
 * The density is an independent variable in the underlying equation of state
 *
 * @return  Returns the density (kg m-3)
 */
doublereal WaterPropsIAPWS::density() const
{
    return (delta * Rho_c);
}

// Returns the temperature (Kelvin)
/*
 * @return  Returns the internally stored temperature
 */
doublereal WaterPropsIAPWS::temperature() const
{
    return (T_c / tau);
}

/*
 * psat_est provides a rough estimate of the saturation
 * pressure given the temperature. This is used as an initial
 * guess for refining the pressure.
 *
 * Input
 *   temperature (kelvin)
 *
 * return:
 *   psat (Pascals)
 */
doublereal WaterPropsIAPWS::psat_est(doublereal temperature) const
{

    static const doublereal A[8] = {
        -7.8889166E0,
        2.5514255E0,
        -6.716169E0,
        33.2239495E0,
        -105.38479E0,
        174.35319E0,
        -148.39348E0,
        48.631602E0
    };
    doublereal ps;
    if (temperature < 314.) {
        doublereal pl = 6.3573118E0 - 8858.843E0 / temperature
                        + 607.56335E0 * pow(temperature, -0.6);
        ps = 0.1 * exp(pl);
    } else {
        doublereal v = temperature / 647.25;
        doublereal w = fabs(1.0-v);
        doublereal b = 0.0;
        for (int i = 0; i < 8; i++) {
            doublereal z = i + 1;
            b += A[i] * pow(w, ((z+1.0)/2.0));
        }
        doublereal q = b / v;
        ps = 22.093*exp(q);
    }
    /*
     * Original correlation was in cgs. Convert to mks
     */
    ps *= 1.0E6;
    return ps;
}

/*
 * Returns the coefficient of isothermal compressibility
 * of temperature and pressure.
 *          kappa = - d (ln V) / dP at constant T.
 */
doublereal WaterPropsIAPWS::isothermalCompressibility() const
{
    doublereal dpdrho_val = dpdrho();
    doublereal dens = delta * Rho_c;
    return (1.0 / (dens * dpdrho_val));
}

// Returns the value of dp / drho at constant T at  the current
// state of the object
/*
 *  units - Joules / kg
 *
 * @return  returns dpdrho
 */
doublereal WaterPropsIAPWS::dpdrho() const
{
    doublereal retn = m_phi->dimdpdrho(tau, delta);
    doublereal temperature = T_c/tau;
    doublereal val = retn * Rgas * temperature / M_water;
    return val;
}

// Returns the isochoric pressure derivative wrt temperature
/*
 *     beta = M / (rho * Rgas) (d (pressure) / dT) at constant rho
 *
 *  Note for ideal gases this is equal to one.
 *
 *    beta = delta (phi0_d() + phiR_d())
 *            - tau delta (phi0_dt() + phiR_dt())
 */
doublereal WaterPropsIAPWS:: coeffPresExp() const
{
    doublereal retn = m_phi->dimdpdT(tau, delta);
    return (retn);
}

// Returns the coefficient of thermal expansion.
/*
 *           alpha = d (ln V) / dT at constant P.
 *
 * @return  Returns the coefficient of thermal expansion
 */
doublereal WaterPropsIAPWS:: coeffThermExp() const
{
    doublereal kappa = isothermalCompressibility();
    doublereal beta = coeffPresExp();
    doublereal dens = delta * Rho_c;
    return (kappa * dens * Rgas * beta / M_water);
}

// Calculate the Gibbs free energy in mks units of J kmol-1 K-1.
// using the last temperature and density
doublereal WaterPropsIAPWS::Gibbs() const
{
    doublereal gRT = m_phi->gibbs_RT();
    doublereal temperature = T_c/tau;
    return (gRT * Rgas * temperature);
}


// Utility routine in the calculation of the saturation pressure
/*
 *  Private routine
 *
 * Calculate the Gibbs free energy in mks units of
 * J kmol-1 K-1.
 *
 * @param temperature    temperature (kelvin)
 * @param pressure       pressure (Pascal)
 * @param densLiq        Output density of liquid
 * @param densGas        output Density of gas
 * @param delGRT         output delGRT
 */
void  WaterPropsIAPWS::
corr(doublereal temperature, doublereal pressure, doublereal& densLiq,
     doublereal& densGas, doublereal& delGRT)
{

    densLiq = density(temperature, pressure, WATER_LIQUID, densLiq);
    if (densLiq <= 0.0) {
        throw Cantera::CanteraError("WaterPropsIAPWS::corr",
                                    "Error occurred trying to find liquid density at (T,P) = "
                                    + Cantera::fp2str(temperature) + "  " + Cantera::fp2str(pressure));
    }
    setState_TR(temperature, densLiq);
    doublereal gibbsLiqRT =  m_phi->gibbs_RT();

    densGas = density(temperature, pressure, WATER_GAS, densGas);
    if (densGas <= 0.0) {
        throw Cantera::CanteraError("WaterPropsIAPWS::corr",
                                    "Error occurred trying to find gas density at (T,P) = "
                                    + Cantera::fp2str(temperature) + "  " + Cantera::fp2str(pressure));
    }
    setState_TR(temperature, densGas);
    doublereal gibbsGasRT = m_phi->gibbs_RT();

    delGRT = gibbsLiqRT - gibbsGasRT;
}

// Utility routine in the calculation of the saturation pressure
/*
 *  Private routine
 *
 * @param temperature    temperature (kelvin)
 * @param pressure       pressure (Pascal)
 * @param densLiq        Output density of liquid
 * @param densGas        output Density of gas
 * @param pcorr          output corrected pressure
 */
void WaterPropsIAPWS::
corr1(doublereal temperature, doublereal pressure, doublereal& densLiq,
      doublereal& densGas, doublereal& pcorr)
{

    densLiq = density(temperature, pressure, WATER_LIQUID, densLiq);
    if (densLiq <= 0.0) {
        throw Cantera::CanteraError("WaterPropsIAPWS::corr1",
                                    "Error occurred trying to find liquid density at (T,P) = "
                                    + Cantera::fp2str(temperature) + "  " + Cantera::fp2str(pressure));
    }
    setState_TR(temperature, densLiq);
    doublereal prL = m_phi->phiR();

    densGas = density(temperature, pressure, WATER_GAS, densGas);
    if (densGas <= 0.0) {
        throw Cantera::CanteraError("WaterPropsIAPWS::corr1",
                                    "Error occurred trying to find gas density at (T,P) = "
                                    + Cantera::fp2str(temperature) + "  " + Cantera::fp2str(pressure));
    }
    setState_TR(temperature, densGas);
    doublereal prG = m_phi->phiR();

    doublereal rhs = (prL - prG) + log(densLiq/densGas);
    rhs /= (1.0/densGas - 1.0/densLiq);

    pcorr = rhs * Rgas * temperature / M_water;
}


// This function returns the saturation pressure given the
// temperature as an input parameter, and sets the internal state to the saturated
// conditions.
/*
 *  Note this function will return the saturation pressure, given the temperature.
 *  It will then set the state of the system to the saturation condition. The input
 *  parameter waterState is used to either specify the liquid state or the
 *  gas state at the desired temperature and saturated pressure.
 *
 *  If the input temperature, T, is above T_c, this routine will set the internal
 *  state to T and the pressure to P_c. Then, return P_c.
 *
 * @param temperature   input temperature (kelvin)
 * @param waterState    integer specifying the water state
 *
 * @return Returns the saturation pressure
 *                units = Pascal
 */
doublereal WaterPropsIAPWS::psat(doublereal temperature, int waterState)
{
    static int method = 1;
    doublereal densLiq = -1.0, densGas = -1.0, delGRT = 0.0;
    doublereal dp, pcorr;
    if (temperature >= T_c) {
        densGas = density(temperature, P_c, WATER_SUPERCRIT);
        setState_TR(temperature, densGas);
        return P_c;
    }
    doublereal p = psat_est(temperature);
    for (int i = 0; i < 30; i++) {
        if (method == 1) {
            corr(temperature, p, densLiq, densGas, delGRT);
            doublereal delV = M_water * (1.0/densLiq - 1.0/densGas);
            dp = - delGRT * Rgas * temperature / delV;
        } else {
            corr1(temperature, p, densLiq, densGas, pcorr);
            dp = pcorr - p;
        }
        p += dp;

        if ((method == 1) && delGRT < 1.0E-8) {
            break;
        } else {
            if (fabs(dp/p) < 1.0E-9) {
                break;
            }
        }
    }
    // Put the fluid in the desired end condition
    if (waterState == WATER_LIQUID) {
        setState_TR(temperature, densLiq);
    } else if (waterState == WATER_GAS) {
        setState_TR(temperature, densGas);
    } else {
        throw Cantera::CanteraError("WaterPropsIAPWS::psat",
                                    "unknown water state input: " + Cantera::int2str(waterState));
    }
    return p;
}

// Returns the Phase State flag for the current state of the object
/*
 * @param checkState If true, this function does a complete check to see where
 *        in parameter space we are
 *
 *  There are three values:
 *     WATER_GAS   below the critical temperature but below the critical density
 *     WATER_LIQUID  below the critical temperature but above the critical density
 *     WATER_SUPERCRIT   above the critical temperature
 */
int WaterPropsIAPWS::phaseState(bool checkState) const
{
    if (checkState) {
        if (tau <= 1.0) {
            iState = WATER_SUPERCRIT;
        } else {
            doublereal T = T_c / tau;
            doublereal rho = delta * Rho_c;
            //doublereal psatTable = psat_est(T);
            doublereal rhoMidAtm = 0.5 * (1.01E5 * M_water / (8314.472 * 373.15) + 1.0E3);
            doublereal rhoMid = Rho_c + (T - T_c) * (Rho_c - rhoMidAtm) / (T_c - 373.15);
            int iStateGuess = WATER_LIQUID;
            if (rho < rhoMid) {
                iStateGuess = WATER_GAS;
            }
            doublereal kappa = isothermalCompressibility();
            if (kappa >= 0.0) {
                iState = iStateGuess;
            } else {
                // When we are here we are between the spinodal curves
                doublereal rhoDel = rho * 1.000001;

                //setState_TR(T, rhoDel);
                doublereal deltaSave = delta;
                doublereal deltaDel = rhoDel / Rho_c;
                delta = deltaDel;
                m_phi->tdpolycalc(tau, deltaDel);

                doublereal kappaDel = isothermalCompressibility();
                doublereal d2rhodp2 = (rhoDel * kappaDel - rho * kappa) / (rhoDel - rho);
                if (d2rhodp2 > 0.0) {
                    iState = WATER_UNSTABLELIQUID;
                } else {
                    iState = WATER_UNSTABLEGAS;
                }
                //setState_TR(T, rho);
                delta = deltaSave;

                m_phi->tdpolycalc(tau, delta);
            }
        }
    }
    return iState;
}

// Return the value of the density at the water spinodal point (on the liquid side)
// for the current temperature.
/*
 * @return returns the density with units of kg m-3
 */
doublereal WaterPropsIAPWS::densSpinodalWater() const
{
    doublereal temperature = T_c/tau;
    doublereal delta_save = delta;
    // return the critical density if we are above or even just a little below
    // the critical temperature. We just don't want to worry about the critical
    // point at this juncture.
    if (temperature >= T_c - 0.001) {
        return Rho_c;
    }
    doublereal p = psat_est(temperature);
    doublereal rho_low = 0.0;
    doublereal rho_high = 1000;

    doublereal densSatLiq = density_const(p, WATER_LIQUID);
    doublereal dens_old = densSatLiq;
    delta = dens_old / Rho_c;
    m_phi->tdpolycalc(tau, delta);
    doublereal dpdrho_old = dpdrho();
    if (dpdrho_old > 0.0) {
        rho_high = std::min(dens_old, rho_high);
    } else {
        rho_low = std::max(rho_low, dens_old);
    }
    doublereal dens_new = densSatLiq* (1.0001);
    delta = dens_new / Rho_c;
    m_phi->tdpolycalc(tau, delta);
    doublereal dpdrho_new = dpdrho();
    if (dpdrho_new > 0.0) {
        rho_high = std::min(dens_new, rho_high);
    } else {
        rho_low = std::max(rho_low, dens_new);
    }
    bool conv = false;

    for (int it = 0; it < 50; it++) {
        doublereal slope = (dpdrho_new - dpdrho_old)/(dens_new - dens_old);
        if (slope >= 0.0) {
            slope = std::max(slope, dpdrho_new *5.0/ dens_new);
        } else {
            slope = -dpdrho_new;
            //slope = MIN(slope, dpdrho_new *5.0 / dens_new);
            // shouldn't be here for liquid spinodal
        }
        doublereal delta_rho =  - dpdrho_new / slope;
        if (delta_rho > 0.0) {
            delta_rho = std::min(delta_rho, dens_new * 0.1);
        } else {
            delta_rho = std::max(delta_rho, - dens_new * 0.1);
        }
        doublereal dens_est =  dens_new + delta_rho;
        if (dens_est < rho_low) {
            dens_est = 0.5 * (rho_low + dens_new);
        }
        if (dens_est > rho_high) {
            dens_est = 0.5 * (rho_high + dens_new);
        }


        dens_old = dens_new;
        dpdrho_old = dpdrho_new;
        dens_new = dens_est;

        delta = dens_new / Rho_c;
        m_phi->tdpolycalc(tau, delta);
        dpdrho_new = dpdrho();
        if (dpdrho_new > 0.0) {
            rho_high = std::min(dens_new, rho_high);
        } else  if (dpdrho_new < 0.0) {
            rho_low = std::max(rho_low, dens_new);
        } else {
            conv = true;
            break;
        }

        if (fabs(dpdrho_new) < 1.0E-5) {
            conv = true;
            break;
        }
    }

    if (!conv) {
        throw Cantera::CanteraError(" WaterPropsIAPWS::densSpinodalWater()",
                                    " convergence failure");
    }
    // Restore the original delta
    delta = delta_save;
    m_phi->tdpolycalc(tau, delta);

    return dens_new;
}

// Return the value of the density at the water spinodal point (on the gas side)
// for the current temperature.
/*
 * @return returns the density with units of kg m-3
 */
doublereal WaterPropsIAPWS::densSpinodalSteam() const
{
    doublereal temperature = T_c/tau;
    doublereal delta_save = delta;
    // return the critical density if we are above or even just a little below
    // the critical temperature. We just don't want to worry about the critical
    // point at this juncture.
    if (temperature >= T_c - 0.001) {
        return Rho_c;
    }
    doublereal p = psat_est(temperature);
    doublereal rho_low = 0.0;
    doublereal rho_high = 1000;

    doublereal densSatGas = density_const(p, WATER_GAS);
    doublereal dens_old = densSatGas;
    delta = dens_old / Rho_c;
    m_phi->tdpolycalc(tau, delta);
    doublereal dpdrho_old = dpdrho();
    if (dpdrho_old < 0.0) {
        rho_high = std::min(dens_old, rho_high);
    } else {
        rho_low = std::max(rho_low, dens_old);
    }
    doublereal dens_new = densSatGas * (0.99);
    delta = dens_new / Rho_c;
    m_phi->tdpolycalc(tau, delta);
    doublereal dpdrho_new = dpdrho();
    if (dpdrho_new < 0.0) {
        rho_high = std::min(dens_new, rho_high);
    } else {
        rho_low = std::max(rho_low, dens_new);
    }
    bool conv = false;

    for (int it = 0; it < 50; it++) {
        doublereal slope = (dpdrho_new - dpdrho_old)/(dens_new - dens_old);
        if (slope >= 0.0) {
            slope = dpdrho_new;
            //slope = MAX(slope, dpdrho_new *5.0/ dens_new);
            // shouldn't be here for gas spinodal
        } else {
            //slope = -dpdrho_new;
            slope = std::min(slope, dpdrho_new *5.0 / dens_new);

        }
        doublereal delta_rho = - dpdrho_new / slope;
        if (delta_rho > 0.0) {
            delta_rho = std::min(delta_rho, dens_new * 0.1);
        } else {
            delta_rho = std::max(delta_rho, - dens_new * 0.1);
        }
        doublereal dens_est =  dens_new + delta_rho;
        if (dens_est < rho_low) {
            dens_est = 0.5 * (rho_low + dens_new);
        }
        if (dens_est > rho_high) {
            dens_est = 0.5 * (rho_high + dens_new);
        }


        dens_old = dens_new;
        dpdrho_old = dpdrho_new;
        dens_new = dens_est;

        delta = dens_new / Rho_c;
        m_phi->tdpolycalc(tau, delta);
        dpdrho_new = dpdrho();
        if (dpdrho_new < 0.0) {
            rho_high = std::min(dens_new, rho_high);
        } else  if (dpdrho_new > 0.0) {
            rho_low = std::max(rho_low, dens_new);
        } else {
            conv = true;
            break;
        }

        if (fabs(dpdrho_new) < 1.0E-5) {
            conv = true;
            break;
        }
    }

    if (!conv) {
        throw Cantera::CanteraError(" WaterPropsIAPWS::densSpinodalSteam()",
                                    " convergence failure");
    }
    // Restore the original delta
    delta = delta_save;
    m_phi->tdpolycalc(tau, delta);

    return dens_new;
}

/*
 * Sets the internal state of the object to the
 * specified temperature and density.
 */
void WaterPropsIAPWS::setState_TR(doublereal temperature, doublereal rho)
{
    calcDim(temperature, rho);
    m_phi->tdpolycalc(tau, delta);
}

/*
 * Calculate the enthalpy in mks units of
 * J kmol-1 K-1.
 */
doublereal WaterPropsIAPWS::enthalpy() const
{
    doublereal temperature = T_c/tau;
    doublereal hRT =  m_phi->enthalpy_RT();
    return (hRT * Rgas * temperature);
}

/*
 * Calculate the internal Energy in mks units of
 * J kmol-1 K-1.
 */
doublereal WaterPropsIAPWS::intEnergy() const
{
    doublereal temperature = T_c / tau;
    doublereal uRT = m_phi->intEnergy_RT();
    return (uRT * Rgas * temperature);
}

/*
 * Calculate the enthalpy in mks units of356
 * J kmol-1 K-1.
 */
doublereal WaterPropsIAPWS::entropy() const
{
    doublereal sR = m_phi->entropy_R();
    return (sR * Rgas);
}

/*
 * Calculate heat capacity at constant volume
 * J kmol-1 K-1.
 */
doublereal WaterPropsIAPWS::cv() const
{
    doublereal cvR = m_phi->cv_R();
    return (cvR * Rgas);
}

// Calculate the constant pressure heat capacity in mks units of J kmol-1 K-1
// at the last temperature and density
doublereal  WaterPropsIAPWS::cp() const
{
    doublereal cpR = m_phi->cp_R();
    return (cpR * Rgas);
}

// Calculate the molar volume (kmol m-3)
// at the last temperature and density
doublereal WaterPropsIAPWS::molarVolume() const
{
    doublereal rho = delta * Rho_c;
    return (M_water / rho);
}

}