Density profiles completely added to the code now.

cluster_toolkit/__init__.py

@@ -20,3 +20,4 @@
 from .averaging import *
 from .massfunction import *
 from .boostfactors import *
+from .density import *

cluster_toolkit/density.py

@@ -0,0 +1,50 @@
+import cluster_toolkit
+from ctypes import c_double, c_int, POINTER
+import numpy as np
+
+def _dcast(x):
+    return cluster_toolkit._ffi.cast('double*', x.ctypes.data)
+
+def rho_nfw_at_R(R, M, c, om, delta=200):
+    """NFW halo density profile.
+
+    Args:
+    R (float or array like): 3d distances from halo center in Mpc/h comoving
+    M (float): Mass in Msun/h
+    c (float): Concentration
+    om (float): Omega_matter, matter fraction of the density
+    delta (int; optional): Overdensity, default is 200
+
+    Returns:
+    rho_nfw (float or array like): NFW halo density profile in Msun h^2/Mpc^3 comoving.
+
+    """
+    if type(R) is list or type(R) is np.ndarray:
+        rho = np.zeros_like(R)
+        cluster_toolkit._lib.calc_rho_nfw(_dcast(R), len(R), M, c, delta, om, _dcast(rho))
+        return rho
+    else:
+        return cluster_toolkit._lib.rho_nfw_at_R(R, M, c, delta, om)
+
+def rho_einasto_at_R(R, M, rs, alpha, om, delta=200, rhos=-1.):
+    """Einasto halo density profile.
+
+    Args:
+    R (float or array like): 3d distances from halo center in Mpc/h comoving
+    M (float): Mass in Msun/h; not used if rhos is specified
+    rhos (float): Scale density in Msun h^2/Mpc^3 comoving; optional
+    rs (float): Scale radius
+    alpha (float): Profile exponent
+    om (float): Omega_matter, matter fraction of the density
+    delta (int): Overdensity, default is 200
+
+    Returns:
+    rho_einasto (float or array like): Einasto halo density profile in Msun h^2/Mpc^3 comoving.
+
+    """
+    if type(R) is list or type(R) is np.ndarray:
+        rho = np.zeros_like(R)
+        cluster_toolkit._lib.calc_rho_einasto(_dcast(R), len(R), M, rhos, rs, alpha, delta, om, _dcast(rho))
+        return rho
+    else:
+        return cluster_toolkit._lib.rho_einasto_at_R(R, M, rhos, rs, alpha, delta, om)

cluster_toolkit/xi.py

@@ -6,7 +6,7 @@ def _dcast(x):
     return cluster_toolkit._ffi.cast('double*', x.ctypes.data)
 
 def xi_nfw_at_R(R, M, c, om, delta=200):
-    """NFW halo profile.
+    """NFW halo profile correlation function.
 
     Args:
     R (float or array like): 3d distances from halo center in Mpc/h comoving
@@ -32,7 +32,7 @@ def xi_einasto_at_R(R, M, rs, alpha, om, delta=200, rhos=-1.):
     Args:
     R (float or array like): 3d distances from halo center in Mpc/h comoving
     M (float): Mass in Msun/h; not used if rhos is specified
-    rhos (float): Scale density in Msun h/pc^2 comoving; optional
+    rhos (float): Scale density in Msun h^2/Mpc^3 comoving; optional
     rs (float): Scale radius
     alpha (float): Profile exponent
     om (float): Omega_matter, matter fraction of the density

include/C_density.h

@@ -0,0 +1,4 @@
+double rho_nfw_at_R(double, double, double, int, double);
+double rho_einasto_at_R(double, double, double, double, double, int, double);
+int calc_rho_nfw(double*, int, double, double, int, double, double*);
+int calc_rho_einasto(double*, int, double, double, double, double, int, double, double*);

src/C_density.c

@@ -0,0 +1,62 @@
+#include "C_xi.h"
+
+#include "gsl/gsl_spline.h"
+#include "gsl/gsl_sf_gamma.h"
+#include <math.h>
+#include <stdio.h>
+
+#define rhomconst 2.775808e+11
+//1e4*3.*Mpcperkm*Mpcperkm/(8.*PI*G); units are Msun h^2/Mpc^3
+
+double rho_nfw_at_R(double R, double Mass, double conc, int delta, double om){
+  double rhom = om*rhomconst;//Msun h^2/Mpc^3
+  double Rdelta = pow(Mass/(1.33333333333*M_PI*rhom*delta), 0.33333333333);
+  double Rscale = Rdelta/conc;
+  double fc = log(1.+conc)-conc/(1.+conc); 
+  return Mass/(4.*M_PI*Rscale*Rscale*Rscale*fc)/(R/Rscale*(1+R/Rscale)*(1+R/Rscale));
+}
+
+int calc_rho_nfw(double*R, int NR, double Mass, double conc, int delta, double om, double*rho_nfw){
+  int i;
+  for(i = 0; i < NR; i++)
+    rho_nfw[i] = rho_nfw_at_R(R[i], Mass, conc, delta, om);
+  return 0;
+}
+
+double rho_einasto_at_R(double R, double Mass, double rhos, double rs, double alpha, int delta, double om){
+  double x;
+  double rhom = om*rhomconst;//Msun h^2/Mpc^3;
+  if (rhos < 0){
+    //Calculate rhos from Mass
+    // Rdelta in Mpc/h comoving
+    double Rdelta = pow(Mass/(1.3333333333333*M_PI*rhom*delta), 0.333333333333);
+    x = 2./alpha * pow(Rdelta/rs, alpha); // Mpc/h comoving
+    double a = 3./alpha;
+    double gam = gsl_sf_gamma(a) - gsl_sf_gamma_inc(a, x);
+    double num = delta * rhom * Rdelta*Rdelta*Rdelta * alpha * pow(2./alpha, a);
+    double den = 3. * rs*rs*rs * gam;
+    rhos = num/den;
+  }
+  x = 2./alpha * pow(R/rs, alpha);
+  return rhos * exp(-x);
+}
+
+int calc_rho_einasto(double*R, int NR, double Mass, double rhos, double rs, double alpha, int delta, double om, double*rho_einasto){
+  if (rhos < 0){
+    //Calculate rhos from Mass
+    double rhom = om*rhomconst;//Msun h^2/Mpc^3
+    // Rdelta in Mpc/h comoving
+    double Rdelta = pow(Mass/(1.3333333333333*M_PI*rhom*delta), 0.333333333333);
+    double x = 2./alpha * pow(Rdelta/rs, alpha); 
+    double a = 3./alpha;
+    double gam = gsl_sf_gamma(a) - gsl_sf_gamma_inc(a, x);
+    double num = delta * rhom * Rdelta*Rdelta*Rdelta * alpha * pow(2./alpha, a);
+    double den = 3. * rs*rs*rs * gam;
+    rhos = num/den;
+  }
+  int i;
+  for(i = 0; i < NR; i++)
+    rho_einasto[i] = rho_einasto_at_R(R[i], Mass, rhos, rs, alpha, delta, om);
+
+  return 0;
+}