<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Setup Cosmology
Step1:  Create Stellar Population
Step2:  Calculate a few things to get going.
Step5:  Define the functions that we'll need
Need to compute the cluster volume...
$M_{vir} = 4/3 \pi r^3_{vir} \rho_c(r<r_{vir}) = 4/3 \pi r^3_{vir} \Delta_c \rho_c$
if we let $\Delta_c = 200$ then 
$M_{200} = 4/3 \pi r^3_{200} 200 \rho_c$ with $\rho_c = \frac{3H(z)^2}{8\pi G}$
or just $M_{200} = V_{200}200\rho_c$. So we'll make a function to calculate $\rho_c$. And we'll make use of the astropy units package to do all the unit analysis for us.
Don't forget that $H(z) = H_0E(z)$ 
We also need to integrate the Schechter luminosity functions..
The Schechter Function
Step6:  Functions to calculate everything at a specific redshift
Step7:  Start Calculating things
Step8:  The figure below is where we let both $\phi^\star$ and $\alpha$ evolve with redshift.
The vertical lines show the faint limit to our integral. In all cases we are integrating to $M=-30$ mag. So in the function above, we are integrating from the left to the right in this figure.
Step9:  The figure below is where we fix both $\phi^\star$ and $\alpha$ at some low redshift value, and do not allow them to evolve with redshift.
Here I've fixed them at $z=0.2$.
In my opinion, this looks a lot more like what we would expect for a LF. It still doesn't fully capture the evolution of the LF that I'd expect. I'd expect the
Step10:  Here is the above figure, with $\phi^\star$ and $\alpha$ fixed to their $z=0.2$ values.
<ASSISTANT_TASK:>
Python Code:
cosmo = LambdaCDM(H0=70, Om0=0.3, Ode0=0.7, Tcmb0=2.725)
Explanation: Setup Cosmology
End of explanation
# check to make sure we have defined the bpz filter path
if not os.getenv('EZGAL_FILTERS'):
    os.environ['EZGAL_FILTERS'] = (f'{os.environ["HOME"]}/Projects/planckClusters/MOSAICpipe/bpz-1.99.3/FILTER/')
model = ezgal.model('bc03_ssp_z_0.02_salp.model')
model = model.make_exponential(1)
model.set_cosmology(Om=cosmo.Om0, Ol=cosmo.Ode0, h=cosmo.h, w=cosmo.w(0))
                                   
model.add_filter('g_MOSAICII.res', name='g')
model.add_filter('r_MOSAICII.res', name='r')
model.add_filter('i_MOSAICII.res', name='i')
model.add_filter('z_MOSAICII.res', name='z')
model.add_filter('K_KittPeak.res', name='K')
# Blanton 2003 Normalization
Mr_star = -20.44 + 5 * np.log10(cosmo.h) # abs mag.
# set the normalization
model.set_normalization('sloan_r', 0.1, Mr_star, vega=False)                                   
Explanation: Create Stellar Population
End of explanation
# desired formation redshift
zf = 6.0
# fetch an array of redshifts out to given formation redshift
zs = model.get_zs(zf)
 
# Calculate some cosmological stuff
DM = cosmo.distmod(zs)
dlum = cosmo.luminosity_distance(zs)
Explanation: Calculate a few things to get going.
End of explanation
def rho_crit(z, cosmo):
    # convert G into better units:
    G = const.G.to(u.km**2 * u.Mpc/(u.M_sun * u.s**2))
    return 3 / (8 * np.pi * G) * cosmo.H0**2 * cosmo.efunc(z)**2 # Mpc^3
def schechterL(luminosity, phiStar, alpha, LStar): 
    Schechter luminosity function. 
    LOverLStar = (luminosity/LStar) 
    return (phiStar/LStar) * LOverLStar**alpha * np.exp(- LOverLStar) 
def schechterM(magnitude, phiStar, alpha, MStar): 
    Schechter luminosity function by magnitudes.
    
    # check to make sure things are all the same size
    if isinstance(phiStar, float) and isinstance(magnitude, np.ndarray):
        phiStar = np.ones_like(magnitude) * phiStar
    if isinstance(alpha , float) and isinstance(magnitude, np.ndarray):
        alpha = np.ones_like(magnitude) * alpha
    if isinstance(MStar, float) and isinstance(magnitude, np.ndarray):
        MStar = np.ones_like(magnitude) * MStar
    
    MStarMinM = 0.4 * (MStar - magnitude)
    return (0.4 * np.log(10) * phiStar * 10.0**(MStarMinM * (alpha + 1.)) * np.exp(-10.**MStarMinM))
Explanation: Define the functions that we'll need
Need to compute the cluster volume...
$M_{vir} = 4/3 \pi r^3_{vir} \rho_c(r<r_{vir}) = 4/3 \pi r^3_{vir} \Delta_c \rho_c$
if we let $\Delta_c = 200$ then 
$M_{200} = 4/3 \pi r^3_{200} 200 \rho_c$ with $\rho_c = \frac{3H(z)^2}{8\pi G}$
or just $M_{200} = V_{200}200\rho_c$. So we'll make a function to calculate $\rho_c$. And we'll make use of the astropy units package to do all the unit analysis for us.
Don't forget that $H(z) = H_0E(z)$ 
We also need to integrate the Schechter luminosity functions..
The Schechter Function:
For Luminosity:
$\Phi(L) = \phi^\star \frac{L}{L_\star}^\alpha e^{-\frac{L}{L_\star}}$
For Magnitudes:
$\Phi(M) = \phi^\star\frac{2}{5}log(10) (10^{\frac{2}{5}(M_\star - M)})^{\alpha+1} e^{-10^{\frac{2}{5}(M_\star - M)}}$
End of explanation
def calc_Maglim(mlim, redshift, filter='i', cosmo=cosmo, model=model):
    return mlim - cosmo.distmod(redshift).value - model.get_kcorrects(zf, filters=filter, zs=redshift)
def calc_phistar(redshift, cosmo=cosmo):
    return 3.6 * cosmo.efunc(redshift)**2
def calc_alpha(redshift):
    return -1.05 * (1 + redshift)**(-2 / 3)
def calc_fr(redshift):
    return 0.8 * (1 + redshift)**(-1 / 2)
Explanation: Functions to calculate everything at a specific redshift
End of explanation
# So now we are going to calculate the volumes as a function of z
#M200 = mass_func(zarr) * u.solMass
M200 = 1e15 * u.solMass
V200 = M200/ (200 * rho_crit(zs, cosmo))
# Calculate the M_star values
Mstar = model.get_absolute_mags(zf, filters='i', zs=zs)
# calculate the abs mag of our limiting magnitude as a function of z
mlim = 23.5
#Mlim = Mstar - 2.5 * np.log10(0.4)
Mlim = mlim - cosmo.distmod(zs).value - model.get_kcorrects(zf, filters='i', zs=zs)
# Here are the Schechter function stuff from Liu et al.
phi_star = 3.6 * cosmo.efunc(zs)**2
alpha = -1.05 * (1 + zs)**(-2/3)
fr = 0.8*(1 + zs)**(-1/2)
LF = []
for phi, a, M_star, M_lim in zip(phi_star, alpha, Mstar, Mlim):
    if M_lim < M_star - 2.5 * np.log10(0.4):
        Mlimit = M_lim
    else:
        Mlimit = M_star - 2.5 * np.log10(0.4)
    y, err = quad(schechterM, -30, Mlimit, args=(phi, a, M_star))
    LF.append(y)
plt.figure()
plt.plot(zs, (LF * V200.value + 1) * fr)
ax = plt.gca()
ax.set_yticks(np.arange(0, 75, 10))
plt.xlim(0.1, 5)
plt.ylim(0, 80)
plt.xlabel('redshift')
plt.ylabel('N (r < r200)')
plt.grid()
calc_Maglim(23.5, 0.2)
magarr = np.arange(-30, -14, 0.2)
redshifts = np.arange(0.2, 2, 0.2) 
# make some colored lines for the plots
cm_subset = np.linspace(0.2, 0.8, len(redshifts))
colors = [cm.bone(x) for x in cm_subset]
Explanation: Start Calculating things
End of explanation
plt.figure(figsize=(8,6))
for z, c in zip(redshifts, colors):
    plt.plot(magarr, schechterM(magarr, calc_phistar(z), calc_alpha(z), model.get_absolute_mags(zf, filters='i', zs=z)), c=c)
    plt.axvline(calc_Maglim(23.5, z), c=c, label=f'{z:.2}')
plt.xlabel('Abs Magnitude')
plt.ylabel('$\Phi(M)$ [Ngal $Mpc^{-3}$]')
plt.semilogy()
plt.ylim(0.01, 20)
plt.legend(loc='upper left')
Explanation: The figure below is where we let both $\phi^\star$ and $\alpha$ evolve with redshift.
The vertical lines show the faint limit to our integral. In all cases we are integrating to $M=-30$ mag. So in the function above, we are integrating from the left to the right in this figure.
End of explanation
plt.figure(figsize=(8,6))
for z, c in zip(redshifts, colors):
    plt.plot(magarr, schechterM(magarr, calc_phistar(0.2), calc_alpha(z), model.get_absolute_mags(zf, filters='i', zs=z)), c=c)
    plt.axvline(calc_Maglim(23.5, z), c=c, label=f'{z:.2}')
plt.xlabel('Abs Magnitude')
plt.ylabel('$\Phi(M)$ [Ngal $Mpc^{-3}$]')
plt.semilogy()
plt.ylim(0.01, 10)
plt.legend(loc='upper left')
LF = []
for phi, a, M_star, M_lim in zip(phi_star, alpha, Mstar, Mlim):
    if M_lim < M_star - 2.5 * np.log10(0.4):
        Mlimit = M_lim
    else:
        Mlimit = M_star - 2.5 * np.log10(0.4)
    y, err = quad(schechterM, -30, Mlimit, args=(7, -0.8, M_star))
    #print(M_star - M_lim, y)
    LF.append(y)
Explanation: The figure below is where we fix both $\phi^\star$ and $\alpha$ at some low redshift value, and do not allow them to evolve with redshift.
Here I've fixed them at $z=0.2$.
In my opinion, this looks a lot more like what we would expect for a LF. It still doesn't fully capture the evolution of the LF that I'd expect. I'd expect the
End of explanation
plt.figure()
plt.plot(zs, (LF * V200.value + 1) * fr)
ax = plt.gca()
#ax.set_yticks(np.arange(0, 75, 10))
plt.xlim(0.1, 2)
#plt.ylim(0, 80)
plt.xlabel('redshift')
plt.ylabel('N (r < r200)')
plt.axhline(3)
plt.grid()
LF * V200.value + 1
Mpiv = 6e14
zpiv = 0.46
fr2 = 0.68 * (M200.value / Mpiv)**-0.1 * ((1 + zs)/ (1 + zpiv))** -0.65
plt.plot(fr, fr2)
plt.plot(fr, fr3)
plt.plot([0,1], [0,1], c='k')
from scipy.special import erf
N200 = LF * V200.value + 1
gz = 0.206 - 0.371 * zs
hz = -3.6 + 25.8 * zs
fr3 = gz * erf(np.log10(N200) - np.log10(hz)) + 0.69
fr3
fr2
calc_phistar(0.2)
phi_star
Explanation: Here is the above figure, with $\phi^\star$ and $\alpha$ fixed to their $z=0.2$ values.
End of explanation
<END_TASK>