Implemented the Y(R) profiles.

cluster_toolkit/__init__.py

@@ -34,4 +34,4 @@ def _dcast(x):
     if isinstance(x, list): x = np.asarray(x, dtype=np.float64, order='C')
     return _ffi.cast('double*', x.ctypes.data)
 
-from . import averaging, bias, boostfactors, concentration, deltasigma, density, exclusion, massfunction, miscentering, peak_height, profile_derivatives, xi
+from . import averaging, bias, boostfactors, concentration, deltasigma, density, exclusion, massfunction, miscentering, peak_height, profile_derivatives, sigma_reconstruction, xi

cluster_toolkit/sigma_reconstruction.py

@@ -0,0 +1,44 @@
+"""Galaxy cluster Sigma(R) reconstruction profiles, also known as 'Y' profiles.
+"""
+
+import cluster_toolkit
+from cluster_toolkit import _dcast
+import numpy as np
+
+def Sigma_REC_from_DeltaSigma(R, DeltaSigma):
+    """Reconstructed Sigma(R) profile, also known as 'Y'
+    in the same units as DeltaSigma.
+
+    Note: R and DeltaSigma must have the same shape,
+    and have more than 1 element. The returned array
+    has one fewer elements.
+
+    Note: R must have constant logarithmic spacing.
+
+    Args:
+        R (array like): Projected radii.
+        DeltaSigma (array like): Differential surface mass density.
+
+    Returns:
+        Reconstructed surface mass density.
+    """
+    R = np.asarray(R)
+    DeltaSigma = np.asarray(DeltaSigma)
+    if R.shape != DeltaSigma.shape:
+        raise Exception("R and DeltaSigma must have the same shape.")
+
+    lnR = np.log(R)
+    for i in range(len(lnR)-2):
+        dlnR0 = lnR[i] - lnR[i+1]
+        dlnR1 = lnR[i+1] - lnR[i+2]
+        if not (np.fabs(dlnR0 - dlnR1) < 1e-6):
+            raise Exception("R must have constant logarithmic spacing.")
+        continue
+
+    #Note the order here, we integrate DOWNWARD
+    dlnR = lnR[0] - lnR[1]
+
+    Sigma = np.zeros(len(R)-1)
+    cluster_toolkit._lib.Sigma_REC_from_DeltaSigma(dlnR, _dcast(DeltaSigma),
+                                                   len(R), _dcast(Sigma))
+    return Sigma

include/C_sigma_reconstruction.h

@@ -0,0 +1,2 @@
+int Sigma_REC_from_DeltaSigma(double dlnR, double*DeltaSigma, int N,
+			      double*Sigma);

src/C_sigma_reconstruction.c

@@ -0,0 +1,54 @@
+/** @file C_sigma_reconstruction.c
+ *  @brief Reconstructed Sigma(R) profiles from large scales.
+ * 
+ *  These experimental routines compute the "Y" 
+ *  cluster profiles. This is also known as the
+ *  Sigma(R) reconstruction profiles.
+ *
+ *  This file will undergo much cleanup in the near future.
+ *  
+ *  @author Tom McClintock (tmcclintock)
+ *  @bug No known bugs.
+ */
+
+#include "C_sigma_reconstruction.h"
+
+#include <stdio.h>
+
+/**
+ * \brief Reconstructed Sigma profiles (i.e. Y)
+ * from a precomputed DeltaSigma profile.
+ * 
+ * Note: output Sigma profile has one less element than 
+ * the input DeltaSigma, also we assume a regular grid
+ * spacing in the ln(R).
+ */
+int Sigma_REC_from_DeltaSigma(double dlnR, double*DeltaSigma, int N,
+			      double*Sigma){
+  /*
+    The transformation is:
+
+    Y = T*DeltaSigma = (2S - I)*DeltaSigma
+
+    Where I is the identity matrix and S contains the Runge-kutta
+    elements for a midpoint integration method (i.e.
+    1/2s on the edges and 1s in the middle)
+   */
+  //Loop over Sigma elements; i.e. rows
+  for(int i = 0; i < N-1; i++){
+    Sigma[i] = -DeltaSigma[i]; //Not sure if this is positive or negative
+
+    //Loop over DeltaSigma elements; i.e. columns
+    //for(int j = N-2-i; j < N; j++){
+    for(int j = i; j < N; j++){
+      Sigma[i] -= 2*dlnR*DeltaSigma[j];
+      //Downweight the first and last contributions (midpoint formula)
+      //if ((j == N-2-i) || (j == N-1)){
+      if ((j == i) || (j == N-1)){
+	Sigma[i] += dlnR*DeltaSigma[j];
+      }
+      //printf("%e\n",Sigma[i]);
+    }
+  }
+  return 0;
+}

tests/test_sigma_reconstruction.py

@@ -0,0 +1,5 @@
+import pytest
+from cluster_toolkit import sigma_reconstruction as SR
+from os.path import dirname, join
+import numpy as np
+import numpy.testing as npt