modified fftlog
/* This code is FFTLog, which is described in arXiv:astro-ph/9905191 */
/* Computes the Gamma function using the Lanczos approximation */
-static double complex gamma_fftlog(double complex z) {
+static double complex gamma_fftlog(double complex z)
+{
/* Lanczos coefficients for g = 7 */
static double p[] = {
0.99999999999980993227684700473478,
@@ -33,22 +34,26 @@ static double complex gamma_fftlog(double complex z) {
return sqrt(2*M_PI) * cpow(t, z+0.5) * cexp(-t) * x;
}
-static double complex polar (double r, double phi) {
+static double complex polar (double r, double phi)
+{
return (r*cos(phi) +I*(r*sin(phi)));
}
-static double complex lngamma_fftlog(double complex z) {
+static double complex lngamma_fftlog(double complex z)
+{
return clog(gamma_fftlog(z));
}
-static void lngamma_4(double x, double y, double* lnr, double* arg) {
+static void lngamma_4(double x, double y, double* lnr, double* arg)
+{
double complex w = lngamma_fftlog(x+y*I);
if(lnr) *lnr = creal(w);
if(arg) *arg = cimag(w);
}
-static double goodkr(int N, double mu, double q, double L, double kr) {
+static double goodkr(int N, double mu, double q, double L, double kr)
+{
double xp = (mu+1+q)/2;
double xm = (mu+1-q)/2;
double y = M_PI*N/(2*L);
@@ -62,7 +67,8 @@ static double goodkr(int N, double mu, double q, double L, double kr) {
return kr;
}
-void compute_u_coefficients(int N, double mu, double q, double L, double kcrc, double complex u[]) {
+void compute_u_coefficients(int N, double mu, double q, double L, double kcrc, double complex u[])
+{
double y = M_PI/L;
double k0r0 = kcrc * exp(-L);
double t = -2*y*log(k0r0/2);
@@ -132,7 +138,9 @@ void fht(int N, const double r[], const double complex a[], double k[], double c
free(ulocal);
}
-void fftlog_ComputeXiLM(double l, double m, int N, const double k[], const double pk[], double r[], double xi[]) {
+void fftlog_ComputeXiLM(double l, double m, int N, const double k[], const double pk[],
+ double r[], double xi[])
+{
double complex* a = malloc(sizeof(complex double)*N);
double complex* b = malloc(sizeof(complex double)*N);
@@ -147,11 +155,13 @@ void fftlog_ComputeXiLM(double l, double m, int N, const double k[], const doubl
free(b);
}
-void pk2xi(int N, const double k[], const double pk[], double r[], double xi[]) {
+void pk2xi(int N, const double k[], const double pk[], double r[], double xi[])
+{
fftlog_ComputeXiLM(0, 2, N, k, pk, r, xi);
}
-void xi2pk(int N, const double r[], const double xi[], double k[], double pk[]) {
+void xi2pk(int N, const double r[], const double xi[], double k[], double pk[])
+{
static const double TwoPiCubed = 8*M_PI*M_PI*M_PI;
fftlog_ComputeXiLM(0, 2, N, r, xi, k, pk);
for(int j = 0; j < N; j++)
