Vectorize p_T and bump patch number (#54)

* Vectorize Saturation Function p_T per optimization on T_p

* Update Patch number to 2.2.1

Author: henningjp

IF97.h

@@ -53,7 +53,7 @@ namespace IF97
     }
 
     // CoolProp-IF97 Version Number
-    static const char IF97VERSION [] = "v2.2.0";
+    static const char IF97VERSION [] = "v2.2.1";
     // Setup Water Constants for Trivial Functions and use in Region Classes
     // Constant values from:
     // Revised Release on the IAPWS Industrial Formulation 1997
@@ -75,7 +75,7 @@ namespace IF97
     // IF97 Constants
     const double Tcrit   = 647.096;             // K
     const double Pcrit   = 22.064*p_fact;       // MPa*
-    const double Rhocrit = 322.0;               // kg/m³
+    const double Rhocrit = 322.0;               // kg/m³
     const double Scrit   = 4.41202148223476*R_fact; // kJ*/kg-K (needed for backward eqn. in Region 3(a)(b)
     const double Ttrip   = 273.16;              // K
     const double Ptrip   = 0.000611657*p_fact;  // MPa*   [Change per IAPWS R7-97(2012), p. 7, Eq. 9]
@@ -2418,7 +2418,7 @@ namespace IF97
         //    The equation is rearranged to solve for rho and turned
         //    into functions f(T,P,rho0) and f'(T,P,rho0) for the
         //    Newton-Raphson technique.  Functions for
-        //    dphi/ddelta and d²phi/ddelta² were also required.  These
+        //    dphi/ddelta and d²phi/ddelta² were also required.  These
         //    additional Taylor functions are defined above.
         //
         double f(double T, double p, double rho0) const{
@@ -2624,10 +2624,26 @@ namespace IF97
             if ( ( T < Tmin ) || ( T > Tcrit ) ){
                 throw std::out_of_range("Temperature out of range");
             }
+            // Initial formulas
             const double theta = T/T_star+n[9]/(T/T_star-n[10]);
-            const double AA =      theta*theta + n[1]*theta + n[2];
-            const double BB = n[3]*theta*theta + n[4]*theta + n[5];
-            const double CC = n[6]*theta*theta + n[7]*theta + n[8];
+            // const double AA =      theta*theta + n[1]*theta + n[2];
+            // const double BB = n[3]*theta*theta + n[4]*theta + n[5];
+            // const double CC = n[6]*theta*theta + n[7]*theta + n[8];
+
+            // First optimization
+            double ABC[3];
+            double& AA = ABC[0];
+            double& BB = ABC[1];
+            double& CC = ABC[2];
+            ABC[0] = 1.0; ABC[1] = n[3]; ABC[2] = n[6];
+            for (int j = 1; j < 3; ++j) {
+                for (int i = 0; i < 3; ++i) {
+                    ABC[i] *= theta;
+                }
+                for (int i = 0; i < 3; ++i) {
+                    ABC[i] += n[i * 3 + j];
+                }
+            }
             return p_star*powi(2*CC/(-BB+sqrt(BB*BB-4*AA*CC)), 4);
         };
         double T_p(double p) const{
@@ -4039,7 +4055,7 @@ namespace IF97
         static const Backwards::Region3aS B3aS;
         static const Backwards::Region3bS B3bS;
 
-        // NOTE: Uncertainty in these reverse functions are ±25 mK, as documented in
+        // NOTE: Uncertainty in these reverse functions are ±25 mK, as documented in
         //       IAPWS R7-97(2012) and IAPWS SR3-03(2014). This can result in temperatures
         //       that are very slightly (within 25 mK) on the opposite side of the
         //       saturation curve from S/H values that are very near saturation.  Use of these
@@ -4062,7 +4078,7 @@ namespace IF97
                                           // When limiting to Tsat, this will keep
                                           // temperature in Region 1 or Region 2.
         if ((p < Pcrit) && Clip) {        // If below Pcrit (where Tsat is available),
-            double Tsat = Tsat97(p);      //     Only calculate Tsat ± eps once and 
+            double Tsat = Tsat97(p);      //     Only calculate Tsat ± eps once and 
             tmin = Tsat + eps;            //     set tmin just above and
             tmax = Tsat - eps;            //     tmax just below saturation.
         }
@@ -4385,7 +4401,7 @@ namespace IF97
             return RegionOutput( IF97_HMASS,RegionOutputBackward(Pmax,s,IF97_SMASS,false,NONE),Pmax, NONE);
         else { 
         // Determining H(s) along Tmax is difficult because there is no direct p(T,s) formulation.
-        // This linear combination fit h(s)=a*ln(s)+b/s+c/s²+d is not perfect, but it's close
+        // This linear combination fit h(s)=a*ln(s)+b/s+c/s²+d is not perfect, but it's close
         // and can serve as a limit along that Tmax boundary. Coefficients in HTmaxdata above.
         // There is a better way to do this using Newton-Raphson on Tmax = T(p,s), but it is iterative and slow.
             double ETA = Hmax_n[0]*log(sigma) + Hmax_n[1]/sigma + Hmax_n[2]/powi(sigma,2) +Hmax_n[3];
@@ -4540,14 +4556,14 @@ namespace IF97
     inline double cvmass_Tp(double T, double p){ return RegionOutput( IF97_CVMASS, T, p, NONE); };
     /// Get the speed of sound [m/s] as a function of T [K] and p [Pa]
     inline double speed_sound_Tp(double T, double p){ return RegionOutput( IF97_W, T, p, NONE); };
-    /// Get the [d(rho)/d(p)]T [kg/m³/Pa] as a function of T [K] and p [Pa]
+    /// Get the [d(rho)/d(p)]T [kg/m³/Pa] as a function of T [K] and p [Pa]
     inline double drhodp_Tp(double T, double p){ return RegionOutput( IF97_DRHODP, T, p, NONE); };
 
     // ******************************************************************************** //
     //                            Transport Properties                                  //
     // ******************************************************************************** //
 
-    /// Get the viscosity [Pa-s] as a function of T [K] and Rho [kg/m³]
+    /// Get the viscosity [Pa-s] as a function of T [K] and Rho [kg/m³]
     /// Used for function verification only
     inline double visc_TRho(double T, double rho) {	
         // Since we have density, we don't need to determine the region for viscosity.
@@ -4557,7 +4573,7 @@ namespace IF97
     /// Get the viscosity [Pa-s] as a function of T [K] and p [Pa]
     inline double visc_Tp(double T, double p) { return RegionOutput(IF97_MU, T, p, NONE); };
 
-    /// Get the thermal conductivity [W/m-K] as a function of T [K], p [Pa], and Rho [kg/m³]
+    /// Get the thermal conductivity [W/m-K] as a function of T [K], p [Pa], and Rho [kg/m³]
     /// Used for function verification only
     inline double tcond_TpRho(double T, double p, double rho) {
         // Since we have density, we don't need to determine the region for viscosity.