# Halo profiles
This example illustrates how to use CCL to calculate halo profiles (and quantities associated to them)


```python
import numpy as np
import pyccl as ccl
import pylab as plt
import matplotlib.cm as cm
%matplotlib inline
```

## Preliminaries
Generate a cosmology object and a couple of $r$ (comoving radius), $k$ (comoving wavenumber) and $M$ (halo mass) arrays.


```python
# Cosmology
cosmo = ccl.Cosmology(Omega_c=0.27, Omega_b=0.045, h=0.67, sigma8=0.83, n_s=0.96)

# Radius, wavenumber and mass
r_arr = np.geomspace(1E-2,1E2,256)
k_arr = np.geomspace(1E-2,1E2,256)
m_arr = np.geomspace(1E11,1E15,32)

# We will calculate everything at z=1 (scale factor = 0.5)
a_sf = 0.5

# We will also use a mass definition with Delta = 200, and the Duffy 2008 concentration-mass relation
hmd_200c = ccl.halos.MassDef200c()
cM = ccl.halos.ConcentrationDuffy08(hmd_200c)
```

## Halo profile classes
All halo profiles are sub-classes of the base `HaloProfile` class. CCL currently supports the following halo profiles:
- `HaloProfileNFW`: the Navarro-Frenk-White density profile.
- `HaloProfileEinasto`: the Einasto profile.
- `HaloProfileHernquist`: the Hernquist profile.
- `HaloProfileGaussian`: a toy Gaussian profile.
- `HaloProfilePowerLaw`: a toy power-law profile.

However, the API is flexible enough for you to implement your own halo profiles. See further down.

For now, let's compare the first three profiles.


```python
profs = {}
profs['NFW'] = ccl.halos.HaloProfileNFW(cM, truncated=False)
profs['Einasto'] = ccl.halos.HaloProfileEinasto(cM, truncated=False)
profs['Hernquist'] = ccl.halos.HaloProfileHernquist(cM, truncated=False)
```

First, let's just plot the real-space profile for all models


```python
# Plot the real-space profile for a fixed mass of M = 1E14 M_sun
plt.figure()
for n,p in profs.items():
    rho_r = p.real(cosmo, r_arr, 1E14, a_sf, mass_def=hmd_200c)
    plt.plot(r_arr, rho_r, label=n)
plt.xscale('log')
plt.yscale('log')
plt.legend(loc='lower left', fontsize=15)
plt.ylabel(r'$\rho(r)\,\,[M_\odot\,{\rm Mpc}^{-3}]$', fontsize=15)
plt.xlabel(r'$r\,\,[{\rm Mpc}]$', fontsize=15);
```

All the profile-related methods can take in an array of halo masses, and output the profile for all of them in a 2D array.


```python
pM = profs['NFW'].real(cosmo,r_arr, m_arr, a_sf, mass_def=hmd_200c)
pM.shape

fig=plt.figure()
cmap = plt.cm.inferno
for i, p in enumerate(pM):
    plt.plot(r_arr, p, c=cmap((i + 0.5) / 32))
plt.xscale('log')
plt.yscale('log')
lm = np.log10(m_arr)
sm = plt.cm.ScalarMappable(cmap=cmap,
                           norm=plt.Normalize(vmin=lm.min(),
                                              vmax=lm.max()))
sm._A = []
plt.ylabel(r'$\rho(r)\,\,[M_\odot\,{\rm Mpc}^{-3}]$', fontsize=15)
plt.xlabel(r'$r\,\,[{\rm Mpc}]$', fontsize=15);
cbar = fig.colorbar(sm)
cbar.set_label(r'$\log_{10}(M/M_\odot)$', fontsize=15)
```

## Derived quantities
Besides the 3D real-space profile, CCL can also calculate the Fourier-space profile ($\rho(k)$), the 2D projected profile ($\Sigma(R)$) and the cumulative surface density ($\Sigma(<R)$). These are defined as:

\begin{align}
   \rho(k) &= 4\pi\int_0^\infty dr\,r^2\,\rho(r)\,j_0(kr),\\
   \Sigma(R) &= \int_{-\infty}^\infty dr_\parallel\,\rho\left(\sqrt{r_\parallel^2+R^2}\right),\\
   \Sigma(<R) &= \frac{2}{R^2}\int_0^R\,dR'\,R'\,\Sigma(R).
\end{align}

CCL can compute these quantities automatically based on the 3D real-space profile (although you can also choose to calculate them analytically in some particular cases). To do so, CCL uses a fast `FFTLog`-based algorithm.


### Fourier profile


```python
# Plot the Fourier-space profile for a fixed mass of M = 1E14 M_sun
plt.figure()
for n,p in profs.items():
    rho_k = p.fourier(cosmo, k_arr, 1E14, a_sf, mass_def=hmd_200c)
    plt.plot(k_arr, rho_k, label=n)
plt.xscale('log')
plt.yscale('log')
plt.legend(loc='lower left', fontsize=15)
plt.ylabel(r'$\rho(k)\,\,[M_\odot]$', fontsize=15)
plt.xlabel(r'$k\,\,[{\rm Mpc}^{-1}]$', fontsize=15);
```

### Projected profile


```python
# Plot the projected profile for a fixed mass of M = 1E14 M_sun
plt.figure()
for n,p in profs.items():
    sig_R = p.projected(cosmo, r_arr, 1E14, a_sf, mass_def=hmd_200c)
    plt.plot(r_arr, sig_R, label=n)
plt.xscale('log')
plt.yscale('log')
plt.legend(loc='lower left', fontsize=15)
plt.ylabel(r'$\Sigma(R)\,\,[M_\odot\,{\rm Mpc}^{-2}]$', fontsize=15)
plt.xlabel(r'$R\,\,[{\rm Mpc}]$', fontsize=15);
```

### Cumulative surface density


```python
# Plot the projected profile for a fixed mass of M = 1E14 M_sun
plt.figure()
for n,p in profs.items():
    sig_lR = p.cumul2d(cosmo, r_arr, 1E14, a_sf, mass_def=hmd_200c)
    plt.plot(r_arr, sig_lR, label=n)
plt.xscale('log')
plt.yscale('log')
plt.legend(loc='lower left', fontsize=15)
plt.ylabel(r'$\Sigma(<R)\,\,[M_\odot\,{\rm Mpc}^{-2}]$', fontsize=15)
plt.xlabel(r'$R\,\,[{\rm Mpc}]$', fontsize=15);
```

## Creating new profiles
Creating new profiles is relatively easy. In the simplest case you only need to create a class that inherits from `HaloProfile` and that contains a method called `_real` with the signature below. It should be able to take radii and masses as arrays or floats, and it should return a float, 1D array or 2D array (with shape `(n_radius, n_mass)` in the latter case) depending on the input.

To demonstrate this, let's define an "Lorentzian" profile. I.e.:

\begin{equation}
  \rho(r) = \rho_0\,\frac{r_s}{r}\,e^{-r/r_s},\hspace{12pt}\rho_0=\frac{M}{4\pi\,r_s^3}
\end{equation}


```python
class HaloProfileLorentz(ccl.halos.HaloProfile):
    def __init__(self):
        super(HaloProfileLorentz, self).__init__()

    def _rs(self, cosmo, M, a, mdef):
        """ This defines the scale radius as a function
        of halo mass. For simplicity we just use the
        overdensity radius.
        """
        return mdef.get_radius(cosmo, M, a) / a

    def _real(self, cosmo, r, M, a, mass_def):
        """ Real space profile.
        
        Args:
            cosmo (obj:`Cosmology`): a cosmology object.
            r (float or array_like): halo radius.
            M (float or array_like): halo mass.
            a (float): scale factor.
            mass_def (obj:`MassDef`): mass definition.

       Returns:
            float or array_like: halo profile. The shape of the
            output will be `(N_M, N_r)` where `N_r` and `N_m` are
            the sizes of `r` and `M` respectively.
        """
        # Generate 2D array by default
        r_use = np.atleast_1d(r)
        M_use = np.atleast_1d(M)

        # Compute profile
        r_s = self._rs(cosmo, M_use, a, mass_def)
        rho0 = M_use / (4 * np.pi * r_s**3)
        x = r[None, :] / r_s[:, None]
        prof = rho0[:, None] * np.exp(-x) / x

        # Make sure the output has the right shape
        if np.ndim(r) == 0:
            prof = np.squeeze(prof, axis=-1)
        if np.ndim(M) == 0:
            prof = np.squeeze(prof, axis=0)
        return prof
```

Let's generate a profile and plot it.


```python
pl = HaloProfileLorentz()
rs = pl._rs(cosmo, 1E14, a_sf, hmd_200c)
rho0 = 1E14 / (4 * np.pi * rs**3)

# Plotting
plt.figure()
plt.plot(r_arr, pl.real(cosmo, r_arr, 1E14, a_sf, hmd_200c),
         'k-', lw=3)
# Let's compare with the expected shape
plt.plot(r_arr, rho0 * np.exp(-r_arr / rs) * rs / r_arr,
         'y--')
plt.ylabel(r'$\rho(r)\,\,[M_\odot\,{\rm Mpc}^{-3}]$', fontsize=15)
plt.xlabel(r'$r\,\,[{\rm Mpc}]$', fontsize=15)
plt.xscale('log')
```

### Fourier profile
CCL can compute the Fourier-space profile automatically using FFTLog. Let us compare it with the analytical solution for this profile:

\begin{equation}
  \rho(k)=\frac{M}{1+(k r_s)^2}
\end{equation}


```python
# FFTLog calculation
pk = pl.fourier(cosmo, k_arr, 1E14, a_sf, hmd_200c)
# Analytical prediction
pk_pred = 1E14 / (1 + (k_arr * rs)**2)

# Plot!
plt.plot(k_arr, pk,'k-', lw=3, label='FFTLog')
plt.plot(k_arr, pk_pred,'y--', label='Analytic')
plt.legend(loc='upper right', fontsize=15)
plt.ylabel(r'$\rho(k)\,\,[M_\odot]$', fontsize=15)
plt.xlabel(r'$k\,\,[{\rm Mpc}^{-1}]$', fontsize=15);
plt.xscale('log')
```

FFTLog is obviously not doing a very good job here. We can improve the result by tweaking the FFTLog precision parameters through the `update_precision_fftlog` method. This will almost always be the case when you create a new profile if you want to rely on FFTLog to produce these predictions. Below we show how we can improve the results significantly doing this (see the documentation to understand what the different precision parameters do).


```python
# Update parameters
pl.update_precision_fftlog(padding_hi_fftlog=1E2,
                           padding_lo_fftlog=1E-2,
                           n_per_decade=1000,
                           plaw_fourier=-2.)

# Recompute fourier profile
pk = pl.fourier(cosmo, k_arr, 1E14, a_sf, hmd_200c)

# Replot!
plt.plot(k_arr, pk,'k-', lw=3, label='FFTLog')
plt.plot(k_arr, pk_pred,'y--', label='Analytic')
plt.legend(loc='upper right', fontsize=15)
plt.ylabel(r'$\rho(k)\,\,[M_\odot]$', fontsize=15)
plt.xlabel(r'$k\,\,[{\rm Mpc}^{-1}]$', fontsize=15);
plt.xscale('log')
```

### Projected profile and cumulative surface density
Let's now do the same with $\Sigma(R)$ and $\Sigma(<R)$. The analytical predictions for the Lorentzian profile in this case are:

\begin{align}
  \Sigma(R)&=\frac{M}{2\pi r_s^2} K_0(R/r_s),\\
  \Sigma(<R)&=\frac{M}{2\pi r_s^2}\,2\frac{r_s^2}{R^2}\left(1-\frac{R}{r_s}\,K_1(R/r_s)\right),
\end{align}

where $K_\mu(x)$ are the modified Bessel functions of the second kind.


```python
from scipy.special import kn

# Predictions
x = r_arr / rs
sR_pred = 1E14 / (2 * np.pi * rs**2) * kn(0, x)
slR_pred = 1E14 / (2 * np.pi * rs**2) * (2 / x**2) * (1 - x * kn(1, x))

# FFTlog calculation
sR = pl.projected(cosmo, r_arr, 1E14, a_sf, hmd_200c)
slR = pl.cumul2d(cosmo, r_arr, 1E14, a_sf, hmd_200c)

# Projected profile
plt.figure()
plt.plot(r_arr, sR,'k-', lw=3, label='FFTLog')
plt.plot(r_arr, sR_pred,'y--', label='Analytic')
plt.legend(loc='upper right', fontsize=15)
plt.ylabel(r'$\Sigma(R)\,\,[M_\odot\,{\rm Mpc}^{-2}]$', fontsize=15)
plt.xlabel(r'$R\,\,[{\rm Mpc}]$', fontsize=15);
plt.xscale('log')

# Cumulative surface density
plt.figure()
plt.plot(r_arr, slR,'k-', lw=3, label='FFTLog')
plt.plot(r_arr, slR_pred,'y--', label='Analytic')
plt.legend(loc='upper right', fontsize=15)
plt.ylabel(r'$\Sigma(<R)\,\,[M_\odot\,{\rm Mpc}^{-2}]$', fontsize=15)
plt.xlabel(r'$R\,\,[{\rm Mpc}]$', fontsize=15);
plt.xscale('log')
```

24e679d7208a098fa95c1576b199cca2f3727445
ipynb
Jupyter Notebook
Halo profiles.ipynb
bjornvz/CCLX
dfb0fba4114dea267dec59ebca57870493f43f57
BSD-3-Clause
2019-12-08T11:05:29.000Z
2022-02-26T19:13:52.000Z
Halo profiles.ipynb
bjornvz/CCLX
dfb0fba4114dea267dec59ebca57870493f43f57
BSD-3-Clause
2019-11-20T02:17:01.000Z
2022-03-11T11:40:10.000Z
Halo profiles.ipynb
bjornvz/CCLX
dfb0fba4114dea267dec59ebca57870493f43f57
BSD-3-Clause
2020-02-14T10:57:19.000Z
2022-03-28T19:21:20.000Z
Qwen/Qwen-72B
1. YES
2. YES