port to C, operate on arrays of r

Author: jhod0

cluster_toolkit/pressure_profile.py

@@ -17,11 +17,136 @@
 obtain a more precise answer.
 '''
 
+from cluster_toolkit import _dcast, _lib
 import numpy as np
 from scipy.integrate import quad
 import scipy.special as spec
 
 
+def P_BBPS(r, M, z, omega_b, omega_m,
+           params_P_0=(18.1, 0.154, -0.758),
+           params_x_c=(0.497, -0.00865, 0.731),
+           params_beta=(4.35, 0.0393, 0.415),
+           alpha=1, gamma=-0.3,
+           delta=200):
+    '''
+    The best-fit pressure profile presented in BBPS2.
+
+    Args:
+        r (float or array): Radii from the cluster center, \
+                            in Mpc :math:`h^{-2/3}`. If an array, an array \
+                            is returned, if a scalar, a scalar is returned.
+        M (float): Cluster :math:`M_{\\Delta}`, in Msun.
+        z (float): Cluster redshift.
+        omega_b (float): Baryon fraction.
+        omega_m (float): Matter fraction.
+        params_P_0 (tuple): 3-tuple of :math:`P_0` mass, redshift dependence \
+                parameters A, :math:`\\alpha_m`, :math:`\\alpha_z`, \
+                respectively. See BBPS2 Equation 11.
+        params_x_c (tuple): 3-tuple of :math:`x_c` mass, redshift dependence, \
+                same as `params_P_0`.
+        params_beta (tuple): 3-tuple of :math:`\\beta` mass, redshift \
+                dependence, same as `params_P_0`.
+
+    Returns:
+        float: Pressure at distance `r` from the cluster, in units of \
+                :math:`h^{8/3} Msun s^{-2} Mpc^{-1}`.
+    '''
+    r = np.asarray(r, dtype=np.double)
+
+    scalar_input = False
+    if r.ndim == 0:
+        scalar_input = True
+        # Convert r to 2d
+        r = r[None]
+    if r.ndim > 1:
+        raise Exception('r cannot be a >1D array.')
+
+    P_out = np.zeros_like(r, dtype=np.double)
+
+    # Set parameters
+    P_0 = _A_BBPS(M, z, *params_P_0)
+    x_c = _A_BBPS(M, z, *params_x_c)
+    beta = _A_BBPS(M, z, *params_beta)
+
+    _lib.P_BBPS(_dcast(P_out),
+                _dcast(r), len(r),
+                M, z,
+                omega_b, omega_m,
+                P_0, x_c, beta,
+                float(alpha), gamma,
+                delta)
+
+    if scalar_input:
+        return np.squeeze(P_out)
+    return P_out
+
+
+def projected_P_BBPS(r, M, z, omega_b, omega_m,
+                     params_P=(18.1, 0.154, -0.758),
+                     params_x_c=(0.497, -0.00865, 0.731),
+                     params_beta=(4.35, 0.0393, 0.415),
+                     alpha=1, gamma=-0.3,
+                     delta=200,
+                     limit=1000,
+                     epsabs=1e-15, epsrel=1e-3,
+                     return_errs=False):
+    '''
+    Computes the projected line-of-sight pressure of a cluster at a radius r
+    from the cluster center.
+
+    Args:
+        r (float or array): Radius from the cluster center, in Mpc.
+        M (float): Cluster :math:`M_{\\Delta}`, in Msun.
+        z (float): Cluster redshift.
+        omega_b (float): Baryon fraction.
+        omega_m (float): Matter fraction.
+
+    Returns:
+        float or array: Integrated line-of-sight pressure at distance `r` from \
+                        the cluster, in units of Msun s^{-2} h^{8/3}.
+    '''
+    r = np.asarray(r, dtype=np.double)
+
+    scalar_input = False
+    if r.ndim == 0:
+        scalar_input = True
+        # Convert r to 2d
+        r = r[None]
+    if r.ndim > 1:
+        raise Exception('r cannot be a >1D array.')
+
+    P_out = np.zeros_like(r, dtype=np.double)
+    P_err_out = np.zeros_like(r, dtype=np.double)
+
+    # Set parameters
+    P_0 = _A_BBPS(M, z, *params_P)
+    x_c = _A_BBPS(M, z, *params_x_c)
+    beta = _A_BBPS(M, z, *params_beta)
+
+    rc = _lib.projected_P_BBPS(_dcast(P_out), _dcast(P_err_out),
+                               _dcast(r), len(r),
+                               M, z,
+                               omega_b, omega_m,
+                               P_0, x_c, beta,
+                               alpha, gamma,
+                               delta,
+                               limit,
+                               epsabs, epsrel)
+
+    if rc != 0:
+        msg = 'C_projected_P_BBPS returned error code: {}'.format(rc)
+        raise RuntimeError(msg)
+
+    if scalar_input:
+        if return_errs:
+            return np.squeeze(P_out), np.squeeze(P_err_out)
+        return np.squeeze(P_out)
+    if return_errs:
+        return P_out, P_err_out
+    return P_out
+
+
 def _rho_crit(z, omega_m):
     '''
     The critical density of the universe :math:`\\rho_{crit}`, in units of
@@ -45,7 +170,7 @@ def R_delta(M, z, omega_m, delta=200):
 
     Args:
         M (float): Halo mass :math:`M_{\\Delta}`, in units of Msun.
-        omega_b (float): The baryon fraction :math:`\\Omega_b`.
+        z (float): Redshift to the cluster center.
         omega_m (float): The matter fraction :math:`\\Omega_m`.
         delta (float): The halo overdensity :math:`\\Delta`.
 
@@ -111,84 +236,24 @@ def P_simple_BBPS(x, M, z):
     return P_simple_BBPS_generalized(x, M, z, P_0, x_c, beta)
 
 
-def P_BBPS(r, M, z, omega_b, omega_m,
-           params_P_0=(18.1, 0.154, -0.758),
-           params_x_c=(0.497, -0.00865, 0.731),
-           params_beta=(4.35, 0.0393, 0.415),
-           alpha=1, gamma=-0.3,
-           delta=200):
-    '''
-    The best-fit pressure profile presented in BBPS2.
-
-    Args:
-        r (float): Radius from the cluster center, in Mpc :math:`h^{-2/3}`.
-        M (float): Cluster :math:`M_{\\Delta}`, in Msun.
-        z (float): Cluster redshift.
-        omega_b (float): Baryon fraction.
-        omega_m (float): Matter fraction.
-        params_P_0 (tuple): 3-tuple of :math:`P_0` mass, redshift dependence \
-                parameters A, :math:`\\alpha_m`, :math:`\\alpha_z`, \
-                respectively. See BBPS2 Equation 11.
-        params_x_c (tuple): 3-tuple of :math:`x_c` mass, redshift dependence, \
-                same as `params_P_0`.
-        params_beta (tuple): 3-tuple of :math:`\\beta` mass, redshift \
-                dependence, same as `params_P_0`.
-
-    Returns:
-        float: Pressure at distance `r` from the cluster, in units of \
-                :math:`h^{8/3} Msun s^{-2} Mpc^{-1}`.
-    '''
-    P_0 = _A_BBPS(M, z, *params_P_0)
-    x_c = _A_BBPS(M, z, *params_x_c)
-    beta = _A_BBPS(M, z, *params_beta)
-
-    x = r / R_delta(M, z, omega_m, delta=delta)
-    Pd = P_delta(M, z, omega_b, omega_m, delta=delta)
-    return Pd * P_simple_BBPS_generalized(x, M, z, P_0, x_c, beta,
-                                          alpha=alpha, gamma=gamma, delta=delta)
-
-
-def projected_P_BBPS(r, M, z, omega_b, omega_m,
-                     dist=8, epsrel=1e-3):
-    '''
-    Computes the projected line-of-sight density of a cluster at a radius r
-    from the cluster center.
-
-    Args:
-        r (float): Radius from the cluster center, in Mpc.
-        M (float): Cluster :math:`M_{\\Delta}`, in Msun.
-        z (float): Cluster redshift.
-        omega_b (float): Baryon fraction.
-        omega_m (float): Matter fraction.
-
-    Returns:
-        float: Integrated line-of-sight pressure at distance `r` from the \
-               cluster, in units of Msun s^{-2}.
-    '''
-    R_del = R_delta(M, z, omega_m)
-    return quad(lambda x: P_BBPS(np.sqrt(x*x + r*r), M, z,
-                                 omega_b, omega_m),
-                -dist * R_del, dist * R_del,
-                epsrel=epsrel)[0] / (1 + z)
-
-
 def projected_y_BBPS(r, M, z, omega_b, omega_m,
-                     Xh=0.76, dist=8, epsrel=1e-3):
+                     Xh=0.76, epsrel=1e-3):
     '''
     Projected Compton-y parameter along the line of sight, at a perpendicular
     distance `r` from a cluster of mass `M` at redshift `z`. All arguments have
     the same meaning as `projected_P_BBPS`.
 
     Args:
-        r (float): Radius from the cluster center, in Mpc.
+        r (float or array): Radius from the cluster center, in Mpc.
         M (float): Cluster :math:`M_{\\Delta}`, in Msun.
         z (float): Cluster redshift.
         omega_b (float): Baryon fraction.
         omega_m (float): Matter fraction.
         Xh (float): Primordial hydrogen mass fraction.
 
     Returns:
-        float: Compton y parameter. Unitless.
+        float or array: Compton y parameter. Units are :math:`h^{8/3}`, so \
+                        multiply by :math:`h^{8/3}` to obtain the true value.
     '''
     # The constant is \sigma_T / (m_e * c^2), i.e. the Thompson cross-section
     # divided by the mass-energy of the electron, in units of s^2 Msun^{-1}.
@@ -198,11 +263,11 @@ def projected_y_BBPS(r, M, z, omega_b, omega_m,
     # equation to do so, see BBPS2 p. 3.
     ch = (2 * Xh + 2) / (5 * Xh + 3)
     return ch * cy * projected_P_BBPS(r, M, z, omega_b, omega_m,
-                                      dist=dist, epsrel=epsrel)
+                                      epsrel=epsrel)
 
 
 def Cxl(M, z, omega_b, omega_m, da, ell,
-        dist=8, epsrel=1e-3):
+        epsrel=1e-3):
     '''
     The Fourier-transform of a cluster's Compton y parameter on the sky. Assumes
     the flat-sky approximation.
@@ -221,22 +286,50 @@ def Cxl(M, z, omega_b, omega_m, da, ell,
     '''
     return quad(lambda theta: 2 * np.pi * theta * spec.j0(ell * theta)
                 * projected_y_BBPS(theta * da, M, z, omega_b, omega_m,
-                                   dist=dist, epsrel=epsrel),
+                                   epsrel=epsrel),
                 0, 2 * np.pi,
                 epsrel=epsrel)[0]
 
 
 def smoothed_xi(theta, M, z, omega_b, omega_m, da,
-                dist=8, epsrel=1e-3, maxl=10000):
+                epsrel=1e-3, maxl=10000):
     # Convert from arcmin to radians
     theta = theta * 60 * np.pi / 180
     return quad(lambda ell: 1 / (2 * np.pi) * ell * spec.j0(ell * theta)
                 * Cxl(M, z, omega_b, omega_m, da, ell,
-                      dist=dist, epsrel=epsrel),
+                      epsrel=epsrel),
                 0, maxl,
                 epsrel=epsrel)[0]
 
 
+##################################################
+# The following functions are for testing only!! #
+##################################################
+
+def py_projected_P_BBPS(r, M, z, omega_b, omega_m,
+                        dist=8, epsrel=1e-3):
+    '''
+    Computes the projected line-of-sight density of a cluster at a radius r
+    from the cluster center.
+
+    Args:
+        r (float): Radius from the cluster center, in Mpc.
+        M (float): Cluster :math:`M_{\\Delta}`, in Msun.
+        z (float): Cluster redshift.
+        omega_b (float): Baryon fraction.
+        omega_m (float): Matter fraction.
+
+    Returns:
+        float: Integrated line-of-sight pressure at distance `r` from the \
+               cluster, in units of Msun s^{-2}.
+    '''
+    R_del = R_delta(M, z, omega_m)
+    return quad(lambda x: P_BBPS(np.sqrt(x*x + r*r), M, z,
+                                 omega_b, omega_m),
+                -dist * R_del, dist * R_del,
+                epsrel=epsrel)[0] / (1 + z)
+
+
 def projected_P_BBPS_real(r, M, z, omega_b, omega_m, chis, zs,
                           dist=8, epsrel=1e-3):
     '''

include/C_pressure_profile.h

@@ -0,0 +1,2 @@
+void P_BBPS(double *P_out, const double *r, unsigned Nr, double M_delta, double z, double omega_b, double omega_m, double P_0, double x_c, double beta, double alpha, double gamma, double delta);
+int projected_P_BBPS(double *P_out, double *P_err_out, const double *r, unsigned Nr, double M_delta, double z, double omega_b, double omega_m, double P_0, double x_c, double beta, double alpha, double gamma, double delta, unsigned limit, double epsabs, double epsrel);