From de6e95e839b635f84b73de1bf99e24f4a6cf16c5 Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.h.bell@gmail.com>
Date: Thu, 4 Dec 2014 11:00:06 -0500
Subject: [PATCH] Allow DY solver for single phase to go to 1.5*Tc

Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
---
 src/Backends/Helmholtz/FlashRoutines.cpp | 24 +++++++++++++++++++++---
 1 file changed, 21 insertions(+), 3 deletions(-)

diff --git a/src/Backends/Helmholtz/FlashRoutines.cpp b/src/Backends/Helmholtz/FlashRoutines.cpp
index 14feccb9..9e2ee7d5 100644
--- a/src/Backends/Helmholtz/FlashRoutines.cpp
+++ b/src/Backends/Helmholtz/FlashRoutines.cpp
@@ -4,6 +4,10 @@
 #include "PhaseEnvelopeRoutines.h"
 #include "Configuration.h"
 
+#if defined(ENABLE_CATCH)
+#include "catch.hpp"
+#endif
+
 namespace CoolProp{
 
 template<class T> T g_RachfordRice(const std::vector<T> &z, const std::vector<T> &lnK, T beta)
@@ -722,7 +726,7 @@ void FlashRoutines::PHSU_D_flash(HelmholtzEOSMixtureBackend &HEOS, parameters ot
                 if (value > Sat->keyed_output(other))
                 {
                     solver_resid resid(&HEOS, HEOS._rhomolar, value, other);
-                    HEOS._T = Brent(resid, Sat->keyed_output(iT), HEOS.Tmax()+1, DBL_EPSILON, 1e-12, 100, errstring);
+                    HEOS._T = Brent(resid, Sat->keyed_output(iT), HEOS.Tmax()*1.5, DBL_EPSILON, 1e-12, 100, errstring);
                     HEOS._Q = 10000;
                     HEOS.calc_pressure();
 					// Update the phase flag
@@ -758,7 +762,7 @@ void FlashRoutines::PHSU_D_flash(HelmholtzEOSMixtureBackend &HEOS, parameters ot
                 {
                     solver_resid resid(&HEOS, HEOS._rhomolar, value, other);
                     HEOS._phase = iphase_gas;
-                    HEOS._T = Brent(resid, TVtriple, HEOS.Tmax()+1, DBL_EPSILON, 1e-12, 100, errstring);
+                    HEOS._T = Brent(resid, TVtriple, HEOS.Tmax()*1.5, DBL_EPSILON, 1e-12, 100, errstring);
                     HEOS._Q = 10000;
                     HEOS.calc_pressure();
                 }
@@ -792,7 +796,7 @@ void FlashRoutines::PHSU_D_flash(HelmholtzEOSMixtureBackend &HEOS, parameters ot
                 {
                     solver_resid resid(&HEOS, HEOS._rhomolar, value, other);
                     HEOS._phase = iphase_liquid;
-                    HEOS._T = Brent(resid, TLtriple, HEOS.Tmax()+1, DBL_EPSILON, 1e-12, 100, errstring);
+                    HEOS._T = Brent(resid, TLtriple, HEOS.Tmax()*1.5, DBL_EPSILON, 1e-12, 100, errstring);
                     HEOS._Q = 10000;
                     HEOS.calc_pressure();
                 }
@@ -1387,4 +1391,18 @@ void FlashRoutines::HS_flash(HelmholtzEOSMixtureBackend &HEOS)
     }
 }
 
+/*
+#define ENABLE_CATCH
+#if defined(ENABLE_CATCH)
+
+TEST_SECTION("PD with T very large should yield error","")
+{
+	shared_ptr<HelmholtzEOSMixtureBackend> HEOS(new HelmholtzEOSMixtureBackend("R134a"));
+	double Tc = HEOS->T_critical();
+	HEOS->update(DmassT, 1.1, 1.5*Tc);
+	CHECK_THROWS(
+}
+#endif
+*/
+
 } /* namespace CoolProp */