Externalise comparison of lambertW functions

src/compute_state.cc

+#include <cassert>
+
 #include "LambertW.h"
 
 #include <gsl/gsl_sf_lambert.h>
@@ -51,3 +53,23 @@ double compute_state_update_lambert_gsl(double h, double unorm, double L,
   double const rhs = unorm / L - old_state;
   return gsl_sf_lambert_W0(h * std::exp(rhs)) - rhs;
 }
+
+double compute_state_update(double h, double unorm, double L,
+                            double old_state) {
+  double ret1 = compute_state_update_bisection(h, unorm, L, old_state);
+  assert(std::abs(1.0 / h * ret1 - (old_state - unorm / L) / h -
+                  std::exp(-ret1)) < 1e-10);
+
+  double ret2 = compute_state_update_lambert(h, unorm, L, old_state);
+  assert(std::abs(1.0 / h * ret2 - (old_state - unorm / L) / h -
+                  std::exp(-ret2)) < 1e-10);
+
+  double ret3 = compute_state_update_lambert_gsl(h, unorm, L, old_state);
+  assert(std::abs(1.0 / h * ret3 - (old_state - unorm / L) / h -
+                  std::exp(-ret3)) < 1e-10);
+
+  assert(std::abs(ret1 - ret2) < 1e-14);
+  assert(std::abs(ret1 - ret3) < 1e-14);
+
+  return ret1;
+}

src/compute_state.hh

 #ifndef COMPUTE_STATE_HH
 #define COMPUTE_STATE_HH
 
-double compute_state_update_bisection(double h, double unorm, double L,
-                                      double old_state);
-double compute_state_update_lambert(double h, double unorm, double L,
-                                    double old_state);
-double compute_state_update_lambert_gsl(double h, double unorm, double L,
-                                        double old_state);
+double compute_state_update(double h, double unorm, double L, double old_state);
 
 #endif

src/one-body-sample.cc

@@ -349,26 +349,7 @@ int main(int argc, char *argv[]) {
               double const L = parset.get<double>("boundary.friction.ruina.L");
               double const unorm = u4_diff[i].two_norm();
 
-              double ret1 =
-                  compute_state_update_bisection(h, unorm, L, s4_old[i][0]);
-              assert(std::abs(1.0 / h * ret1 - (s4_old[i] - unorm / L) / h -
-                              std::exp(-ret1)) < 1e-10);
-
-              double ret2 =
-                  compute_state_update_lambert(h, unorm, L, s4_old[i][0]);
-              assert(std::abs(1.0 / h * ret2 - (s4_old[i] - unorm / L) / h -
-                              std::exp(-ret2)) < 1e-10);
-
-              double ret3 =
-                  compute_state_update_lambert_gsl(h, unorm, L, s4_old[i][0]);
-              assert(std::abs(1.0 / h * ret3 - (s4_old[i] - unorm / L) / h -
-                              std::exp(-ret3)) < 1e-10);
-
-              assert(std::abs(ret1 - ret2) < 1e-14);
-              assert(std::abs(ret1 - ret3) < 1e-14);
-              assert(std::abs(ret2 - ret3) < 1e-14);
-
-              (*s4_new)[i][0] = ret1;
+              (*s4_new)[i][0] = compute_state_update(h, unorm, L, s4_old[i][0]);
             }
           }
         }
@@ -421,26 +402,7 @@ int main(int argc, char *argv[]) {
               double const L = parset.get<double>("boundary.friction.ruina.L");
               double const unorm = u5_diff[i].two_norm();
 
-              double ret1 =
-                  compute_state_update_bisection(h, unorm, L, s5_old[i][0]);
-              assert(std::abs(1.0 / h * ret1 - (s5_old[i] - unorm / L) / h -
-                              std::exp(-ret1)) < 1e-12);
-
-              double ret2 =
-                  compute_state_update_lambert(h, unorm, L, s5_old[i][0]);
-              assert(std::abs(1.0 / h * ret2 - (s5_old[i] - unorm / L) / h -
-                              std::exp(-ret2)) < 1e-12);
-
-              double ret3 =
-                  compute_state_update_lambert_gsl(h, unorm, L, s4_old[i][0]);
-              assert(std::abs(1.0 / h * ret3 - (s4_old[i] - unorm / L) / h -
-                              std::exp(-ret3)) < 1e-10);
-
-              assert(std::abs(ret1 - ret2) < 1e-14);
-              assert(std::abs(ret1 - ret3) < 1e-14);
-              assert(std::abs(ret2 - ret3) < 1e-14);
-
-              (*s5_new)[i][0] = ret1;
+              (*s5_new)[i][0] = compute_state_update(h, unorm, L, s5_old[i][0]);
             }
           }
         }