/** This is an attempt to use the water.frink libraries to obtain a
    high-accuracy boiling point of water at a specific ambient atmospheric
    pressure.

    It is primarily based on The International Association for the Properties
    of Water and Steam (IAPWS, http://iapws.org/ ) document IAPWS R6-95(2018),
    a PDF of which is usually available from:

    http://www.iapws.org/relguide/IAPWS-95.html

    Specifically, the file
    http://www.iapws.org/relguide/IAPWS95-2018.pdf
*/

use water.frink
use secant.frink

ambientPressure = 150 psi
println["Ambient pressure = " + format[ambientPressure, "atm", 6]]

Torig = Water.saturatedTemperature[ambientPressure]
println["Original temperature guess = " + format[C[Torig], 1, 6] + " C"]

rhoL = Water.saturatedLiquidDensity[Torig]
println["Original liquid density guess = " + format[rhoL, "g/cm^3", 6]]

rhoV = Water.saturatedVaporDensity[Torig]
println["Original vapor  density guess = " + format[rhoV, "g/cm^3", 6]]

T = Torig
T1 = T * 1.000001

/* This next bit is trying to apply the phase-equilibrium condition
   (Maxwell criterion) to solve for the unknowns.  See top of p. 12 in IAPWS-95
*/
deltaV = rhoV / Water.rhoc
tau = Water.Tc / T
psigmaV = Water.R T rhoV (1 + deltaV Water.dphirdeltaDimensioned[rhoV, T])
println["psigmaV = " + format[psigmaV, "atm", 6]]

deltaL = rhoL / Water.rhoc
tau = Water.Tc / T
psigmaL = Water.R T rhoL (1 + deltaL Water.dphirdeltaDimensioned[rhoL, T])
println["psigmaL = " + format[psigmaL, "atm", 6]]

println["Actual boiling point: " + (format[C[Water.boilingPoint[ambientPressure]], 1, 6] + " C")]