euclidean_action.cpp

Go to the documentation of this file.

/*
 * This file is part of BubbleProfiler.
 *
 * BubbleProfiler is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * BubbleProfiler is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with BubbleProfiler.  If not, see <http://www.gnu.org/licenses/>.
 */

#include "euclidean_action.hpp"
#include "error.hpp"
#include "field_profiles.hpp"
#include "numeric.hpp"
#include "potential.hpp"

#include <cmath>

namespace BubbleProfiler {

double calculate_kinetic_action(
   const Potential& /* potential */, const Field_profiles& profiles,
   std::size_t max_intervals, double rel_tol, double abs_tol,
   Integration_rule rule)
{
   const int n_fields = profiles.get_number_of_fields();
   const int d = profiles.get_number_of_dimensions();
   const double domain_start = profiles.get_domain_start();
   const double domain_end = profiles.get_domain_end();

   // finite domain correction
   double correction = 0;
   for (int i = 0; i < n_fields; ++i) {
      const double deriv = profiles.derivative_at(i, 1, domain_end);
      correction += 0.5 * deriv * deriv;
   }

   const auto integrand
      = [&profiles, n_fields, d, correction](double rho)
   {
      double result = 0.;
      for (int i = 0; i < n_fields; ++i) {
         const double deriv = profiles.derivative_at(i, 1, rho);
         result += 0.5 * deriv * deriv;
      }

      result -= correction;
      result *= std::pow(rho, d - 1);

      return result;
   };

   const auto result = integrate_gsl_qag(integrand, domain_start, domain_end,
                                         abs_tol, rel_tol, max_intervals, rule);

   const auto error = std::get<0>(result);
   if (error) {
      throw Numerical_error("calculation of action failed with error code: "
                            + std::to_string(error));
   }

   return 2. * area_n_sphere(d - 1) * std::get<1>(result) / d;
}

double calculate_potential_action(
   const Potential& potential, const Field_profiles& profiles,
   std::size_t max_intervals, double rel_tol, double abs_tol,
   Integration_rule rule)
{
   const int d = profiles.get_number_of_dimensions();
   const double domain_start = profiles.get_domain_start();
   const double domain_end = profiles.get_domain_end();

   // finite domain correction
   double correction = potential(profiles.evaluate_at(domain_end));

   const auto integrand
      = [&potential, &profiles, d, correction](double rho)
   {
      double result = potential(profiles.evaluate_at(rho));
      result -= correction;
      result *= std::pow(rho, d - 1);

      return result;
   };

   const auto result = integrate_gsl_qag(integrand, domain_start, domain_end,
                                         abs_tol, rel_tol, max_intervals, rule);

   const auto error = std::get<0>(result);
   if (error) {
      throw Numerical_error("calculation of action failed with error code: "
                            + std::to_string(error));
   }

   return 2. * area_n_sphere(d - 1) * std::get<1>(result) / (2. - d);
}

double calculate_full_action(
   const Potential& potential, const Field_profiles& profiles,
   std::size_t max_intervals, double rel_tol,
   double abs_tol, Integration_rule rule)
{
   const int n_fields = profiles.get_number_of_fields();
   const int d = profiles.get_number_of_dimensions();
   const double domain_start = profiles.get_domain_start();
   const double domain_end = profiles.get_domain_end();

   // finite domain correction
   double correction = 0.;
   for (int i = 0; i < n_fields; ++i) {
      const double deriv = profiles.derivative_at(i, 1, domain_end);
      correction += 0.5 * deriv * deriv;
   }
   correction += potential(profiles.evaluate_at(domain_end));

   const auto integrand
      = [&potential, &profiles, n_fields, d, correction](double rho)
      {
         double result = 0.;
         for (int i = 0; i < n_fields; ++i) {
            const double deriv = profiles.derivative_at(i, 1, rho);
            result += 0.5 * deriv * deriv;
         }

         result += potential(profiles.evaluate_at(rho));
         result -= correction;
         result *= pow(rho, d - 1);

         return result;
      };

   const auto result = integrate_gsl_qag(integrand, domain_start, domain_end,
                                         abs_tol, rel_tol, max_intervals, rule);

   const auto error = std::get<0>(result);
   if (error) {
      throw Numerical_error("calculation of action failed with error code: "
                            + std::to_string(error));
   }

   return area_n_sphere(d - 1) * std::get<1>(result);
}

double calculate_action(
   const Potential& potential, const Field_profiles& profiles,
   std::size_t max_intervals, double rel_tol,
   double abs_tol, Integration_rule rule,
   bool use_kinetic)
{
   return use_kinetic ?
      calculate_kinetic_action(potential, profiles, max_intervals,
                               rel_tol, abs_tol, rule) :
      calculate_potential_action(potential, profiles, max_intervals,
                                 rel_tol, abs_tol, rule);
}

} // namespace BubbleProfiler