53#include <gsl/gsl_integration.h>
54#include <gsl/gsl_sf_bessel.h>
55#include <gsl/gsl_sf_erf.h>
56#include <gsl/gsl_sf_gamma.h>
64#include <lal/LALConstants.h>
65#include <lal/XLALError.h>
66#include <lal/LALMarcumQ.h>
99static double Q(
double mu,
double x)
101 return gsl_sf_gamma_inc_Q(
mu,
x);
119 return gsl_sf_bessel_Inu_scaled(
mu, xi) / gsl_sf_bessel_Inu_scaled(
mu - 1., xi);
131 return gsl_sf_lngamma(
mu + n) - gsl_sf_lngamma(
mu - n) - n * log(2.) - gsl_sf_lngamma(n + 1);
143 if(fabs(
y -
x - 1.) < 1e-3) {
158 +((72. *
x + 42.) *
x + 7.) / 36.,
159 -(((2700. *
x + 2142.) *
x + 657.) *
x + 73.) / 540.,
160 +((((181440. *
x + 177552.) *
x + 76356.) *
x + 15972.) *
x + 1331.) / 12960.
162 double z = (
y -
x - 1.) / pow(2. *
x + 1., 2.);
163 double z_to_k_p_1 = z;
164 double sum = z_to_k_p_1;
166 for(k = 1; k < 5; k++) {
168 sum += c[k] * z_to_k_p_1;
170 return pow(2. *
x + 1., 3.) * sum * sum / 2.;
172 double root_1_plus_4xy = sqrt(1. + 4. *
x *
y);
173 return x +
y - root_1_plus_4xy + log((1. + root_1_plus_4xy) / (2. *
y));
200 return theta < 1e-4 ? theta * theta / 6. + 1. : theta / sin(theta);
218 return theta < 1e-4 ? (theta * theta / 45. + 1. / 3.) * theta * theta : 1. - theta / tan(theta);
222static double rho(
double theta_over_sin_theta_,
double xi)
224 return sqrt(theta_over_sin_theta_ * theta_over_sin_theta_ + xi * xi);
228static double r(
double theta,
double y,
double xi)
232 xi /= theta_over_sin_theta_;
233 return (1. + sqrt(1. + xi * xi)) * theta_over_sin_theta_ / (2. *
y);
241 xi /= theta_over_sin_theta_;
258static double f(
double theta,
double y,
double xi)
260 double r_ =
r(theta,
y, xi);
261 double r_minus_cos_theta = r_ - cos(theta);
262 double sin_theta = sin(theta);
264 return (
r_primed_sin_theta(theta,
y, xi) - r_minus_cos_theta * r_) / (r_minus_cos_theta * r_minus_cos_theta + sin_theta * sin_theta);
273static double psi(
double theta,
double xi)
276 double rho_ =
rho(theta_over_sin_theta_, xi);
277 double root_xi2_plus_1 = sqrt(1. + xi * xi);
279 return cos(theta) * rho_ - root_xi2_plus_1 - log((theta_over_sin_theta_ + rho_) / (1. + root_xi2_plus_1));
288static void f1_f2(
double x,
double M,
double *f1,
double *f2)
291 double b = sqrt(4. *
x + 2. *
M);
314 double x_to_n_over_n_factorial = 1.;
318 double term = x_to_n_over_n_factorial *
Q(
M + n,
y);
320 if(term <= 1e-17 * sum)
323 x_to_n_over_n_factorial *=
x / n;
326 return exp(-
x) * sum;
338 double root_y_minus_root_x = sqrt(
y) - sqrt(
x);
341 double sigma = root_y_minus_root_x * root_y_minus_root_x / xi;
342 double rho_ = sqrt(
y /
x);
344 double rho_to_M_over_root_8_pi = pow(
y /
x,
M / 2.) / sqrt(8. *
LAL_PI);
345 double e_to_neg_sigma_xi_over_xi_to_n_minus_half = exp(-root_y_minus_root_x * root_y_minus_root_x) * sqrt(xi);
374 double Phi = fabs(root_y_minus_root_x) < 1e-5 ? sigma == 0. ? 0. : sqrt(
LAL_PI) / sqrt(sigma) - 2. * sqrt(xi) : sqrt(
LAL_PI / sigma) * gsl_sf_erfc(fabs(root_y_minus_root_x));
413 double Psi = rho_ == 1. ? .5 : copysign(pow(rho_,
M - .5) * gsl_sf_erfc(fabs(root_y_minus_root_x)) / 2., rho_ - 1.);
415 double sum =
x >
y ? 1. : 0.;
421 if(fabs(Psi) <= 1e-17 * fabs(sum))
427 rho_to_M_over_root_8_pi = -rho_to_M_over_root_8_pi;
428 e_to_neg_sigma_xi_over_xi_to_n_minus_half /= xi;
431 Phi = (e_to_neg_sigma_xi_over_xi_to_n_minus_half - sigma * Phi) / (n - .5);
434 double lnA_n_M_minus_1 =
lnA(n,
M - 1.);
435 Psi = rho_to_M_over_root_8_pi * exp(lnA_n_M_minus_1) * (1. - exp(
lnA(n,
M) - lnA_n_M_minus_1) / rho_) * Phi;
443 if(-1e-200 < sum && sum < 0.)
458 double root_y_over_x = sqrt(
y /
x);
461 double mu = sqrt(2. * xi) - 1.;
469 double cmu = root_y_over_x *
c_mu(
mu, xi);
472 double Qmu_saved = Qmu;
473 Qmu = (1. + cmu) * Qmu - cmu * Qmu_minus_1;
474 Qmu_minus_1 = Qmu_saved;
478 if(fabs(
mu -
M) /
M > 1e-15)
493 double zeta_ =
zeta(
x,
y);
499 Psi[0] = sqrt(
LAL_PI / (2. *
M)) * gsl_sf_erfc(-zeta_ * sqrt(
M / 2.));
500 Psi[1] = e_to_neg_half_M_zeta_2 /
M;
509 for(j = 0; j <= k; j++)
510 Bk += Psi[j] / pow(
M, k - j);
513 if(Bk <= 1e-17 * sum)
525 Psi[k] = (k - 1) /
M * Psi[k - 2] + pow(-zeta_, k - 1) /
M * e_to_neg_half_M_zeta_2;
529 return gsl_sf_erfc(-zeta_ * sqrt(
M / 2.)) / 2. - sqrt(
M /
LAL_TWOPI) * sum;
565 return exp(params.
M *
psi(theta, params.
xi)) *
f(theta, params.
y, params.
xi);
576 gsl_function integrand_ = {
580 gsl_integration_workspace *workspace = gsl_integration_workspace_alloc(20);
581 double integral = 0.;
586 gsl_integration_qag(&integrand_, 0,
LAL_PI - 1./512., abserr, 1e-12, 20, GSL_INTEG_GAUSS41, workspace, &integral, &abserr);
587 gsl_integration_workspace_free(workspace);
654 xi = 2. * sqrt(
x *
y);
667 }
else if(
xi > 30. &&
M *
M < 2. *
xi) {
673 }
else if(f1 <
y &&
y < f2 &&
M < 135.) {
679 }
else if(f1 <
y && y < f2 && M >= 135.) {
703 if(1. < Q_ && Q_ < 1. + 1e-12)
706 if(Q_ < 0. || Q_ > 1.)
static double theta_over_sin_theta(double theta)
double XLALMarcumQmodified(double M, double x, double y)
The modified form of the Marcum Q function.
static void f1_f2(double x, double M, double *f1, double *f2)
static double MarcumQ_quadrature(double M, double x, double y, double xi)
static double f(double theta, double y, double xi)
static double r(double theta, double y, double xi)
static double integrand(double theta, void *integrand_params)
static double MarcumQ_recurrence(double M, double x, double y, double xi)
static double r_primed_sin_theta(double theta, double y, double xi)
static double MarcumQ_large_M(double M, double x, double y)
static double rho(double theta_over_sin_theta_, double xi)
static double theta_over_sin_theta_primed_sin_theta(double theta)
static double MarcumQ_small_x(double M, double x, double y)
static double c_mu(double mu, double xi)
static double half_zeta_2(double x, double y)
static double psi(double theta, double xi)
static double MarcumQ_large_xy(double M, double x, double y, double xi)
static double Q(double mu, double x)
double XLALMarcumQ(double M, double a, double b)
The function defined by J. Marcum,.
static double zeta(double x, double y)
static double lnA(int n, double mu)
#define LAL_PI
Archimedes's constant, pi.
#define LAL_TWOPI
2*pi is circumference of a circle divided by its radius
#define XLAL_ERROR_REAL8(...)
Macro to invoke a failure from a XLAL routine returning a REAL8.
@ XLAL_ERANGE
Output range error.
@ XLAL_EMAXITER
Exceeded maximum number of iterations.
@ XLAL_ELOSS
Loss of accuracy.
@ XLAL_EDOM
Input domain error.