Merge pull request #356 from LSSTDESC/correlation_simplify

simplified use of fftlog for correlation functions

include/fftlog.h

@@ -42,7 +42,13 @@ void xi2pk(int N,  const double r[],  const double xi[], double k[], double pk[]
 void fftlog_ComputeXiLM(double l, double m, int N, const double k[],  const double pk[],
 			double r[], double xi[]);
 
-
+/* Compute the function
+ *   \xi_\alpha(\theta) = \int_0^\infty \frac{d\ell}{2\pi} \ell J_\alpha(\ell\theta) C_\ell
+ * The input l-values must be logarithmically spaced.  The
+ * resulting xi_alpha(th) will be evaluated at the dual th-values
+ *   th[0] = 1/l[N-1], ..., th[N-1] = 1/l[0]. */
+void fftlog_ComputeXi2D(double bessel_order,int N,const double l[],const double cl[],
+			double th[], double xi[]);
 #include <complex.h>
 
 /* Compute the discrete Hankel transform of the function a(r).  See the FFTLog

src/ccl_correlation.c

@@ -126,19 +126,12 @@ static void ccl_tracer_corr_fftlog(ccl_cosmology *cosmo,
     th_arr[i]=0;
   //Although set here to 0, theta is modified by FFTlog to obtain the correlation at ~1/l
 
-  // FFTlog uses spherical bessel functions, j_n, but for projected
-  // correlations we need bessel functions of first order, J_n.
-  // To compensate for the difference, we use the relation
-  // j_n(x) = sqrt(Pi/2x)J_{n+1/2}(x)
-  // J_{m}(x) = sqrt(2x/Pi) j_{m-1/2}(x)
   int i_bessel=0;
   if(corr_type==CCL_CORR_GG) i_bessel=0;
   if(corr_type==CCL_CORR_GL) i_bessel=2;
   if(corr_type==CCL_CORR_LP) i_bessel=0;
   if(corr_type==CCL_CORR_LM) i_bessel=4;
-  fftlog_ComputeXiLM(i_bessel-0.5,1.5,N_ELL_FFTLOG,l_arr,cl_arr,th_arr,wth_arr);
-  for(i=0;i<N_ELL_FFTLOG;i++)
-    wth_arr[i]*=sqrt(th_arr[i]*2.0*M_PI);
+  fftlog_ComputeXi2D(i_bessel,N_ELL_FFTLOG,l_arr,cl_arr,th_arr,wth_arr);
 
   // Interpolate to output values of theta
   SplPar *wth_spl=ccl_spline_init(N_ELL_FFTLOG,th_arr,wth_arr,wth_arr[0],0);

src/fftlog.c

@@ -138,6 +138,22 @@ void fht(int N, const double r[], const double complex a[], double k[], double c
   free(ulocal);
 }
 
+void fftlog_ComputeXi2D(double bessel_order,int N,const double l[],const double cl[],
+			double th[], double xi[])
+{
+  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[])
 {
@@ -148,8 +164,7 @@ void fftlog_ComputeXiLM(double l, double m, int N, const double k[], const doubl
     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]);
-    xi[i] = creal(pow(2*M_PI*r[i], -1.5) * b[i]);
+    xi[i] = creal(pow(2*M_PI*r[i], -(m-0.5)) * b[i]);
 
   free(a);
   free(b);