From b2506e404ba2ceeae108fec9357bc4aea7dd208d Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.h.bell@gmail.com>
Date: Wed, 24 Feb 2016 21:44:41 -0700
Subject: [PATCH] Fix some edge cases in critical point calculation

---
 src/Backends/Cubics/CubicBackend.cpp             |  4 ++--
 .../Helmholtz/HelmholtzEOSMixtureBackend.cpp     | 16 +++++++++++++---
 src/Solvers.cpp                                  |  8 ++++++--
 3 files changed, 21 insertions(+), 7 deletions(-)

diff --git a/src/Backends/Cubics/CubicBackend.cpp b/src/Backends/Cubics/CubicBackend.cpp
index 88e19e8e..712bc50f 100644
--- a/src/Backends/Cubics/CubicBackend.cpp
+++ b/src/Backends/Cubics/CubicBackend.cpp
@@ -38,8 +38,8 @@ void CoolProp::AbstractCubicBackend::get_critical_point_search_radii(double &R_d
     // Now we scale them to get the appropriate search radii
     double Tr_GERGlike, rhor_GERGlike;
     get_linear_reducing_parameters(rhor_GERGlike, Tr_GERGlike);
-    R_delta *= rhor_GERGlike/rhomolar_reducing();
-    R_tau *= T_reducing()/Tr_GERGlike;
+    R_delta *= rhor_GERGlike/rhomolar_reducing()*5;
+    R_tau *= T_reducing()/Tr_GERGlike*5;
 }
 
 bool CoolProp::AbstractCubicBackend::get_critical_is_terminated(double &delta, double &tau)
diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp
index d449bea8..0e9cb90f 100644
--- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp
+++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp
@@ -3117,7 +3117,7 @@ CoolProp::CriticalState HelmholtzEOSMixtureBackend::calc_critical_point(double r
     tau_delta[0] = T_reducing()/T0; 
     tau_delta[1] = rho0/rhomolar_reducing(); 
     std::string errstr;
-    x = NDNewtonRaphson_Jacobian(&resid, tau_delta, 1e-11, 100, &errstr);
+    x = NDNewtonRaphson_Jacobian(&resid, tau_delta, 1e-10, 100, &errstr);
     _critical.T = T_reducing()/x[0];
     _critical.rhomolar = x[1]*rhomolar_reducing();
     _critical.p = calc_pressure_nocache(_critical.T, _critical.rhomolar);
@@ -3247,11 +3247,21 @@ public:
                 // have an analytic solution for the derivative
                 R_tau = R_tau_tracer; R_delta = R_delta_tracer;
                 theta = Newton(this, theta_last, 1e-10, 100, errstr);
+                
                 // If the solver takes a U-turn, going in the opposite direction of travel
                 // this is not a good thing, and force a Brent's method solver call to make sure we keep
                 // tracing in the same direction
                 if (std::abs(angle_difference(theta, theta_last)) > M_PI/2.0){
-                    theta = Brent(this, theta_last-M_PI/3.5, theta_last+M_PI/3.5, DBL_EPSILON, 1e-10, 100, errstr);
+                    // We have found at least one critical point and we are now well above the density
+                    // associated with the first critical point
+                    if (N_critical_points > 0 && delta > 2.5*critical_points[0].rhomolar/HEOS.rhomolar_reducing()){
+                        // Stopping the search - probably we have a kink in the L1*=0 contour
+                        // caused by poor low-temperature behavior
+                        continue;
+                    }
+                    else{
+                        theta = Brent(this, theta_last-M_PI/3.5, theta_last+M_PI/3.5, DBL_EPSILON, 1e-10, 100, errstr);
+                    }
                 }
             }
             
@@ -3272,7 +3282,7 @@ public:
             // you have bracketed a critical point, run the full critical point solver to
             // find the critical point and store it
             if (i > 0 && M1*M1_last < 0){
-                double rhomolar = HEOS.rhomolar_reducing()*delta_new, T = HEOS.T_reducing()/tau_new;
+                double rhomolar = HEOS.rhomolar_reducing()*(delta+delta_new)/2.0, T = HEOS.T_reducing()/((tau+tau_new)/2.0);
                 CoolProp::CriticalState crit = HEOS.calc_critical_point(rhomolar, T);
                 critical_points.push_back(crit);
                 N_critical_points++;
diff --git a/src/Solvers.cpp b/src/Solvers.cpp
index cce22a8a..e7493c08 100644
--- a/src/Solvers.cpp
+++ b/src/Solvers.cpp
@@ -75,6 +75,11 @@ std::vector<double> NDNewtonRaphson_Jacobian(FuncWrapperND *f, std::vector<doubl
 
         // Update the guess
         for (std::size_t i = 0; i<x0.size(); i++){ x0[i] += v(i);}
+        
+        // Stop if the solution is not changing by more than numerical precision
+        if (v.cwiseAbs().maxCoeff() < DBL_EPSILON*100){
+            return x0;
+        }
         error = root_sum_square(f0);
         if (iter>maxiter){
             *errstring=std::string("reached maximum number of iterations");
@@ -112,8 +117,7 @@ double Newton(FuncWrapper1DWithDeriv* f, double x0, double ftol, int maxiter, st
 
         x += dx;
 
-        if (std::abs(dx/x) < 10*DBL_EPSILON)
-        {
+        if (std::abs(dx/x) < 1e-11){
             return x;
         }