reverted a changed made to the Takahashi method. small adjustments to…

… the code style.

src/transport/HighPressureGasTransport.cpp

@@ -34,7 +34,7 @@ double takahashi_correction_factor(double Pr, double Tr)
         return 1.0;
     }
 
-    // Data from Table 2 of Takahashi 1975 paper:
+    // Data from Table 2 of Takahashi 1975 paper
     const static double Pr_lookup[17] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0,
         1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 4.0, 5.0};
     const static double DP_Rt_lookup[17] = {1.01, 1.01, 1.01, 1.01, 1.01, 1.01,
@@ -45,36 +45,53 @@ double takahashi_correction_factor(double Pr, double Tr)
     const static double B_ij_lookup[17] = {1.52267, 2.16794, 2.42910, 2.77605,
         2.98256, 3.11384, 3.50264, 3.07773, 3.54744, 3.61216, 3.41882, 3.18415,
         3.37660, 3.27984, 3.39031, 3.23513, 3.13001};
-    const static double C_ij_lookup[17] = {0., 0., 0., 0., 0., 0., 0., 0.141211,
-        0.278407, 0.372683, 0.504894, 0.678469, 0.665702, 0., 0.602907, 0., 0.};
-    const static double E_ij_lookup[17] = {1., 1., 1., 1., 1., 1., 1., 13.45454,
-        14., 10.00900, 8.57519, 10.37483, 11.21674, 1., 6.19043, 1., 1.};
+    const static double C_ij_lookup[17] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.141211,
+        0.278407, 0.372683, 0.504894, 0.678469, 0.665702, 0.0, 0.602907, 0.0, 0.0};
+    const static double E_ij_lookup[17] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 13.45454,
+        14.0, 10.00900, 8.57519, 10.37483, 11.21674, 1.0, 6.19043, 1.0, 1.0};
 
     // Interpolate to obtain the value of the constants (DP)_R, A, B, C, E at
     // the provided value of the reduced pressure (Pr).
-    int Pr_i = 0;
+    int Pr_lower = 0; // Index of the lower bounding value of Pr
+    int Pr_upper = 0; // Index of the upper bounding value of Pr
     double frac = 0.0;
-    double A, B, C,E, DP_Rt;
+    double A, B, C, E, DP_Rt;
 
+    bool found = false; // Flag to indicate if a bounding value of Pr has been found
     for (int j = 1; j < 17; j++){
         if (Pr_lookup[j] > Pr) {
             frac = (Pr - Pr_lookup[j-1])/(Pr_lookup[j] - Pr_lookup[j-1]);
+            found = true;
+            Pr_lower = j-1;
+            Pr_upper = j;
             break;
         }
-        Pr_i++;
     }
+
     // If this loop completes without finding a bounding value of Pr, use
     // the final table value.
-    frac = 1.0;
-
-
-    DP_Rt = DP_Rt_lookup[Pr_i]*(1.0 - frac) + DP_Rt_lookup[Pr_i+1]*frac;
-    A = A_ij_lookup[Pr_i]*(1.0 - frac) + A_ij_lookup[Pr_i+1]*frac;
-    B = B_ij_lookup[Pr_i]*(1.0 - frac) + B_ij_lookup[Pr_i+1]*frac;
-    C = C_ij_lookup[Pr_i]*(1.0 - frac) + C_ij_lookup[Pr_i+1]*frac;
-    E = E_ij_lookup[Pr_i]*(1.0 - frac) + E_ij_lookup[Pr_i+1]*frac;
+    if (!found)
+    {
+        frac = 1.0;
+        Pr_lower = 16;
+        Pr_upper = 16;
+    }
 
-    double DP_R = DP_Rt*(1.0 - A*pow(Tr,-B))*(1.0 - C*pow(Tr,-E));
+    double DP_Rt_lower = DP_Rt_lookup[Pr_lower];
+    double A_lower = A_ij_lookup[Pr_lower];
+    double B_lower = B_ij_lookup[Pr_lower];
+    double C_lower = C_ij_lookup[Pr_lower];
+    double E_lower = E_ij_lookup[Pr_lower];
+    double DP_R_lower = DP_Rt_lower*(1.0 - A_lower/pow(Tr,B_lower))*(1.0 - C_lower/pow(Tr,E_lower));
+
+    double DP_Rt_upper = DP_Rt_lookup[Pr_upper];
+    double A_upper = A_ij_lookup[Pr_upper];
+    double B_upper = B_ij_lookup[Pr_upper];
+    double C_upper = C_ij_lookup[Pr_upper];
+    double E_upper = E_ij_lookup[Pr_upper];
+    double DP_R_upper = DP_Rt_upper*(1.0 - A_upper/pow(Tr,B_upper))*(1.0 - C_upper/pow(Tr,E_upper));
+
+    double DP_R = DP_R_lower*(1.0 - frac) + DP_R_upper*frac;
     return DP_R;
 }