fftlog.c File Reference

#include <stdlib.h>
#include <math.h>
#include <fftw3.h>
#include "fftlog.h"

Include dependency graph for fftlog.c:

Go to the source code of this file.

Functions
static double complex	gamma_fftlog (double complex z)
static double complex	polar (double r, double phi)
static double complex	lngamma_fftlog (double complex z)
static void	lngamma_4 (double x, double y, double *lnr, double *arg)
static double	goodkr (int N, double mu, double q, double L, double kr)
void	compute_u_coefficients (int N, double mu, double q, double L, double kcrc, double complex u[])
void	fht (int N, const double r[], const double complex a[], double k[], double complex b[], double mu, double q, double kcrc, int noring, double complex *u)
void	fftlog_ComputeXi2D (double bessel_order, int N, const double l[], const double cl[], double th[], double xi[])
void	fftlog_ComputeXiLM (double l, double m, 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[])
void	xi2pk (int N, const double r[], const double xi[], double k[], double pk[])

Function Documentation

void compute_u_coefficients	(	int	N,
		double	mu,
		double	q,
		double	L,
		double	kcrc,
		double complex	u[]
	)

Definition at line 73 of file fftlog.c.

References lngamma_4(), m, M_PI, polar(), and x.

Referenced by fht().

{
  double y = M_PI/L;
  double k0r0 = kcrc * exp(-L);
  double t = -2*y*log(k0r0/2);
  
  if(q == 0) {
    double x = (mu+1)/2;
    double lnr, phi;
    for(int m = 0; m <= N/2; m++) {
      lngamma_4(x, m*y, &lnr, &phi);
      u[m] = polar(1.0,m*t + 2*phi);
    }
  }
  else {
    double xp = (mu+1+q)/2;
    double xm = (mu+1-q)/2;
    double lnrp, phip, lnrm, phim;
    for(int m = 0; m <= N/2; m++) {
      lngamma_4(xp, m*y, &lnrp, &phip);
      lngamma_4(xm, m*y, &lnrm, &phim);
      u[m] = polar(exp(q*log(2) + lnrp - lnrm), m*t + phip - phim);
    }
  }
  
  for(int m = N/2+1; m < N; m++)
    u[m] = conj(u[N-m]);
  if((N % 2) == 0)
    u[N/2] = (creal(u[N/2]) + I*0.0);
}

void fftlog_ComputeXi2D	(	double	bessel_order,
		int	N,
		const double	l[],
		const double	cl[],
		double	th[],
		double	xi[]
	)

Definition at line 144 of file fftlog.c.

References fht(), and M_PI.

Referenced by ccl_tracer_corr_fftlog().

{
  double complex* a = malloc(sizeof(complex double)*N);
  double complex* b = malloc(sizeof(complex double)*N);
  
  for(int i=0;i<N;i++)
    a[i]=l[i]*cl[i];
  fht(N,l,a,th,b,bessel_order,0,1,1,NULL);
  for(int i=0;i<N;i++)
    xi[i]=creal(b[i]/(2*M_PI*th[i]));
  
  free(a);
  free(b);
}

void fftlog_ComputeXiLM	(	double	l,
		double	m,
		int	N,
		const double	k[],
		const double	pk[],
		double	r[],
		double	xi[]
	)

Definition at line 160 of file fftlog.c.

References fht(), M_PI, and pow().

Referenced by pk2xi(), and xi2pk().

{
  double complex* a = malloc(sizeof(complex double)*N);
  double complex* b = malloc(sizeof(complex double)*N);
  
  for(int i = 0; i < N; i++)
    a[i] = pow(k[i], m - 0.5) * pk[i];
  fht(N, k, a, r, b, l + 0.5, 0, 1, 1, NULL);
  for(int i = 0; i < N; i++)
    xi[i] = creal(pow(2*M_PI*r[i], -(m-0.5)) * b[i]);
  
  free(a);
  free(b);
}

void fht	(	int	N,
		const double	r[],
		const double complex	a[],
		double	k[],
		double complex	b[],
		double	mu,
		double	q,
		double	kcrc,
		int	noring,
		double complex *	u
	)

Definition at line 104 of file fftlog.c.

References compute_u_coefficients(), goodkr(), and m.

Referenced by fftlog_ComputeXi2D(), and fftlog_ComputeXiLM().

{
  double L = log(r[N-1]/r[0]) * N/(N-1.);
  double complex* ulocal = NULL;
  if(u == NULL) {
    if(noring)
      kcrc = goodkr(N, mu, q, L, kcrc);
    ulocal = malloc (sizeof(complex double)*N); 
    compute_u_coefficients(N, mu, q, L, kcrc, ulocal);
    u = ulocal;
  }
  
  /* Compute the convolution b = a*u using FFTs */
  fftw_plan forward_plan = fftw_plan_dft_1d(N, (fftw_complex*) a, (fftw_complex*) b,  -1, FFTW_ESTIMATE);
  fftw_plan reverse_plan = fftw_plan_dft_1d(N, (fftw_complex*) b, (fftw_complex*) b, +1, FFTW_ESTIMATE);
  fftw_execute(forward_plan);
  for(int m = 0; m < N; m++)
    b[m] *= u[m] / (double)(N);       // divide by N since FFTW doesn't normalize the inverse FFT
  fftw_execute(reverse_plan);
  fftw_destroy_plan(forward_plan);
  fftw_destroy_plan(reverse_plan);
  
  /* Reverse b array */
  double complex tmp;
  for(int n = 0; n < N/2; n++) {
    tmp = b[n];
    b[n] = b[N-n-1];
    b[N-n-1] = tmp;
  }
  
  /* Compute k's corresponding to input r's */
  double k0r0 = kcrc * exp(-L);
  k[0] = k0r0/r[0];
  for(int n = 1; n < N; n++)
    k[n] = k[0] * exp(n*L/N);
  
  free(ulocal);
}

static double complex gamma_fftlog ( double complex  z )

static

Definition at line 15 of file fftlog.c.

References M_PI, p, sqrt(), and x.

Referenced by lngamma_fftlog().

{
  /* Lanczos coefficients for g = 7 */
  static double p[] = {
    0.99999999999980993227684700473478,
    676.520368121885098567009190444019,
    -1259.13921672240287047156078755283,
    771.3234287776530788486528258894,
    -176.61502916214059906584551354,
    12.507343278686904814458936853,
    -0.13857109526572011689554707,
    9.984369578019570859563e-6,
    1.50563273514931155834e-7
  };
  
  if(creal(z) < 0.5)
    return M_PI / (sin(M_PI*z)*gamma_fftlog(1. - z));
  z -= 1;
  double complex x = p[0];
  for(int n = 1; n < 9; n++)
    x += p[n] / (z + (double)(n));
  double complex t = z + 7.5;
  return sqrt(2*M_PI) * cpow(t, z+0.5) * cexp(-t) * x;
}

static double goodkr ( int  N, double  mu, double  q, double  L, double  kr  )

static

Definition at line 58 of file fftlog.c.

References lngamma_4(), and M_PI.

Referenced by fht().

{
  double xp = (mu+1+q)/2;
  double xm = (mu+1-q)/2;
  double y = M_PI*N/(2*L);
  double lnr, argm, argp;
  lngamma_4(xp, y, &lnr, &argp);
  lngamma_4(xm, y, &lnr, &argm);
  double arg = log(2/kr) * N/L + (argp + argm)/M_PI;
  double iarg = round(arg);
  if(arg != iarg)
    kr *= exp((arg - iarg)*L/N);
  return kr;
}

static void lngamma_4 ( double  x, double  y, double *  lnr, double *  arg  )

static

Definition at line 51 of file fftlog.c.

References lngamma_fftlog(), and w.

Referenced by compute_u_coefficients(), and goodkr().

{
  double complex w = lngamma_fftlog(x+y*I);
  if(lnr) *lnr = creal(w);
  if(arg) *arg = cimag(w);
}

static double complex lngamma_fftlog ( double complex  z )

static

Definition at line 45 of file fftlog.c.

References Catch::clog(), and gamma_fftlog().

Referenced by lngamma_4().

{
  return clog(gamma_fftlog(z));
}

void pk2xi	(	int	N,
		const double	k[],
		const double	pk[],
		double	r[],
		double	xi[]
	)

Definition at line 176 of file fftlog.c.

References fftlog_ComputeXiLM().

Referenced by ccl_correlation_3d().

{
  fftlog_ComputeXiLM(0, 2, N, k, pk, r, xi);
}

static double complex polar ( double  r, double  phi  )

static

Definition at line 40 of file fftlog.c.

Referenced by compute_u_coefficients().

{
  return (r*cos(phi) +I*(r*sin(phi)));
}

void xi2pk	(	int	N,
		const double	r[],
		const double	xi[],
		double	k[],
		double	pk[]
	)

Definition at line 181 of file fftlog.c.

References fftlog_ComputeXiLM(), and M_PI.

{
  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++)
    pk[j] *= TwoPiCubed;
}