23#include <gsl/gsl_errno.h>
24#include <gsl/gsl_integration.h>
25#include <gsl/gsl_sf_bessel.h>
26#include <gsl/gsl_sf_expint.h>
30#define M_EULER_MASCHERONI 0.5772156649015328606065120
64static const double tol = 1e-7;
65static double zero = 1e-100;
67static double f(
double w,
void *p) {
69 const double kx = (params->
kx);
70 double kr = (params->
kr);
71 double x0 = (params->
x0);
72 double r0 = (params->
r0);
90static double f_kx_zero(
double w,
void *p) {
92 double kr = (params->
kr);
93 double x0 = (params->
x0);
99 result = 2.0 * M_PI * w * gsl_sf_bessel_J0(
kr * w) *
117 double result, error;
118 double xinf = 100.0 /
r0;
122 const size_t wrk_size = 5000;
127 const double epsilon_abs = 1e-3;
128 const double epsilon_rel = 1e-3;
131 gsl_integration_workspace *ws = gsl_integration_workspace_alloc(wrk_size);
136 assert(which_int == PFC_F_0 || which_int == PFC_F_KR ||
137 which_int == PFC_F_KX);
151 gsl_integration_qag(&F, -xinf, xinf, epsilon_abs, epsilon_rel, wrk_size, 3,
152 ws, &result, &error);
156 gsl_integration_qag(&F, 0.0,
r0, epsilon_abs, epsilon_rel, wrk_size, 3, ws,
161 gsl_integration_qag(&F, 0.0,
x0, epsilon_abs, epsilon_rel, wrk_size, 3, ws,
167 gsl_integration_workspace_free(ws);
178 4.0 *
kr *
r0 * gsl_sf_bessel_Jn(1,
kr *
r0) *
180 4.0 * gsl_sf_bessel_Jn(0,
kr *
r0) *
189 result = -2.0 * M_PI *
194 return result / (4.0 * M_PI);
198double FC_FUNC_(c_poisson_cutoff_3d_1d_finite,
199 C_POISSON_CUTOFF_3D_1D_FINITE)(
double *gx,
double *gperp,
200 double *xsize,
double *rsize) {
209static double bessel_J0(
double w,
void *p) {
return gsl_sf_bessel_J0(w); }
212double FC_FUNC_(c_poisson_cutoff_2d_0d, C_POISSON_CUTOFF_2D_0D)(
double *x,
214 double result, error;
215 const size_t wrk_size = 500;
216 const double epsilon_abs = 1e-3;
217 const double epsilon_rel = 1e-3;
219 gsl_integration_workspace *ws = gsl_integration_workspace_alloc(wrk_size);
223 gsl_integration_qag(&F, *x, *y, epsilon_abs, epsilon_rel, wrk_size, 3, ws,
225 gsl_integration_workspace_free(ws);
237double FC_FUNC(intcoslog, INTCOSLOG)(
double *mu,
double *a,
double *b) {
238 if (
fabs(*a) > 0.0) {
239 return (1.0 / (*a)) *
240 (
log((*b) * (*mu)) *
sin((*a) * (*mu)) - gsl_sf_Si((*a) * (*mu)));
242 return (*mu) * (
log((*mu) * (*b)) - 1.0);
260 double gx = (params->
gx);
261 double gy = (params->
gy);
265 k0 = -(
log(k0arg / 2) + M_EULER_MASCHERONI);
266 }
else if (k0arg < 50.0) {
267 k0 = gsl_sf_bessel_K0(k0arg);
269 k0 =
sqrt(2.0 * M_PI / k0arg) *
exp(-k0arg);
271 return 4.0 *
cos(w *
gy) * k0;
275double FC_FUNC_(c_poisson_cutoff_2d_1d,
276 C_POISSON_CUTOFF_2D_1D)(
double *
gy,
double *
gx,
double *
rc) {
277 double result, error, res, b;
279 const size_t wrk_size = 5000;
280 const double epsilon_abs = 1e-3;
281 const double epsilon_rel = 1e-3;
284 gsl_integration_workspace *ws = gsl_integration_workspace_alloc(wrk_size);
297 res = -4.0 * FC_FUNC(intcoslog, INTCOSLOG)(&mu,
gy, &b);
300 res = res - (4.0 * M_EULER_MASCHERONI / (*gy)) *
sin((*
gy) * mu);
302 res = res - 4.0 * M_EULER_MASCHERONI * mu;
305 gsl_integration_qag(&F, mu, (*
rc), epsilon_abs, epsilon_rel, wrk_size, 3, ws,
309 gsl_integration_workspace_free(ws);
326 double g = (params->
g);
327 double a = (params->
a);
328 return 2.0 * (
cos(w) /
sqrt(w * w +
a *
a *
g *
g));
332double FC_FUNC_(c_poisson_cutoff_1d_0d,
333 C_POISSON_CUTOFF_1D_0D)(
double *
g,
double *
a,
double *rc) {
334 double result, error;
336 const size_t wrk_size = 5000;
337 const double epsilon_abs = 1e-3;
338 const double epsilon_rel = 1e-3;
340 gsl_integration_workspace *ws;
344 return 2.0 *
asinh((*rc) / (*
a));
349 if ((*
g) * (*rc) > 100.1 * M_PI) {
350 return 2.0 * gsl_sf_bessel_K0((*
a) * (*
g));
353 gsl_set_error_handler_off();
355 ws = gsl_integration_workspace_alloc(wrk_size);
363 status = gsl_integration_qag(&F, 0.0, (*rc) * (*
g), epsilon_abs, epsilon_rel,
364 wrk_size, 3, ws, &result, &error);
366 gsl_integration_workspace_free(ws);
double pow(double __x, double __y) __attribute__((__nothrow__
double log(double __x) __attribute__((__nothrow__
double exp(double __x) __attribute__((__nothrow__
double sin(double __x) __attribute__((__nothrow__
double fabs(double __x) __attribute__((__nothrow__
double sqrt(double __x) __attribute__((__nothrow__
double asinh(double __x) __attribute__((__nothrow__
double cos(double __x) __attribute__((__nothrow__
const int *restrict index
static double cutoff_1d_0d(double w, void *p)
static double f_kx_zero(double w, void *p)
static double int_aux(int index, double kx, double kr, double x0, double r0, int which_int)
static double cutoff_2d_1d(double w, void *p)
double poisson_finite_cylinder(double kx, double kr, double x0, double r0)
static double f_kr_zero(double w, void *p)
static double bessel_J0(double w, void *p)
static double f(double w, void *p)