Improve names of some variables, added support for T = 0 bounce solve…

…r and added unit tests for it.

include/BSMPT/bounce_solution/action_calculation.h

@@ -73,7 +73,7 @@ class BounceActionInt
    * @brief l0 - Initial_lmin for solutions starting very near the true vacuum
    *
    */
-  double l0Initiallmin;
+  double l0_minus_lmin;
 
   /**
    * @brief True if path deformation reached the desired results without solving
@@ -276,6 +276,12 @@ class BounceActionInt
    */
   tk::spline RasterizeddVdl; // RasterizeddVdl
 
+  /**
+   * @brief Default constructor (unit tests)
+   *
+   */
+  BounceActionInt();
+
   /**
    * @brief Construct a new Bounce Action Int object
    *

src/bounce_solution/action_calculation.cpp

@@ -12,6 +12,8 @@
 namespace BSMPT
 {
 
+BounceActionInt::BounceActionInt(){};
+
 BounceActionInt::BounceActionInt(
     std::vector<std::vector<double>> InitPath_In,
     std::vector<double> TrueVacuum_In,
@@ -418,7 +420,7 @@ void BounceActionInt::RK5_step(std::vector<double> y,
   return;
 }
 
-double BounceActionInt::BesselI(double Alpha, double x, int terms)
+double BounceActionInt::BesselI(double alpha, double x, int terms)
 {
   // This implementation seems to converge quite quicly
   // https://en.wikipedia.org/wiki/Bessel_function#:~:text=Modified%20Bessel%20functions%3A%20I%CE%B1%2C%20K%CE%B1%5B,first%20and%20second%20kind%20and%20are%20defined%20as%5B19%5D
@@ -428,8 +430,8 @@ double BounceActionInt::BesselI(double Alpha, double x, int terms)
   while ((m < terms) && (abs((r - r0) / r) > 1e-15))
   {
     r0 = r; // Save step
-    r += 1 / (tgamma(m + Alpha + 1) * tgamma(m + 1)) *
-         std::pow(x / 2.0, 2.0 * m + Alpha); // One more term
+    r += 1 / (tgamma(m + alpha + 1) * tgamma(m + 1)) *
+         std::pow(x / 2.0, 2.0 * m + alpha); // One more term
     m++; // Update later to not mess up summation
   }
   return r;
@@ -479,19 +481,37 @@ std::vector<double> BounceActionInt::ExactSolutionFromMinimum(double l0,
   {
     LinearSolution = [=](const double rho_in)
     {
-      return l - (Initial_lmin + l0Initiallmin *
+      return l - (Initial_lmin + l0_minus_lmin *
                                      sinh(sqrt(TrueVacuumHessian) * rho_in) /
                                      (sqrt(TrueVacuumHessian) * rho_in));
     };
     LinearSolutionDerivative = [=](const double rho_in)
     {
-      return ((l0Initiallmin * cosh(rho_in * std::sqrt(TrueVacuumHessian))) /
+      return ((l0_minus_lmin * cosh(rho_in * std::sqrt(TrueVacuumHessian))) /
               rho_in) -
-             (l0Initiallmin * sinh(rho_in * std::sqrt(TrueVacuumHessian))) /
+             (l0_minus_lmin * sinh(rho_in * std::sqrt(TrueVacuumHessian))) /
                  (pow(rho_in, 2) * std::sqrt(TrueVacuumHessian));
     };
   }
 
+  // T = 0
+  if (Alpha == 3)
+  {
+    LinearSolution = [=](const double rho_in)
+    {
+      return l -
+             (Initial_lmin + 2 * l0_minus_lmin *
+                                 BesselI(1, sqrt(TrueVacuumHessian) * rho_in) /
+                                 (sqrt(TrueVacuumHessian) * rho_in));
+    };
+    LinearSolutionDerivative = [=](const double rho_in)
+    {
+      return (
+          (2 * l0_minus_lmin * BesselI(2, rho_in * sqrt(TrueVacuumHessian))) /
+          rho_in);
+    };
+  }
+
   // Check if lower limit is viable
   assert(LinearSolution(rho_down) > 0);
 
@@ -1022,7 +1042,7 @@ void BounceActionInt::Solve1DBounce(
     if (mode == 0)
     {
       l0            = (lmax + lmin) / 2.0; // Perform binary search
-      l0Initiallmin = l0 - Initial_lmin;
+      l0_minus_lmin = l0 - Initial_lmin;
       IntegrateBounce(
           l0, conv, rho, l, dl_drho, d2l_drho2, 100000, error, error * 0.0015);
       if (StateOfBounceActionInt != ActionStatus::NotCalculated) return;
@@ -1083,8 +1103,8 @@ void BounceActionInt::Solve1DBounce(
     if (mode == 1)
     {
       mu_middle     = (mu_min + mu_max) / 2;
-      l0Initiallmin = exp(mu_middle);
-      l0 = Initial_lmin + l0Initiallmin; // Perform binary search in log space
+      l0_minus_lmin = exp(mu_middle);
+      l0 = Initial_lmin + l0_minus_lmin; // Perform binary search in log space
 
       IntegrateBounce(
           l0, conv, rho, l, dl_drho, d2l_drho2, 100000, error, error * 0.0015);
@@ -1632,34 +1652,67 @@ double BounceActionInt::CalculateKineticTermAction(std::vector<double> &rho,
 {
   double integral  = 0;
   double int_delta = rho[rho.size() - 1] / 2000;
-  for (double r = 0.0; r <= rho[rho.size() - 1]; r += int_delta)
+  if (Alpha == 2)
   {
-    // Simpson Integration (1 + 4 + 1)/ 6 * step
-    integral += r * r * (0.5 * std::pow(dl_drho_spl(r), 2));
-    integral +=
-        4 * r * r * (0.5 * std::pow(dl_drho_spl(r + int_delta / 2.0), 2));
-    integral += r * r * (0.5 * std::pow(dl_drho_spl(r + int_delta), 2));
+    for (double r = 0.0; r <= rho[rho.size() - 1]; r += int_delta)
+    {
+      // Simpson Integration (1 + 4 + 1)/ 6 * step
+      integral += r * r * (0.5 * std::pow(dl_drho_spl(r), 2));
+      integral +=
+          4 * r * r * (0.5 * std::pow(dl_drho_spl(r + int_delta / 2.0), 2));
+      integral += r * r * (0.5 * std::pow(dl_drho_spl(r + int_delta), 2));
+    }
+    integral = integral * 4 * M_PI * int_delta /
+               6.0; // Angular integration and Simpson step
+    return integral;
+  }
+  else if (Alpha == 3)
+  {
+    for (double r = 0.0; r <= rho[rho.size() - 1]; r += int_delta)
+    {
+      // Simpson Integration (1 + 4 + 1)/ 6 * step
+      integral += r * r * r * (0.5 * std::pow(dl_drho_spl(r), 2));
+      integral +=
+          4 * r * r * r * (0.5 * std::pow(dl_drho_spl(r + int_delta / 2.0), 2));
+      integral += r * r * r * (0.5 * std::pow(dl_drho_spl(r + int_delta), 2));
+    }
+    integral = integral * 2 * M_PI * M_PI * int_delta /
+               6.0; // Angular integration and Simpson step
+    return integral;
   }
-  integral = integral * 4 * M_PI * int_delta /
-             6.0; // Angular integration and Simpson step
-  return integral;
 }
 
 double BounceActionInt::CalculatePotentialTermAction(std::vector<double> &rho,
                                                      tk::spline &l_rho_spl)
 {
   double integral  = 0;
   double int_delta = rho[rho.size() - 1] / 2000;
-  for (double r = 0.0; r <= rho[rho.size() - 1]; r += int_delta)
+  if (Alpha == 2)
   {
-    // Simpson Integration (1 + 4 + 1)/ 6 * step
-    integral += r * r * (V(Spline(l_rho_spl(r))));
-    integral += 4 * r * r * (V(Spline(l_rho_spl(r + int_delta / 2.0))));
-    integral += r * r * (V(Spline(l_rho_spl(r + int_delta))));
+    for (double r = 0.0; r <= rho[rho.size() - 1]; r += int_delta)
+    {
+      // Simpson Integration (1 + 4 + 1)/ 6 * step
+      integral += r * r * (V(Spline(l_rho_spl(r))));
+      integral += 4 * r * r * (V(Spline(l_rho_spl(r + int_delta / 2.0))));
+      integral += r * r * (V(Spline(l_rho_spl(r + int_delta))));
+    }
+    integral = integral * 4 * M_PI * int_delta /
+               6.0; // Angular integration and Simpson step
+    return integral;
+  }
+  else if (Alpha == 3)
+  {
+    for (double r = 0.0; r <= rho[rho.size() - 1]; r += int_delta)
+    {
+      // Simpson Integration (1 + 4 + 1)/ 6 * step
+      integral += r * r * r * (V(Spline(l_rho_spl(r))));
+      integral += 4 * r * r * r * (V(Spline(l_rho_spl(r + int_delta / 2.0))));
+      integral += r * r * r * (V(Spline(l_rho_spl(r + int_delta))));
+    }
+    integral = integral * 2 * M_PI * M_PI * int_delta /
+               6.0; // Angular integration and Simpson step
+    return integral;
   }
-  integral = integral * 4 * M_PI * int_delta /
-             6.0; // Angular integration and Simpson step
-  return integral;
 }
 
 void BounceActionInt::CalculateAction(

tests/unittests/Test-gw.cpp

@@ -19,6 +19,18 @@ using Approx = Catch::Approx;
 #include <BSMPT/utility/Logger.h> // for Logger Class
 #include <fstream>
 
+TEST_CASE("Test I_\alpha", "[gw]")
+{
+  // Tests bounce solver with analytical derivative
+  using namespace BSMPT;
+
+  BounceActionInt BACalc;
+
+  REQUIRE(BACalc.BesselI(3, 1) == Approx(0.0221684249).epsilon(5e-2));
+  REQUIRE(BACalc.BesselI(1, 3) == Approx(3.953370217).epsilon(5e-2));
+  REQUIRE(BACalc.BesselI(1, 1.5) == Approx(0.9816664).epsilon(5e-2));
+}
+
 TEST_CASE("Solve bounce equation with analytical derivative", "[gw]")
 {
   // Tests bounce solver with analytical derivative
@@ -61,6 +73,51 @@ TEST_CASE("Solve bounce equation with analytical derivative", "[gw]")
   REQUIRE(bc.Action == Approx(4.5011952256).epsilon(5e-2));
 }
 
+TEST_CASE(
+    "Solve bounce equation with analytical derivative and Alpha = 3 (T = 0)",
+    "[gw]")
+{
+  // Tests bounce solver with analytical derivative
+  using namespace BSMPT;
+  std::function<double(std::vector<double>)> V = [&](std::vector<double> x)
+  {
+    double c  = 5;
+    double fx = 0;
+    double fy = 80;
+
+    double r1 = x[0] * x[0] + c * x[1] * x[1];
+    double r2 = c * pow(x[0] - 1, 2) + pow(x[1] - 1, 2);
+    double r3 = fx * (0.25 * pow(x[0], 4) - pow(x[0], 3) / 3.);
+    r3 += fy * (0.25 * pow(x[1], 4) - pow(x[1], 3) / 3.);
+
+    return (r1 * r2 + r3);
+  };
+
+  std::function<std::vector<double>(std::vector<double>)> dV =
+      [&](std::vector<double> l0)
+  {
+    int dim = 2;
+    std::vector<double> result(dim);
+    result = {2 * l0[0] * (5 * pow(-1 + l0[0], 2) + pow(-1 + l0[1], 2)) +
+                  10 * (-1 + l0[0]) * (pow(l0[0], 2) + 5 * pow(l0[1], 2)),
+              10 * (5 * pow(-1 + l0[0], 2) + pow(-1 + l0[1], 2)) * l0[1] +
+                  2 * (-1 + l0[1]) * (pow(l0[0], 2) + 5 * pow(l0[1], 2)) +
+                  80 * (-1. * pow(l0[1], 2) + 1. * pow(l0[1], 3))};
+    return result;
+  };
+
+  std::vector<double> FalseVacuum = {0, 0};
+  std::vector<double> TrueVacuum  = {1, 1};
+
+  std::vector<std::vector<double>> path = {TrueVacuum, FalseVacuum};
+
+  BounceActionInt bc(path, TrueVacuum, FalseVacuum, V, 0, 6);
+  bc.Alpha = 3;
+  bc.CalculateAction();
+
+  REQUIRE(bc.Action == Approx(13.3129767888).epsilon(5e-2));
+}
+
 TEST_CASE("Solve bounce equation with numerical derivative", "[gw]")
 {
   // Tests bounce solver with numerical derivative
@@ -226,6 +283,7 @@ TEST_CASE("Solve bounce equation with numerical derivative thin walled", "[gw]")
 {
   // Tests bounce solver with numerical derivative
   using namespace BSMPT;
+  SetLogger({"--logginglevel::complete"});
 
   std::function<double(std::vector<double>)> V = [&](std::vector<double> x)
   {