diff --git a/benchmarks/data/curved_sky_cocoa.npz b/benchmarks/data/curved_sky_cocoa.npz
new file mode 100644
index 000000000..2bf23e001
Binary files /dev/null and b/benchmarks/data/curved_sky_cocoa.npz differ
diff --git a/benchmarks/test_correlation.py b/benchmarks/test_correlation.py
index cdc74efb3..8ded491dc 100644
--- a/benchmarks/test_correlation.py
+++ b/benchmarks/test_correlation.py
@@ -229,3 +229,47 @@ def test_xi(set_up, corr_method, t1, t2, bm, er, kind, pref):
 
     xi *= pref
     assert np.all(np.fabs(xi - bms[bm]) < ers[er] * errfac)
+
+
+def test_correlation_curvedsky():
+
+    data = np.load(os.path.dirname(__file__) +
+                   "/data/curved_sky_cocoa.npz", allow_pickle=True)
+    COSMO = ccl.Cosmology(**data['cosmo_setup'][()])
+
+    nsrcs = 5
+    nlens = 10
+
+    theta_b = data['theta_b']
+    theta_u = theta_b[1:]
+    theta_l = theta_b[:-1]
+    ell_hr = np.arange(60000)
+
+    inp = dict(ell=ell_hr, theta=theta_l, theta_max=theta_u,
+               method='Legendre')
+
+    xip, xim, gammat, wtheta = [], [], [], []
+
+    for i in range(nsrcs*(nsrcs+1)//2):
+        C_ss_i = np.interp(ell_hr, data['ell'], data['C_ss'][i])
+        xip.append(ccl.correlation(COSMO, C_ell=C_ss_i, type='GG+', **inp))
+        xim.append(ccl.correlation(COSMO, C_ell=C_ss_i, type='GG-', **inp))
+
+    for i in range(nsrcs*nlens):
+        C_gs_i = np.interp(ell_hr, data['ell'], data['C_gs'][i])
+        gammat.append(ccl.correlation(COSMO, C_ell=C_gs_i, type='NG', **inp))
+
+    for i in range(nlens):
+        C_gg_i = np.interp(ell_hr, data['ell'], data['C_gg'][i])
+        wtheta.append(ccl.correlation(COSMO, C_ell=C_gg_i, type='NN', **inp))
+
+    xip = np.concatenate(xip)
+    xim = np.concatenate(xim)
+    gammat = np.concatenate(gammat)
+    wtheta = np.concatenate(wtheta)
+    dv_ccl = np.concatenate((xip, xim, gammat, wtheta))
+    dv_cca = data['dv_cocoa']
+
+    diff = dv_ccl/dv_cca - 1
+
+    assert np.max(np.abs(diff)) < 0.01
diff --git a/include/ccl_correlation.h b/include/ccl_correlation.h
index c5ab8b7ca..8603824f3 100644
--- a/include/ccl_correlation.h
+++ b/include/ccl_correlation.h
@@ -41,6 +41,36 @@ void ccl_correlation(ccl_cosmology *cosmo,
 		     int corr_type,int do_taper_cl,double *taper_cl_limits,int flag_method,
 		     int *status);
 
+
+             /**
+ * Computes the correlation function (wrapper)
+ * @param cosmo :Cosmological parameters
+ * @param n_ell : number of multipoles in the input power spectrum
+ * @param ell : multipoles at which the power spectrum is evaluated
+ * @param cls : input power spectrum
+ * @param n_theta : number of output values of the separation angle (theta)
+ * @param theta_min : Min/left  edge of the separation bin in degrees.
+ * @param theta_max : Max/right edge of the separation bin in degrees.
+ * @param wtheta : the values of the correlation function at the angles above will be returned in this array, which should be pre-allocated
+ * @param do_taper_cl :
+ * @param taper_cl_limits
+ * @param flag_method : method to compute the correlation function. Choose between:
+ *  - CCL_CORR_FFTLOG : fast integration with FFTLog
+ *  - CCL_CORR_BESSEL : direct integration over the Bessel function
+ *  - CCL_CORR_LGNDRE : brute-force sum over legendre polynomials
+ * @param corr_type : type of correlation function. Choose between:
+ *  - CCL_CORR_GG : spin0-spin0
+ *  - CCL_CORR_GL : spin0-spin2
+ *  - CCL_CORR_LP : spin2-spin2 (xi+)
+ *  - CCL_CORR_LM : spin2-spin2 (xi-)
+ * Currently supported spin-0 fields are number counts and CMB lensing. The only spin-2 is currently shear.
+ */
+void ccl_correlation_binned(ccl_cosmology *cosmo,
+		     int n_ell,double *ell,double *cls,
+		     int n_theta,double *theta_min,double *theta_max,double *wtheta,
+		     int corr_type,int do_taper_cl,double *taper_cl_limits,int flag_method,
+		     int *status);
+
 /**
  * Computes the 3dcorrelation function (wrapper)
  * @param cosmo :Cosmological parameters
diff --git a/pyccl/ccl_correlation.i b/pyccl/ccl_correlation.i
index 4e975a00b..3db796d02 100644
--- a/pyccl/ccl_correlation.i
+++ b/pyccl/ccl_correlation.i
@@ -11,6 +11,8 @@
     (double* larr, int nlarr),
     (double* clarr, int nclarr),
     (double* theta, int nt),
+    (double* theta_min, int nt),
+    (double* theta_max, int ntmax),
     (double* r, int nr),
     (double* s, int ns),
     (double* sig, int nsig)}
@@ -27,6 +29,17 @@
         raise CCLError("Input shape for `theta` must match `(nout,)`!")
 %}
 
+%feature("pythonprepend") correlation_vec_binned %{
+    if numpy.shape(larr) != numpy.shape(clarr):
+        raise CCLError("Input shape for `larr` must match `clarr`!")
+
+    if numpy.shape(theta_min) != (nout,):
+        raise CCLError("Input shape for `theta_min` must match `(nout,)`!")
+    
+    if numpy.shape(theta_max) != (nout,):
+        raise CCLError("Input shape for `theta_max` must match `(nout,)`!")
+%}
+
 %feature("pythonprepend") correlation_3d_vec %{
     if numpy.shape(r) != (nxi,):
         raise CCLError("Input shape for `r` must match `(nxi,)`!")
@@ -63,6 +76,30 @@ void correlation_vec(ccl_cosmology *cosmo, double* larr, int nlarr,
         output, corr_type, 0, NULL, method, status);
 }
 
+void correlation_vec_binned(ccl_cosmology *cosmo,
+                            double* larr, int nlarr,
+                            double* clarr, int nclarr,
+                            double* theta_min, int nt,
+                            double* theta_max, int ntmax,
+                            int corr_type, int method,
+                            int nout,
+                            double* output,
+                            int *status) {
+    if (ntmax != nt) {
+        *status = CCL_ERROR_INCONSISTENT;
+        return;
+    }
+
+    if (nout != nt) {
+        *status = CCL_ERROR_INCONSISTENT;
+        return;
+    }
+
+    ccl_correlation_binned(
+        cosmo, nlarr, larr, clarr, nt, theta_min, theta_max,
+        output, corr_type, 0, NULL, method, status);
+}
+
 void correlation_3d_vec(ccl_cosmology *cosmo,ccl_f2d_t *psp,
                         double a, double* r, int nr,
                         int nxi, double* xi, int *status) {
diff --git a/pyccl/correlations.py b/pyccl/correlations.py
index 8102239f3..e93b66349 100644
--- a/pyccl/correlations.py
+++ b/pyccl/correlations.py
@@ -42,7 +42,8 @@ class CorrelationTypes(Enum):
 }
 
 
-def correlation(cosmo, *, ell, C_ell, theta, type='NN', method='fftlog'):
+def correlation(cosmo, *, ell, C_ell, theta,
+                type='NN', method='fftlog', theta_max=None):
     r"""Compute the angular correlation function.
 
     .. math::
@@ -71,6 +72,12 @@ def correlation(cosmo, *, ell, C_ell, theta, type='NN', method='fftlog'):
         * :math:`s_a=2`, :math:`s_b=0` e.g. galaxy-shear, and :math:`\kappa`-shear
         * :math:`s_a=s_b=2` e.g. shear-shear.
 
+    Bin-averaging, where we predict the average correlation within a given theta-bin,
+    is also implemented. This is only available for the Legendre method, as the 
+    bin-averaging can be done analytically in that case. The method we use for the
+    NN, NG, and GG+/- correlations follow Eq 65, B1, and B5, respectively, from 
+    arxiv:2012.08568
+
     .. note::
         For scales smaller than :math:`\sim 0.1^{\circ}`, the input power
         spectrum should be sampled to sufficienly high :math:`\ell` to ensure
@@ -87,7 +94,9 @@ def correlation(cosmo, *, ell, C_ell, theta, type='NN', method='fftlog'):
                           spectrum.
         C_ell (array): Input angular power spectrum.
         theta (:obj:`float` or `array`): Angular separation(s) at which to
-            calculate the angular correlation function (in degrees).
+            calculate the angular correlation function (in degrees). If theta_max
+            is passed then this is interpreted as the minimum/left edge in a given
+            separation bin.
         type (:obj:`str`): Type of correlation function. Choices: ``'NN'`` (0x0),
             ``'NG'`` (0x2), ``'GG+'`` (2x2, :math:`\xi_+`),
             ``'GG-'`` (2x2, :math:`\xi_-`), where numbers refer to the spins
@@ -99,6 +108,9 @@ def correlation(cosmo, *, ell, C_ell, theta, type='NN', method='fftlog'):
             Choices: ``'Bessel'`` (direct integration over Bessel function),
             ``'FFTLog'`` (fast integration with FFTLog), ``'Legendre'``
             (brute-force sum over Legendre polynomials).
+        theta_max (:obj:`float` or `array1): Maximum angular separation(s) for the given
+            separation bin (in degrees). If passed, then the `theta` input is interpreted
+            as a minimum/left edge of the bin.
 
     Returns:
         (:obj:`float` or `array`): Value(s) of the correlation function at the
@@ -119,15 +131,40 @@ def correlation(cosmo, *, ell, C_ell, theta, type='NN', method='fftlog'):
     if scalar := isinstance(theta, (int, float)):
         theta = np.array([theta, ])
 
+        if theta_max is not None:
+            if isinstance(theta_max, (int, float)):
+                theta_max = np.array([theta_max, ])
+
+            assert theta.size == theta_max.size, \
+                (f"theta and theta_max have different sizes"
+                 f"({theta.size} != {theta_max.size})")
+
+            assert not np.allclose(theta, theta_max), \
+                "theta_min cannot be the same as theta_max"
+
     if np.all(np.array(C_ell) == 0):
         # short-cut and also avoid integration errors
         wth = np.zeros_like(theta)
     else:
-        # Call correlation function
-        wth, status = lib.correlation_vec(cosmo, ell, C_ell, theta,
-                                          correlation_types[type],
-                                          correlation_methods[method],
-                                          len(theta), status)
+        if theta_max is None:
+            wth, status = lib.correlation_vec(
+                cosmo, ell, C_ell, theta,
+                correlation_types[type],
+                correlation_methods[method],
+                len(theta), status)
+        else:
+            wth, status = lib.correlation_vec_binned(
+                cosmo,
+                ell,
+                C_ell,
+                theta,
+                theta_max,
+                correlation_types[type],
+                correlation_methods[method],
+                len(theta),
+                status,
+            )
+
     check(status, cosmo_in)
     if scalar:
         return wth[0]
diff --git a/pyccl/tests/test_correlations.py b/pyccl/tests/test_correlations.py
index c421ec5b7..615d4105d 100644
--- a/pyccl/tests/test_correlations.py
+++ b/pyccl/tests/test_correlations.py
@@ -29,6 +29,20 @@ def test_correlation_smoke(method):
         assert np.all(np.isfinite(corr))
         assert np.shape(corr) == np.shape(tval)
 
+    if method == 'legendre':
+        for tval in [t_arr, t_lst, t_scl, t_int]:
+
+            if isinstance(tval, list):
+                tmax = [t*2 for t in tval]
+            else:
+                tmax = tval*2
+
+            corr = ccl.correlation(
+                COSMO, ell=ell, C_ell=cl, theta=tval, type='NN', method=method,
+                theta_max=tmax)
+            assert np.all(np.isfinite(corr))
+            assert np.shape(corr) == np.shape(tval)
+
 
 @pytest.mark.parametrize(
     'rval',
diff --git a/src/ccl_correlation.c b/src/ccl_correlation.c
index 626cede74..330bf5af8 100644
--- a/src/ccl_correlation.c
+++ b/src/ccl_correlation.c
@@ -270,26 +270,223 @@ static void ccl_tracer_corr_bessel(ccl_cosmology *cosmo,
 TASK: Compute input factor for ccl_tracer_corr_legendre
 INPUT: tracer 1, tracer 2, i_bessel, theta array, n_theta, L_max, output Pl_theta
  */
-static void ccl_compute_legendre_polynomial(int corr_type,double theta,int ell_max,double *Pl_theta)
+static void ccl_compute_legendre_polynomial(int corr_type,
+                                            double theta,
+                                            int ell_max,
+                                            double *Pl_theta)
 {
-  int j;
-  double cth=cos(theta*M_PI/180);
-
-  //Initialize Pl_theta
-  for (j=0;j<=ell_max;j++)
-      Pl_theta[j]=0.;
-
-  if(corr_type==CCL_CORR_GG) {
-    gsl_sf_legendre_Pl_array(ell_max,cth,Pl_theta);
-    for (j=0;j<=ell_max;j++)
-      Pl_theta[j]*=(2*j+1);
-  }
-  else if(corr_type==CCL_CORR_GL) {
-    for (j=2;j<=ell_max;j++) {//https://arxiv.org/pdf/1007.4809.pdf
-      Pl_theta[j]=gsl_sf_legendre_Plm(j,2,cth);
-      Pl_theta[j]*=(2*j+1.)/((j+0.)*(j+1.));
+  int ell;
+  double x = cos(theta * M_PI / 180.0);
+  const double eps = 1e-14;
+
+  /* Avoid exact endpoints because G^{+/-}_{l,2} contains 1/(1-x^2) */
+  if (x >  1.0 - eps) x =  1.0 - eps;
+  if (x < -1.0 + eps) x = -1.0 + eps;
+
+  for (ell = 0; ell <= ell_max; ell++)
+    Pl_theta[ell] = 0.0;
+
+  if (corr_type == CCL_CORR_GG) {
+    gsl_sf_legendre_Pl_array(ell_max, x, Pl_theta);
+    for (ell = 0; ell <= ell_max; ell++)
+      Pl_theta[ell] *= (2.0 * ell + 1.0);
+  }
+
+  else if (corr_type == CCL_CORR_GL) {
+    gsl_sf_legendre_Plm_array(ell_max, 2, x, &(Pl_theta[2]));
+    for (ell = 2; ell <= ell_max; ell++)
+      Pl_theta[ell] *= (2.0 * ell + 1.0) / (ell * (ell + 1.0));
+  }
+
+  else if ((corr_type == CCL_CORR_LP) || (corr_type == CCL_CORR_LM)) {
+    double *P2;
+    double omx2 = 1.0 - x * x; //1 - x^2, with x = cos(theta)
+    int m       = 2.0;
+
+    P2 = malloc((ell_max + 1) * sizeof(double));
+    if (P2 == NULL)
+      return ;
+
+    for (ell = 0; ell <= ell_max; ell++)
+      P2[ell] = 0.0;
+
+    /* GSL fills P_l^2(x) for l=2..ell_max into P2[2..ell_max] */
+    gsl_sf_legendre_Plm_array(ell_max, (int) m, x, &(P2[2]));
+
+    for (ell = 2; ell <= ell_max; ell++) {
+      double Pl2   = P2[ell];
+      double Plm12 = (ell > 2) ? P2[ell - 1] : 0.0;
+      double Gplus, Gminus;
+      double pref;
+
+      /* Eq. 4.19 from https://arxiv.org/pdf/astro-ph/9609149 */
+      Gplus  = -(((ell - m * m) / omx2) + 0.5 * ell * (ell - 1.0)) * Pl2 + ((ell + m) * x / omx2) * Plm12;
+      Gminus = m * (((ell - 1.0) * x / omx2) * Pl2 - ((ell + m) / omx2) * Plm12);
+
+      /* ccl_tracer_corr_legendre divides the sum by 4*pi afterward,
+       * so we store an extra factor of 2 here to reproduce Eq. 19 in https://arxiv.org/abs/2105.13548
+       */
+      pref = 2.0 * (2.0 * ell + 1.0) / ((double) ell * ell * (ell + 1.0) * (ell + 1.0));
+
+      if (corr_type == CCL_CORR_LP)
+        Pl_theta[ell] = pref * (Gplus + Gminus);
+      else
+        Pl_theta[ell] = pref * (Gplus - Gminus);
     }
+
+    free(P2);
+  }
+}
+
+
+/* A convenience function for the binning kernels */
+static inline double dpl_from_array(int n, double x, const double *P)
+{
+  if (n == 0) {
+    return 0.0;
   }
+  return ((double)n) * (P[n - 1] - x * P[n]) / (1.0 - x * x);
+}
+
+static void ccl_compute_legendre_polynomial_binned(int corr_type,
+                                                   double theta_min,
+                                                   double theta_max,
+                                                   int ell_max,
+                                                   double *Pl_theta)
+{
+  int ell;
+  const double eps = 1e-14;
+
+  /* x_hi corresponds to the lower angular edge theta_min,
+   * x_lo corresponds to the upper angular edge theta_max.
+   */
+  double x_hi = cos(theta_min * M_PI / 180.0);
+  double x_lo = cos(theta_max * M_PI / 180.0);
+  double dx = x_hi - x_lo;
+
+  double *P_hi = NULL;
+  double *P_lo = NULL;
+
+  //Initialize with zeros
+  for (ell = 0; ell <= ell_max; ell++)
+    Pl_theta[ell] = 0.0;
+
+  /* Clamp endpoints slightly away from |x|=1 for derivative evaluations. */
+  if (x_hi >  1.0 - eps) x_hi =  1.0 - eps;
+  if (x_hi < -1.0 + eps) x_hi = -1.0 + eps;
+  if (x_lo >  1.0 - eps) x_lo =  1.0 - eps;
+  if (x_lo < -1.0 + eps) x_lo = -1.0 + eps;
+
+  dx = x_hi - x_lo;
+
+  /* Need up to ell_max + 1 because some expressions involve P_{ell+1}. */
+  P_hi = malloc((ell_max + 2) * sizeof(double));
+  P_lo = malloc((ell_max + 2) * sizeof(double));
+  if ((P_hi == NULL) || (P_lo == NULL)) {
+    free(P_hi);
+    free(P_lo);
+    return;
+  }
+
+  gsl_sf_legendre_Pl_array(ell_max + 1, x_hi, P_hi);
+  gsl_sf_legendre_Pl_array(ell_max + 1, x_lo, P_lo);
+
+  if (corr_type == CCL_CORR_GG) {
+    /* Eq. 20 of https://arxiv.org/abs/2105.13548:
+     * (2 ell + 1) * \bar{P_ell} =
+     *   ([P_{ell+1} - P_{ell-1}]_{x_lo}^{x_hi}) / (x_hi - x_lo)
+     */
+    Pl_theta[0] = 1.0;
+    for (ell = 1; ell <= ell_max; ell++) {
+      Pl_theta[ell] =
+        ((P_hi[ell + 1] - P_hi[ell - 1]) -
+         (P_lo[ell + 1] - P_lo[ell - 1])) / dx;
+    }
+  }
+
+  else if (corr_type == CCL_CORR_GL) {
+    /* Eq. B2 of https://arxiv.org/abs/2012.08568:
+     * average of P_ell^2 over the angular bin.
+     * CCL stores:
+     *   (2 ell + 1) / [ell (ell + 1)] * <P_ell^2>_bin
+     */
+    for (ell = 2; ell <= ell_max; ell++) {
+      double Ibin_hi, Ibin_lo, avgP2;
+
+      Ibin_hi =
+        (((double)ell + 2.0 / (2.0 * ell + 1.0)) * P_hi[ell - 1]) +
+        ((2.0 - ell) * x_hi * P_hi[ell]) -
+        (2.0 / (2.0 * ell + 1.0)) * P_hi[ell + 1];
+
+      Ibin_lo =
+        (((double)ell + 2.0 / (2.0 * ell + 1.0)) * P_lo[ell - 1]) +
+        ((2.0 - ell) * x_lo * P_lo[ell]) -
+        (2.0 / (2.0 * ell + 1.0)) * P_lo[ell + 1];
+
+      avgP2 = (Ibin_hi - Ibin_lo) / dx;
+
+      Pl_theta[ell] =
+        ((2.0 * ell + 1.0) / ((double)ell * (ell + 1.0))) * avgP2;
+    }
+  }
+
+  else if ((corr_type == CCL_CORR_LP) || (corr_type == CCL_CORR_LM)) {
+    /* Eq. B5:
+     * average of G^+_{ell,2} +/- G^-_{ell,2} over the angular bin.
+     * CCL stores an extra factor of 2 so that the later 1/(4 pi)
+     * gives the Eq. 19 normalization for xi_+ / xi_-.
+     */
+    for (ell = 2; ell <= ell_max; ell++) {
+      double dPl_hi   = dpl_from_array(ell,     x_hi, P_hi);
+      double dPl_lo   = dpl_from_array(ell,     x_lo, P_lo);
+      double dPlm1_hi = dpl_from_array(ell - 1, x_hi, P_hi);
+      double dPlm1_lo = dpl_from_array(ell - 1, x_lo, P_lo);
+
+      double term_hi, term_lo;
+      double avgGpm;
+      double pref;
+
+      term_hi =
+        -(0.5 * ell * (ell - 1.0)) *
+          (ell + 2.0 / (2.0 * ell + 1.0)) * P_hi[ell - 1]
+        -0.5 * ell * (ell - 1.0) * (2.0 - ell) * (x_hi * P_hi[ell])
+        +(ell * (ell - 1.0) / (2.0 * ell + 1.0)) * P_hi[ell + 1]
+        +(4.0 - ell) * dPl_hi
+        +(ell + 2.0) * (x_hi * dPlm1_hi - P_hi[ell - 1]);
+
+      term_lo =
+        -(0.5 * ell * (ell - 1.0)) *
+          (ell + 2.0 / (2.0 * ell + 1.0)) * P_lo[ell - 1]
+        -0.5 * ell * (ell - 1.0) * (2.0 - ell) * (x_lo * P_lo[ell])
+        +(ell * (ell - 1.0) / (2.0 * ell + 1.0)) * P_lo[ell + 1]
+        +(4.0 - ell) * dPl_lo
+        +(ell + 2.0) * (x_lo * dPlm1_lo - P_lo[ell - 1]);
+
+      if (corr_type == CCL_CORR_LP) {
+        term_hi +=  2.0 * (ell - 1.0) * (x_hi * dPl_hi - P_hi[ell])
+                   -2.0 * (ell + 2.0) * dPlm1_hi;
+        term_lo +=  2.0 * (ell - 1.0) * (x_lo * dPl_lo - P_lo[ell])
+                 -  2.0 * (ell + 2.0) * dPlm1_lo;
+      }
+      else {
+        term_hi += -2.0 * (ell - 1.0) * (x_hi * dPl_hi - P_hi[ell])
+                   +2.0 * (ell + 2.0) * dPlm1_hi;
+        term_lo += -2.0 * (ell - 1.0) * (x_lo * dPl_lo - P_lo[ell])
+                   +2.0 * (ell + 2.0) * dPlm1_lo;
+      }
+
+      avgGpm = (term_hi - term_lo) / dx;
+
+      pref = 2.0 * (2.0 * ell + 1.0) /
+             ((double)ell * ell * (ell + 1.0) * (ell + 1.0));
+
+      Pl_theta[ell] = pref * avgGpm;
+    }
+  }
+
+  free(P_hi);
+  free(P_lo);
+
 }
 
 /*--------ROUTINE: ccl_tracer_corr_legendre ------
@@ -306,13 +503,6 @@ static void ccl_tracer_corr_legendre(ccl_cosmology *cosmo,
   double *l_arr = NULL, *cl_arr = NULL, *Pl_theta = NULL;
   ccl_f1d_t *cl_spl;
 
-  if(corr_type==CCL_CORR_LM || corr_type==CCL_CORR_LP){
-    *status=CCL_ERROR_NOT_IMPLEMENTED;
-    ccl_cosmology_set_status_message(cosmo,
-                                     "ccl_correlation.c: ccl_tracer_corr_legendre(): "
-                                     "CCL does not support full-sky xi+- calcuations.\nhttps://arxiv.org/abs/1702.05301 indicates flat-sky to be sufficient.\n");
-  }
-
   if(*status==0) {
     l_arr=malloc(((int)(cosmo->spline_params.ELL_MAX_CORR)+1)*sizeof(double));
     if(l_arr==NULL) {
@@ -395,6 +585,103 @@ static void ccl_tracer_corr_legendre(ccl_cosmology *cosmo,
   free(cl_arr);
 }
 
+
+/*--------ROUTINE: ccl_tracer_corr_legendre_binned ------
+TASK: Compute correlation function via Legendre polynomials, including angular binning
+INPUT: cosmology, number of theta bins, theta array, tracer 1, tracer 2, i_bessel, boolean
+       for tapering, vector of tapering limits, correlation vector, angular_cl function.
+ */
+static void ccl_tracer_corr_legendre_binned(ccl_cosmology *cosmo,
+                                     int n_ell,double *ell,double *cls,
+                                     int n_theta,double *theta_min,double *theta_max,double *wtheta,
+                                     int corr_type,int do_taper_cl,double *taper_cl_limits,
+                                     int *status) {
+  int i;
+  double *l_arr = NULL, *cl_arr = NULL, *Pl_theta = NULL;
+  ccl_f1d_t *cl_spl;
+
+  if(*status==0) {
+    l_arr=malloc(((int)(cosmo->spline_params.ELL_MAX_CORR)+1)*sizeof(double));
+    if(l_arr==NULL) {
+      *status=CCL_ERROR_MEMORY;
+      ccl_cosmology_set_status_message(cosmo,
+                                       "ccl_correlation.c: ccl_tracer_corr_legendre(): "
+                                       "ran out of memory\n");
+    }
+  }
+
+  if(*status==0) {
+    cl_arr=malloc(((int)(cosmo->spline_params.ELL_MAX_CORR)+1)*sizeof(double));
+    if(cl_arr==NULL) {
+      *status=CCL_ERROR_MEMORY;
+      ccl_cosmology_set_status_message(cosmo,
+                                       "ccl_correlation.c: ccl_tracer_corr_legendre(): "
+                                       "ran out of memory\n");
+    }
+  }
+
+  if(*status==0) {
+    //Interpolate input Cl into
+    cl_spl=ccl_f1d_t_new(n_ell,ell,cls,cls[0],0,
+			 ccl_f1d_extrap_const,
+			 ccl_f1d_extrap_logx_logy, status);
+    if(cl_spl==NULL) {
+      *status=CCL_ERROR_MEMORY;
+      ccl_cosmology_set_status_message(cosmo,
+                                       "ccl_correlation.c: ccl_tracer_corr_legendre(): "
+                                       "ran out of memory\n");
+    }
+  }
+
+  if(*status==0) {
+    for(i=0;i<=(int)(cosmo->spline_params.ELL_MAX_CORR);i++) {
+      double l=(double)i;
+      l_arr[i]=l;
+      cl_arr[i]=ccl_f1d_t_eval(cl_spl,l);
+    }
+    ccl_f1d_t_free(cl_spl);
+
+    if (do_taper_cl)
+      *status=taper_cl((int)(cosmo->spline_params.ELL_MAX_CORR)+1,l_arr,cl_arr,taper_cl_limits);
+  }
+
+  int local_status, i_L;
+#pragma omp parallel default(none) \
+                     shared(cosmo, theta_min, theta_max, cl_arr, wtheta, n_theta, status, corr_type) \
+                     private(Pl_theta, i, i_L, local_status)
+  {
+    Pl_theta = NULL;
+    local_status = *status;
+
+    if (local_status == 0) {
+      Pl_theta = malloc(sizeof(double)*((int)(cosmo->spline_params.ELL_MAX_CORR)+1));
+      if (Pl_theta == NULL) {
+        local_status = CCL_ERROR_MEMORY;
+      }
+    }
+
+    #pragma omp for schedule(dynamic)
+    for (int i=0; i < n_theta; i++) {
+      if (local_status == 0) {
+        wtheta[i] = 0;
+        ccl_compute_legendre_polynomial_binned(corr_type, theta_min[i], theta_max[i], (int)(cosmo->spline_params.ELL_MAX_CORR), Pl_theta);
+        for (i_L=1; i_L < (int)(cosmo->spline_params.ELL_MAX_CORR); i_L+=1)
+          wtheta[i] += cl_arr[i_L]*Pl_theta[i_L];
+        wtheta[i] /= (M_PI*4);
+      }
+    }
+
+    if (local_status) {
+      #pragma omp atomic write
+      *status = local_status;
+    }
+
+    free(Pl_theta);
+  }
+  free(l_arr);
+  free(cl_arr);
+}
+
 /*--------ROUTINE: ccl_tracer_corr ------
 TASK: For a given tracer, get the correlation function. Do so by running
       ccl_angular_cls. If you already have Cls calculated, go to the next
@@ -427,6 +714,43 @@ void ccl_correlation(ccl_cosmology *cosmo,
 
 }
 
+
+/*--------ROUTINE: ccl_correlation_binned ------
+TASK: For a given tracer, get the correlation function. Include angular
+      binning directly in the kernels
+INPUT: cosmology, number of theta values to evaluate = NL, theta vector,
+       tracer 1, tracer 2, i_bessel, key for tapering, limits of tapering
+       correlation function.
+ */
+void ccl_correlation_binned(ccl_cosmology *cosmo,
+                     int n_ell,double *ell,double *cls,
+                     int n_theta,double *theta_min,double *theta_max,double *wtheta,
+                     int corr_type,int do_taper_cl,double *taper_cl_limits,int flag_method,
+                     int *status) {
+  switch(flag_method) {
+  case CCL_CORR_FFTLOG :
+    *status=CCL_ERROR_INCONSISTENT;
+    ccl_cosmology_set_status_message(cosmo, 
+                                     "ccl_correlation.c: ccl_correlation_binned(): "
+                                     "Cannot use method=FFTLOG for binned correlations, use LEGENDRE\n");
+    break;
+  case CCL_CORR_LGNDRE :
+    ccl_tracer_corr_legendre_binned(cosmo,n_ell,ell,cls,n_theta,theta_min,theta_max,wtheta,
+                                    corr_type, do_taper_cl,taper_cl_limits,status);
+    break;
+  case CCL_CORR_BESSEL :
+    *status=CCL_ERROR_INCONSISTENT;
+    ccl_cosmology_set_status_message(cosmo, 
+                                     "ccl_correlation.c: ccl_correlation_binned(): "
+                                     "Cannot use method=BESSEL for binned correlations, use LEGENDRE\n");
+    break;
+  default :
+    *status=CCL_ERROR_INCONSISTENT;
+    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_correlation_binned(): Unknown algorithm\n");
+  }
+
+}
+
 /*--------ROUTINE: ccl_correlation_3d ------
 TASK: Calculate the 3d-correlation function. Do so by using FFTLog.