Merge pull request #303 from LSSTDESC/fix_correlation

Disable full-sky shear correlation function calculation
LSSTDESC · Jan 24, 2018 · 1f68a8b · 1f68a8b
2 parents 9b011b0 + 4e76cbc
commit 1f68a8b
Show file tree

Hide file tree

Showing 4 changed files with 21 additions and 28 deletions.
diff --git a/include/ccl_error.h b/include/ccl_error.h
@@ -22,6 +22,7 @@ extern "C" {
 #define CCL_ERROR_HMF_INTERP 1035
 #define CCL_ERROR_PARAMETERS 1036
 #define CCL_ERROR_NU_INT 1037
+#define CCL_ERROR_NOT_IMPLEMENTED 1038
 
 typedef enum {
   CCL_ERROR_POLICY_EXIT = 0,

diff --git a/pyccl/ccl_wrap.c b/pyccl/ccl_wrap.c
@@ -19267,6 +19267,17 @@ SWIGINTERN PyObject *CCL_ERROR_NU_INT_swigconstant(PyObject *SWIGUNUSEDPARM(self
 }
 
 
+SWIGINTERN PyObject *CCL_ERROR_NOT_IMPLEMENTED_swigconstant(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *module;
+  PyObject *d;
+  if (!PyArg_ParseTuple(args,(char*)"O:swigconstant", &module)) return NULL;
+  d = PyModule_GetDict(module);
+  if (!d) return NULL;
+  SWIG_Python_SetConstant(d, "CCL_ERROR_NOT_IMPLEMENTED",SWIG_From_int((int)(1038)));
+  return SWIG_Py_Void();
+}
+
+
 SWIGINTERN PyObject *CCL_ERROR_POLICY_EXIT_swigconstant(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *module;
   PyObject *d;
@@ -20561,6 +20572,7 @@ static PyMethodDef SwigMethods[] = {
 	 { (char *)"CCL_ERROR_HMF_INTERP_swigconstant", CCL_ERROR_HMF_INTERP_swigconstant, METH_VARARGS, NULL},
 	 { (char *)"CCL_ERROR_PARAMETERS_swigconstant", CCL_ERROR_PARAMETERS_swigconstant, METH_VARARGS, NULL},
 	 { (char *)"CCL_ERROR_NU_INT_swigconstant", CCL_ERROR_NU_INT_swigconstant, METH_VARARGS, NULL},
+	 { (char *)"CCL_ERROR_NOT_IMPLEMENTED_swigconstant", CCL_ERROR_NOT_IMPLEMENTED_swigconstant, METH_VARARGS, NULL},
 	 { (char *)"CCL_ERROR_POLICY_EXIT_swigconstant", CCL_ERROR_POLICY_EXIT_swigconstant, METH_VARARGS, NULL},
 	 { (char *)"CCL_ERROR_POLICY_CONTINUE_swigconstant", CCL_ERROR_POLICY_CONTINUE_swigconstant, METH_VARARGS, NULL},
 	 { (char *)"raise_exception", _wrap_raise_exception, METH_VARARGS, (char *)"raise_exception(int err, char * msg)"},

diff --git a/pyccl/ccllib.py b/pyccl/ccllib.py
@@ -1432,6 +1432,9 @@ def __init__(self):
 _ccllib.CCL_ERROR_NU_INT_swigconstant(_ccllib)
 CCL_ERROR_NU_INT = _ccllib.CCL_ERROR_NU_INT
 
+_ccllib.CCL_ERROR_NOT_IMPLEMENTED_swigconstant(_ccllib)
+CCL_ERROR_NOT_IMPLEMENTED = _ccllib.CCL_ERROR_NOT_IMPLEMENTED
+
 _ccllib.CCL_ERROR_POLICY_EXIT_swigconstant(_ccllib)
 CCL_ERROR_POLICY_EXIT = _ccllib.CCL_ERROR_POLICY_EXIT
 

diff --git a/src/ccl_correlation.c b/src/ccl_correlation.c
@@ -291,32 +291,6 @@ static void ccl_compute_legendre_polynomial(int corr_type,int n_theta,double *th
       }
     }
   }
-  else if(corr_type==CCL_CORR_LP) {
-    for (int i=0;i<n_theta;i++) {
-      gsl_sf_legendre_Pl_array(ell_max,cos(theta[i]*M_PI/180),Pl_theta[i]);
-      for (int j=0;j<=ell_max;j++) {
-	Pl_theta[i][j]*=(2*j+1);
-      }
-    }
-  }
-  else if(corr_type==CCL_CORR_LM) {
-    for (int i=0;i<n_theta;i++) {
-      for (int j=0;j<=ell_max;j++) {
-	if(j>1e4) {///////////Some theta points thrown away for speed
-	  Pl_theta[i][j]=0;
-	  continue;
-	}
-	if (j<4) { 
-	  Pl_theta[i][j]=0;
-	  continue;
-	}
-	Pl_theta[i][j]=gsl_sf_legendre_Plm(j,4,cos(theta[i]*M_PI/180));
-	Pl_theta[i][j]*=(2*j+1)*pow(j,4);//approximate.. Using relation between bessel and legendre functions from Steibbens96.
-	for (k=-3;k<=4;k++)
-	  Pl_theta[i][j]/=(j+k);
-      }
-    }
-  }
 }
 
 /*--------ROUTINE: ccl_tracer_corr_legendre ------
@@ -333,6 +307,11 @@ static void ccl_tracer_corr_legendre(ccl_cosmology *cosmo,
   int i;
   double *l_arr,*cl_arr;
 
+  if(corr_type==CCL_CORR_LM || corr_type==CCL_CORR_LP){
+    *status=CCL_ERROR_NOT_IMPLEMENTED;
+    strcpy(cosmo->status_message,"ccl_correlation.c: CCL does not support full-sky xi+- calcuations.\nhttps://arxiv.org/abs/1702.05301 indicates flat-sky to be sufficient.\n");
+    return;
+  }
   l_arr=malloc((ELL_MAX_FFTLOG+1)*sizeof(double));
   if(l_arr==NULL) {
     *status=CCL_ERROR_MEMORY;
@@ -347,8 +326,6 @@ static void ccl_tracer_corr_legendre(ccl_cosmology *cosmo,
     return;
   }
 
-  if(corr_type==CCL_CORR_LM)
-    printf("WARNING: legendre sum for xi- is still not correctly implemented.\n");
 
   //Interpolate input Cl into 
   SplPar *cl_spl=ccl_spline_init(n_ell,ell,cls,cls[0],0);