<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 ClusterEnsemble demonstration
Some example usage of how to build up a dataframe of galaxy cluster properties, including NFW halo profiles. Each cluster is treated as an individual, meaning we track its individual mass and redshift, and other properties. This is useful for fitting a stacked weak lensing profile, for example, where you want to avoid fitting a single average cluster mass. See the project documentation for more details.
Step1:  Create a ClusterEnsemble object by passing in a numpy array (or list) of redshifts
Step2:  Display what we have so far
Below the DataFrame (which so far only contains the cluster redshifts), we see the default assumptions for the power-law slope and normalization that will be used to convert richness $N_{200}$ to mass $M_{200}$. We'll see how to change those parameters below.
Step3:  Add richness values to the dataframe
This step will also generate $M_{200}$, $r_{200}$, $c_{200}$, scale radius $r_s$, and other parameters, assuming the scaling relation given below.
Step4:  Access any column of the dataframe as an array
Notice that astropy units are present for the appropriate columns.
Step5:  If you don't want units, you can get just the values
Step6:  Or access the Pandas DataFrame directly
Step7:  Change the redshifts
These changes will propogate to all redshift-dependant cluster attributes.
Step8:  Change the mass or richness values
Changing mass will affect richness and vice-versa, through the mass-richness scaling relation. These changes will propogate to all mass-dependant cluster attributes, as appropriate.
Step9:  Change the parameters in the mass-richness relation
The mass-richness slope and normalization can both be changed. The new parameters will be applied to the current n200, and will propagate to mass and other dependant quantities.
Step10:  Calculate $\Sigma(r)$ and $\Delta\Sigma(r)$ for NFW model
First select the radial bins in units of Mpc.
Step11:  Calculate Miscentered NFW Profiles
When the true underlying dark matter distribution is offset from the assumed cluster "center" (such as a BCG or some other center proxy), the weak lensing profiles measured around the assumed centers will be different than for the perfectly centered case. One way to account for this is to describe the cluster centroid offsets as a Gaussian distribution around the true centers. We say the probability of an offset is given by
$P(R_\mathrm{off}) = \frac{R_\mathrm{off}}{\sigma_\mathrm{off}^2}e^{-\frac{1}{2}\left(\frac{R_\mathrm{off}}{\sigma_\mathrm{off}}\right)^2}$,
which is parameterized by the width of the 2D offset distribution $\sigma_\mathrm{off} = \sqrt{\sigma_x^2 + \sigma_y^2}$. Then the measured surface mass density is given by 
$\Sigma^\mathrm{sm}(R) = \int_0^\infty \Sigma(R | R_\mathrm{off})\  P(R_\mathrm{off})\ \mathrm{d}R_\mathrm{off}$,
where
$\Sigma(R | R_\mathrm{off}) = \frac{1}{2\pi} \int_0^{2\pi} \Sigma(r')\ \mathrm{d}\theta$, 
and
$r' = \sqrt{R^2 + R_\mathrm{off}^2 - 2 R R_\mathrm{off} \cos{\theta}}$.
More details on the cluster miscentering problem can be found in Ford et al 2015, George et al 2012, and Johnston et al 2007.
To calculate the miscentered profiles, simply create an array of offsets in units of Mpc, and pass it to the calc_nfw method
<ASSISTANT_TASK:>
Python Code:
from __future__ import print_function
import numpy as np
from astropy import units
from matplotlib import pyplot as plt
%matplotlib inline
from clusterlensing import ClusterEnsemble
Explanation: ClusterEnsemble demonstration
Some example usage of how to build up a dataframe of galaxy cluster properties, including NFW halo profiles. Each cluster is treated as an individual, meaning we track its individual mass and redshift, and other properties. This is useful for fitting a stacked weak lensing profile, for example, where you want to avoid fitting a single average cluster mass. See the project documentation for more details.
End of explanation
z = [0.1,0.2,0.3]
c = ClusterEnsemble(z)
c.describe
Explanation: Create a ClusterEnsemble object by passing in a numpy array (or list) of redshifts
End of explanation
c.show()
Explanation: Display what we have so far
Below the DataFrame (which so far only contains the cluster redshifts), we see the default assumptions for the power-law slope and normalization that will be used to convert richness $N_{200}$ to mass $M_{200}$. We'll see how to change those parameters below.
End of explanation
n200 = np.ones(3)*20.
c.n200 = n200
c.show()
Explanation: Add richness values to the dataframe
This step will also generate $M_{200}$, $r_{200}$, $c_{200}$, scale radius $r_s$, and other parameters, assuming the scaling relation given below.
End of explanation
print('z: \t', c.z)
print('n200: \t', c.n200)
print('r200: \t', c.r200)
print('m200: \t', c.m200)
print('c200: \t', c.c200)
print('rs: \t', c.rs)
Explanation: Access any column of the dataframe as an array
Notice that astropy units are present for the appropriate columns.
End of explanation
c.r200.value
Explanation: If you don't want units, you can get just the values
End of explanation
c.dataframe
Explanation: Or access the Pandas DataFrame directly
End of explanation
c.z = np.array([0.4,0.5,0.6])
c.dataframe
Explanation: Change the redshifts
These changes will propogate to all redshift-dependant cluster attributes.
End of explanation
c.m200 = [3e13,2e14,1e15]
c.show()
c.n200 = [20,30,40]
c.show()
Explanation: Change the mass or richness values
Changing mass will affect richness and vice-versa, through the mass-richness scaling relation. These changes will propogate to all mass-dependant cluster attributes, as appropriate.
End of explanation
c.massrich_slope = 1.5
c.show()
# Show basic table without Pandas formatting
c.show(notebook = False)
Explanation: Change the parameters in the mass-richness relation
The mass-richness slope and normalization can both be changed. The new parameters will be applied to the current n200, and will propagate to mass and other dependant quantities.
End of explanation
rmin, rmax = 0.1, 5. # Mpc
nbins = 20
rbins = np.logspace(np.log10(rmin), np.log10(rmax), num = nbins)
print('rbins range from', rbins.min(), 'to', rbins.max(), 'Mpc')
c.calc_nfw(rbins)    # calculate the profiles
sigma = c.sigma_nfw  # access the profiles
deltasigma = c.deltasigma_nfw
sigma[0]
fig = plt.figure(figsize=(12,5))
fig.suptitle('Centered NFW Cluster Profiles', size=30)
first = fig.add_subplot(1,2,1)
second = fig.add_subplot(1,2,2)
for rich, profile in zip(c.n200,deltasigma):
    first.plot(rbins, profile, label='$N_{200}=$ '+str(rich))
first.set_xscale('log')
first.set_xlabel('$r\ [\mathrm{Mpc}]$', fontsize=20)
first.set_ylabel('$\Delta\Sigma(r)\ [\mathrm{M}_\mathrm{sun}/\mathrm{pc}^2]$', 
                 fontsize=20)
first.set_xlim(rbins.min(), rbins.max())
first.legend(fontsize=20)
for rich, profile in zip(c.n200,sigma):
    second.plot(rbins, profile, label='$N_{200}=$ '+str(rich))
second.set_xscale('log')
second.set_xlabel('$r\ [\mathrm{Mpc}]$', fontsize=20)
second.set_ylabel('$\Sigma(r)\ [\mathrm{M}_\mathrm{sun}/\mathrm{pc}^2]$', 
                 fontsize=20)
second.set_xlim(0.05, 1.)
second.set_xlim(rbins.min(), rbins.max())
second.legend(fontsize=20)
fig.tight_layout()
plt.subplots_adjust(top=0.88)
Explanation: Calculate $\Sigma(r)$ and $\Delta\Sigma(r)$ for NFW model
First select the radial bins in units of Mpc.
End of explanation
offsets = np.array([0.1,0.1,0.1]) #same length as number of clusters
c.calc_nfw(rbins, offsets=offsets)
deltasigma_off = c.deltasigma_nfw
sigma_off = c.sigma_nfw
fig = plt.figure(figsize=(12,5))
fig.suptitle('Miscentered NFW Cluster Profiles', size=30)
first = fig.add_subplot(1,2,1)
second = fig.add_subplot(1,2,2)
for rich, profile in zip(c.n200,deltasigma_off):
    first.plot(rbins, profile, label='$N_{200}=$ '+str(rich))
first.set_xscale('log')
first.set_xlabel('$r\ [\mathrm{Mpc}]$', fontsize=20)
ytitle = '$\Delta\Sigma^\mathrm{sm}(r)\ [\mathrm{M}_\mathrm{sun}/\mathrm{pc}^2]$'
first.set_ylabel(ytitle, fontsize=20)
first.set_xlim(rbins.min(), rbins.max())
first.legend(fontsize=20)
for rich, profile in zip(c.n200,sigma_off):
    second.plot(rbins, profile, label='$N_{200}=$ '+str(rich))
second.set_xscale('log')
second.set_xlabel('$r\ [\mathrm{Mpc}]$', fontsize=20)
ytitle = '$\Sigma^\mathrm{sm}(r)\ [\mathrm{M}_\mathrm{sun}/\mathrm{pc}^2]$'
second.set_ylabel(ytitle, fontsize=20)
second.set_xlim(rbins.min(), rbins.max())
second.legend(fontsize=20)
fig.tight_layout()
plt.subplots_adjust(top=0.88)
Explanation: Calculate Miscentered NFW Profiles
When the true underlying dark matter distribution is offset from the assumed cluster "center" (such as a BCG or some other center proxy), the weak lensing profiles measured around the assumed centers will be different than for the perfectly centered case. One way to account for this is to describe the cluster centroid offsets as a Gaussian distribution around the true centers. We say the probability of an offset is given by
$P(R_\mathrm{off}) = \frac{R_\mathrm{off}}{\sigma_\mathrm{off}^2}e^{-\frac{1}{2}\left(\frac{R_\mathrm{off}}{\sigma_\mathrm{off}}\right)^2}$,
which is parameterized by the width of the 2D offset distribution $\sigma_\mathrm{off} = \sqrt{\sigma_x^2 + \sigma_y^2}$. Then the measured surface mass density is given by 
$\Sigma^\mathrm{sm}(R) = \int_0^\infty \Sigma(R | R_\mathrm{off})\  P(R_\mathrm{off})\ \mathrm{d}R_\mathrm{off}$,
where
$\Sigma(R | R_\mathrm{off}) = \frac{1}{2\pi} \int_0^{2\pi} \Sigma(r')\ \mathrm{d}\theta$, 
and
$r' = \sqrt{R^2 + R_\mathrm{off}^2 - 2 R R_\mathrm{off} \cos{\theta}}$.
More details on the cluster miscentering problem can be found in Ford et al 2015, George et al 2012, and Johnston et al 2007.
To calculate the miscentered profiles, simply create an array of offsets in units of Mpc, and pass it to the calc_nfw method:
End of explanation
<END_TASK>