FGRadPot.hpp

#pragma once
#include <array>
#include <cassert>
#include <cmath>
#include <functional>
#include <vector>
 
 
namespace FGRP {
 
// Fine structure constant:
constexpr double alpha = 1.0 / 137.035999084;
// electron reduced compton wavelength (atomic units)
constexpr double lam_c = alpha;
 
//==============================================================================
double V_Uehling(double Z, double r, double rN = 0.0);
 
double H_Magnetic(double Z, double r, double rN = 0.0);
 
double V_SEh(double Z, double r, double rN = 0.0);
 
double V_SEl(double Z, double r, double rN = 0.0);
 
double V_WK(double Z, double r, double);
 
//==============================================================================
namespace Fit {
// static const auto a01 = std::vector{1.071, 0.0, -1.976, -2.128, 0.169};
static const auto a01 = std::vector{1.071, 0.0, -1.976, -2.128, 0.169};
static const auto b01 = std::vector{0.074, 0.35};
static const auto b2 = std::vector{0.056, 0.050, 0.195};
 
double A(double Z, int l = 0);
double B(double Z, int l = 0);
 
} // namespace Fit
 
//==============================================================================
double t_integral(double (*f)(double, void *), std::vector<double> params,
                  double eps = 1.0e-6);
 
enum class IntType { Linear, Log };
 
template <IntType IT = IntType::Linear>
double r_integral(std::function<double(double)> f, double a, double b,
                  long unsigned n_pts = 1000);
 
//==============================================================================
// Helper functions:
 
namespace Uehling {
double G_Ueh(double xn, double x);
// Form of function required by GSL integrations; p is {r, Rn}
double J_Ueh_gsl(double t, void *p);
} // namespace Uehling
 
namespace Magnetic {
double G_mag(double xn, double x);
// Form of function required by GSL integrations; p is {r, Rn}
double J_mag_gsl(double t, void *p);
} // namespace Magnetic
 
namespace SE {
double xi(double x);
double F_SEl(double Z, double r, double rN);
double I1(double Z, double t);
double I2(double Z, double r, double rn, double t);
// Form of function required by GSL integrations; p is {r, Rn, Z}
double J_SE_gsl(double t, void *p);
} // namespace SE
 
//==============================================================================
// Implementation:
template <IntType IT>
double r_integral(std::function<double(double)> f, double a, double b,
                  long unsigned n_pts) {
 
  static constexpr std::array w{475.0 / 1440, 1902.0 / 1440, 1104.0 / 1440,
                                1586.0 / 1440, 1413.0 / 1440};
 
  const auto dt = (IT == IntType::Linear) ? (b - a) / double(n_pts - 1) :
                                            std::log(b / a) / double(n_pts - 1);
 
  const auto x = [=](auto i) {
    if constexpr (IT == IntType::Linear) {
      return a + double(i) * dt;
    } else {
      assert(a != 0.0);
      return a * std::exp(double(i) * dt);
    }
  };
  const auto dxdt = [=](auto i) {
    if constexpr (IT == IntType::Linear) {
      (void)i; // suppress unused variable warning on older gcc compilers
      return 1.0;
    } else {
      return x(i); // dxdt = x(t) for log
    }
  };
 
  double Rint = 0.0;
  for (long unsigned i = 0; i < w.size(); i++) {
    Rint += w[i] * (f(x(i)) + f(x(n_pts - 1 - i))) * dxdt(i);
  }
 
  for (auto i = w.size(); i < n_pts - w.size(); i++) {
    Rint += f(x(i)) * dxdt(i);
  }
 
  return Rint * dt;
}
 
//==============================================================================
//==============================================================================
double Q_MLVP(double r, double rN = 0.0);
namespace Helper {
// Magnetic-loop vacuum polarisation, Q = Int(Qt, {t,1,infty}).
double Qt_MLVP(double t, void *p);
} // namespace Helper
 
// struct to pass the parameters to the GSL function
struct MLVP_params {
  // simple struct that only stores current point on the radial grid and the
  // nuclear radius
  double r, rN;
};
 
} // namespace FGRP