Purpose: Conveniently obtaining the linear power spectrum (currently only from) CAMB outputs:

  • Method: "CAMBoutfile"
    • from power spectrum output at desired redshift
    • from transfer function output at desired redshift and a cosmological model with Amplitude specified
    • from transfer function output at desired redshift and a cosmological model with sigma8 specified at the same time
  • Method: "CAMBoutgrowth"
    • from transfer function output at z = 0 and a cosmological model with sigma8 at z =0, obtain power spectrum at z < 0.
    • from power spectrum output at z = 0, and a cosmological model, obtain the power spectrum at z < 0.
    • from the transfer function output at z =0, and a cosmological model , As, obtain power spectrum at z < 0.

Some Options:

  • matter power spectrum: matter = CDM, Baryons, Massive neutrinos
  • cb power spectrum = CDM , baryons only.
  • Normalization: is normalized to sigma8 for the cb part or the total matter fluctuations
  • Obtain the logarithmically or inearly interpolated power spectrum at a list of wave numbers. For wavenumbers supplied outside the range of data points, the extrpolated values are given by np.nan, by default.

In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import camb_utils.cambio as cio
import psutils as psu
from interfaces import FCPL
import utils.plotutils as pu
import matplotlib.gridspec as gridspec
import matplotlib.ticker as ticker
#plt.ticklabel_format(style='sci',axis='x',scilimits=(0,0))
#majorformatter = ticker.ScalarFormatter(useOffset =False)

In [2]:
# INPUT FILES (CAMB outputs)
dirn  = "./"
#M000 (LCDM)
tkfile  = dirn +"example_data/M000/m000n0_transfer_fin_out.dat"
pkfile = dirn + "example_data/M000/m000n0_matterpower_fin.dat"
#M000n1 (same LCDM, with some fraction of CDM replaced by massive neutrinos)
ntkfile = dirn +"example_data/M000n1/m000n1_transfer_fin_out.dat"
npkfile = dirn +"example_data/M000n1/m000n1_matterpower_fin.dat"

In [9]:
#Cosmological Models:
#               Same cosmology represented differently
M000 = FCPL(H0 = 71. ,Om0  = 0.26479, Ob0 = 0.044793, ns = 0.963, As = 2.16e-9)
M000s  = FCPL(H0 = 71., Om0 = 0.26479, Ob0 = 0.044793, ns = 0.963, sigma8 = 0.8)
#
M000n1 = FCPL(H0 = 71., Om0 = 0.26479, Ob0 = 0.044793, ns = 0.963, sigma8 = 0.8, sigmamnu=0.94)

Results for M000 (LCDM)


In [10]:
#powerspectrum 
#use power spectrum from CAMB power spectrum output
psfrompk = psu.powerspectrum(koverh = None, asciifile =pkfile)
#use matter power spectrum from CAMB transfer function output, with As
psfromtkm = psu.powerspectrum(koverh = None, pstype ="matter", asciifile =tkfile, cosmo = M000)
#use matter power spectrum from CAMB transfer function output, with sigma8
psfromtkms = psu.powerspectrum(koverh = None, pstype ="matter", sigma8type = "matter", asciifile =tkfile, cosmo = M000s)
#use cb power spectrum from CAMB transfer function output, normalized to matter
psfromtkcbs = psu.powerspectrum(koverh = None, pstype ="cb", sigma8type = "matter", asciifile =tkfile, cosmo = M000s)
#use cb power spectrum from CAMB transfer function output, normalized to cb
psfromtkcbscb = psu.powerspectrum(koverh = None, pstype ="cb", sigma8type = "cb", asciifile =tkfile, cosmo = M000s)
#
#psfrompkn = psu.powerspectrum(koverh = None, asciifile  = npkfile)
#psfromtknm = psu.powerspectrum(koverh = None, asciifile  = ntkfile, pstype ="matter", cosmo = M000n1)

In [11]:
#Want to check results to 1 percent 
majorformatter = ticker.ScalarFormatter(useOffset =False)
pscompfig, pscomp_ax0, pscomp_ax1 = pu.settwopanel(setdifflimits=  [0.9995,1.0005])
pscomp_ax0.loglog(psfromtkm[0], psfromtkm[1], label = "TKm, As = 2.16e-9")
pscomp_ax0.loglog(psfromtkms[0], psfromtkms[1], label = r'TKm, $\sigma_8 =0.8$')
pscomp_ax0.loglog(psfromtkcbs[0], psfromtkcbs[1], label = "TKcb, $\sigma_8 =0.8$, matter")
pscomp_ax0.loglog(psfromtkcbs[0], psfromtkcbscb[1], label = "TKcb, $\sigma_8 =0.8$, cb")
pscomp_ax0.loglog(psfrompk[0],psfrompk[1],label="PS File")
pscomp_ax0.legend(loc="lower center")
pscomp_ax1.plot(psfromtkm[0],psfromtkm[1]/psfrompk[1])
pscomp_ax1.plot(psfromtkm[0],psfromtkms[1]/psfrompk[1])
pscomp_ax1.plot(psfromtkcbs[0], psfromtkcbs[1]/psfrompk[1], label = "TKcb, $\sigma_8 =0.8$, matter")
pscomp_ax1.plot(psfromtkcbs[0], psfromtkcbscb[1]/psfrompk[1], label = "TKcb, $\sigma_8 =0.8$, cb")
pscomp_ax1.yaxis.set_major_formatter(majorformatter)
pscomp_ax1.set_xscale('log')
#ppscomp_ax1.get_yticklabels()
#pscomp_ax1.set_yscale('log')


Caption: Comparison of the different modes of computing power spectrum for LCDM models. These should be identical, and are very nearly so, as shown by the next plot in greater detail.

Check only the cb components:


In [6]:
cbfig, cb_ax0, cb_ax1 = pu.settwopanel(setdifflimits = [0.997,1.003])
cb_ax0.loglog(psfromtkcbs[0], psfromtkcbs[1], label = "TKcb, $\sigma_8 =0.8$, matter")
cb_ax0.loglog(psfromtkcbs[0], psfromtkcbscb[1], label = "TKcb, $\sigma_8 =0.8$, cb")
cb_ax0.loglog(psfrompk[0],psfrompk[1],label="PS File")
cb_ax0.legend(loc="lower center")
cb_ax0.set_ylabel("PK ")
cb_ax1.plot(psfromtkcbs[0], psfromtkcbs[1]/psfrompk[1], label = "TKcb, $\sigma_8 =0.8$, matter")
cb_ax1.plot(psfromtkcbs[0], psfromtkcbscb[1]/psfrompk[1],ls='dashed', label = "TKcb, $\sigma_8 =0.8$, cb")
cb_ax1.set_xscale('log')
plt.setp(cb_ax1.get_yticklabels()[1::2],visible=False)
cb_ax1.set_xlabel(r'Mass ($h^{-1} M_\odot)$)')
#cb_ax1.set_ylabel("PK / PK (from pk)")


Out[6]:
<matplotlib.text.Text at 0x10fec9450>

In [12]:
#Only the case where As is provided is at 1.0000, all other cases where sigma8 is calculated
#are lower 
plt.plot(psfromtkcbs[0], psfromtkcbscb[1]/psfrompk[1],'k',ls='dashed', label = "TKcb, $\sigma_8 =0.8$, cb")
plt.plot(psfromtkcbs[0], psfromtkms[1]/psfrompk[1],'r',ls='dotted', label = "TKm, $\sigma_8 =0.8$")
plt.plot(psfromtkcbs[0], psfromtkcbs[1]/psfrompk[1],'b-' ,label = "TKcb, $\sigma_8 =0.8$")
plt.plot(psfromtkcbs[0], psfromtkm[1]/psfrompk[1],'g-' ,label = "TKm, As")
plt.xscale('log')
plt.legend(loc="best")
plt.grid(True)
plt.ylim(0.9998,1.0002)
import matplotlib.ticker as ticker
#plt.ticklabel_format(style='sci',axis='x',scilimits=(0,0))
majorformatter = ticker.ScalarFormatter(useOffset =False)
ax = plt.gca()
ax.yaxis.set_major_formatter(majorformatter)

ax.xaxis.set_major_formatter(majorformatter)
#plt.ticklabel_format(useOffset =False)



In [13]:
psfrompkz3 = psu.powerspectrum(koverh= None, asciifile= pkfile3)
psfromtkz3 = psu.powerspectrum(koverh= None, asciifile= tkfile3,cosmo = M000)
psfromtksz3 = psu.powerspectrum(koverh= None, asciifile= tkfile3,cosmo = M000s)
plt.plot(psfrompkz3[0],psfrompkz3[1]/psfromtkz3[1])
plt.plot(psfrompkz3[0],psfrompkz3[1]/psfromtksz3[1])
plt.ylim(0.999,1.001)


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-13-b45dc7748ee3> in <module>()
----> 1 psfrompkz3 = psu.powerspectrum(koverh= None, asciifile= pkfile3)
      2 psfromtkz3 = psu.powerspectrum(koverh= None, asciifile= tkfile3,cosmo = M000)
      3 psfromtksz3 = psu.powerspectrum(koverh= None, asciifile= tkfile3,cosmo = M000s)
      4 plt.plot(psfrompkz3[0],psfrompkz3[1]/psfromtkz3[1])
      5 plt.plot(psfrompkz3[0],psfrompkz3[1]/psfromtksz3[1])

NameError: name 'pkfile3' is not defined
ERROR: NameError: name 'pkfile3' is not defined [IPython.core.interactiveshell]

Caption: For LCDM models, cb and matter power spectrum are the same thing. Reproducing the CAMB power spectra is limited only by our $\sigma_8$ integrals, which matches CAMB power spectra to an order similar to numerical noise. Whenever the amplitude of the power spectrum is supplied (ie. we do not have to calculate $\sigma_8$ to fix the amplitude, we get something that lies under the green curve. If we have to obtain $\sigma_8,$ we get something that lies under the blue curve.

Model with neutrinos (M000n1)


In [14]:
#Transfer function routines:

ntf = cio.cbtransfer(transferfile = ntkfile, Omegacdm = M000n1.Oc0 , Omegab = M000n1.Ob0)
ps = cio.matterpowerfromtransfersforsinglespecies(koverh = None, transfer = ntf, h = M000n1.h ,As = 4.1e-9 , ns = M000n1.ns)
#print np.shape(psfromtkncbs[1])
#print ps[:,1]/psfromtkncbs[1]

In [50]:
khval = np.logspace(-6,0,20)
print khval


[  1.00000000e-06   2.06913808e-06   4.28133240e-06   8.85866790e-06
   1.83298071e-05   3.79269019e-05   7.84759970e-05   1.62377674e-04
   3.35981829e-04   6.95192796e-04   1.43844989e-03   2.97635144e-03
   6.15848211e-03   1.27427499e-02   2.63665090e-02   5.45559478e-02
   1.12883789e-01   2.33572147e-01   4.83293024e-01   1.00000000e+00]

In [63]:
arr = np.zeros(len(khval))
for i, khv in enumerate(khval):
    arr[i] = psu.sigmasq(R =20.86 , ps = psfrompk, bgtype = "cb", cosmo = M000, khmin = khv)**0.5
print arr


[ 0.38570044  0.38570044  0.38570044  0.38570044  0.38570043  0.38570044
  0.38570043  0.38570044  0.38570044  0.38570041  0.38570009  0.38569443
  0.3856082   0.38451367  0.37504026  0.32968158  0.17704242  0.05089352
  0.01417638  0.00474895]

In [68]:
plt.plot(khval, arr/psu.sigmasq(R =20.86 , ps = psfrompk, bgtype = "cb", cosmo = M000,khmin = 5e-5)**0.5,'b')
plt.xscale('log')
plt.axvline(2*np.pi/1491, color= "k", label = "L = 1491Mpc/h")
plt.axvline(2*np.pi/256, color ="r",label = "L =256 Mpc/h")
plt.grid(True)
plt.ylim(0.999,1.01)
plt.xlabel("kmin (h/Mpc)")
plt.legend(loc="best")
plt.ylabel(r'$\frac{\sigma(R=20.86, kmin(Lbox))}{\sigma (R=20.86, kmin=0.97e-5)}$')


Out[68]:
<matplotlib.text.Text at 0x1078bc950>

In [25]:
psu.filterradiusformass(M = 1.0e15/0.71, bgtype = "cb", cosmo = M000)


Out[25]:
20.861478160859075

In [7]:
print psu.sigmaM(M = 8.453e+14/0.71 , ps = psfromtkncbscb, bgtype = "cb", cosmo = M000n1)

print psu.filterradiusformass(M = 3.47e12/0.71, bgtype = "cb", cosmo = M000n1)*0.71


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-7-7f48b65ecc21> in <module>()
----> 1 print psu.sigmaM(M = 8.453e+14/0.71 , ps = psfromtkncbscb, bgtype = "cb", cosmo = M000n1)
      2 
      3 print psu.filterradiusformass(M = 3.47e12/0.71, bgtype = "cb", cosmo = M000n1)*0.71

NameError: name 'psfromtkncbscb' is not defined
ERROR: NameError: name 'psfromtkncbscb' is not defined [IPython.core.interactiveshell]

In [17]:
print psu.sigmaM(M = 3.46882614991e+12/0.71 , ps = psfromtkncbscb, bgtype = "cb", cosmo = M000n1)*M000n1.growth(z= 2.018)[0]
print "growth ", M000n1.growth(z= 2.018) [0]


RinMpcoverh ***************
2.3011101463
0.699435820676
growth  0.447448519906

In [16]:
print M000n1.growth(z=2.018)[0]


0.447448519906

In [8]:
#Model With Neutrinos
psfrompkn = psu.powerspectrum(koverh = None, asciifile  = npkfile)#power spectrum output 
psfromtknm = psu.powerspectrum(koverh = None, asciifile  = ntkfile, pstype ="matter", cosmo = M000n1)#transfer output, tot
psfromtkncbs = psu.powerspectrum(koverh = None, asciifile  = ntkfile, pstype ="cb",sigma8type="matter", cosmo = M000n1) #transfer output cb, notmakized to sigma8
psfromtkncbscb = psu.powerspectrum(koverh = None, asciifile  = ntkfile, pstype ="cb",sigma8type="cb", cosmo = M000n1)#transfer function output cb, normalized to isgma8 for cb
psfromtkncbms = psu.powerspectrum(koverh = None, asciifile  = ntkfile, pstype ="cbmatter",sigma8type="matter", cosmo = M000n1)#transfer function output cb, normalized to isgma8 for cb

In [11]:
plt.plot(psfrompkn[0],psfromtknm[1]/psfrompkn[1],'k-',lw=2.0, label = "TKm, sigma8_m= 0.8/PS" )
plt.plot(psfrompkn[0],psfromtkncbs[1]/psfrompkn[1],'r-',label ="Tkcb, sigma8_m = 0.8/PS")
#plt.plot(psfrompkn[0],psfromtkncbs[1]/psfrompkn[1]/f,'b-',label ="TkcbA, sigma8_m = 0.8/PS")
plt.plot(psfrompkn[0],psfromtkncbms[1]/psfrompkn[1],'r--',lw =2.0, label ="Tkcbmatter, sigma8_m = 0.8/PS")
plt.plot(psfrompkn[0],psfromtkncbscb[1]/psfrompkn[1],'b',label ="Tkcb, sigma8_cb = 0.8/PS", lw = 2.0)
plt.plot(psfrompkn[0],3.704/4.107*psfromtkncbs[1]/psfrompkn[1],"r",ls= "dotted",lw =4.0,label ="Tkcb, sigma8_cb, amp")
#Cannot predict the change in As, but show that the red curve is expected if I know what As is
M000n1_rhom = M000n1.Oc0 +M000n1.Ob0 + M000n1.On0
M000n1_rhocb = M000n1.Oc0 + M000n1.Ob0
f = M000n1_rhom**2/M000n1_rhocb**2
plt.plot(psfrompkn[0],np.ones(len(psfrompkn[0]))*f,'g',ls='dashed', lw = 2)
#Assume rho_nu totally unclustered, then only contribution comes from 
#the denominator in delta rho/rho 
#Shown in green dashed line
plt.legend(loc="upper left")
plt.xscale('log')
plt.grid(True)
plt.xlabel(r'k $(h/Mpc)$')
plt.ylabel(r'ratio')


Out[11]:
<matplotlib.text.Text at 0x110090290>

Caption: Ratio of power spectra for M000n1 model at z = 0. The denominator is always the CAMB power spectrum output (matter power spectrum) at z = 0 of M000n1. Tkm (thick black curve) is the power spectrum calculated from the total transfer function and a sigma8 = 0.8, so this ratio should be 1. Tkcb, sigma8_m = 0.8 (thin red solid) is the cb power spectrum, normalized so that the total power spectrum has sigma8 = 0.8. The amplitude of fluctuations is the same as for the matter power spectrum, and the curve should be the same as the matter power spectrum at low k. At high k where the neutrinos are essentially unclustered, this cb power spectrum is expected to be ((rhocb + rhonu)/rhocb)^2 (green dashed line), and works asymptotically. The cb power spectrum (solid blue line), normalized so that the cb components have a $\sigma_8$ of 0.8 have a different amplitude of fluctuations from the matter power spectrum. Since Pcb ($\sigma_8(matter) =0.8$) is above the matter power spectrum, this amplitude is smaller than the amplitude of fluctuations for M000n1, but should have the same shape as Pcb ($\sigma_8 (matter) = 0.8)$ Pcb. The (blue-dashed) curve is obtained from the Pcb (sigma_8(m) =0.8) curve by multiplying the curve by ratio of the amplitude of the power spectrum for the solid blue curve (calculated) to the amplitude for the total matter power spectrum (calculated). Therefore we understand this entire set of plots

Suppression of Power spectrum due to Neutrinos


In [12]:
psfrompkint = psu.powerspectrum(koverh = psfromtknm[0], asciifile =tkfile,cosmo =M000s,sigmatype="cb")
plt.loglog(psfrompkint[0],psfrompkint[1],'k-',lw =2.0, label = "Matter, M000")
plt.loglog(psfromtknm[0],psfromtknm[1]*2.16/4.074,'b--',lw = 2.0, label = "M000n1, Amp =M000_{Amp}")
plt.grid(True)
plt.legend(loc="best")
plt.xlabel("k (h/Mpc)")
plt.ylabel("PK ($(Mpc/h)^{3}$)")
#print len(psfrompkint[0])
#print len(psfromtknm[0])
#print psfrompkint


Out[12]:
<matplotlib.text.Text at 0x105b0c790>
/usr/local/lib/python2.7/site-packages/matplotlib/scale.py:90: RuntimeWarning: invalid value encountered in less_equal
  mask = a <= 0.0

Caption : Shows the matter power spectrum computed for M000 and M00n1 with the amplitude adjusted to be the same as the amplitude for M000.


In [13]:
ncomp_ax1 = plt.subplot()
ncomp_ax1.plot(psfromtknm[0],psfromtknm[1]/psfrompkint[1], label   ="$\sigma_8=0.8$")
ncomp_ax1.plot(psfromtknm[0],psfromtkncbs[1]/psfrompkint[1], label ="$\sigma_8 (m) = 0.8$")
ncomp_ax1.plot(psfromtknm[0],psfromtkncbscb[1]/psfrompkint[1], label = "$\sigma_8 (cb) =0.8$")
ncomp_ax1.plot(psfromtknm[0],np.ones(len(psfromtknm[0]))*3.704/2.16,"g", ls = 'dashed')
ncomp_ax1.plot(psfromtknm[0],np.ones(len(psfromtknm[0]))*4.074/2.16,"g", ls = 'dashed')
#ncomp_ax1.plot(psfromtknm[0],np.ones(len(psfromtknm[0]))*0.8*f,"r", ls = 'dashed')
ncomp_ax1.legend(loc="lower left")
ncomp_ax1.set_ylabel("PK(M000n1)/PK(M000)")
ncomp_ax1.set_xlabel("k (h/Mpc)")
plt.grid(True)
ncomp_ax1.set_xscale('log')


Caption: Ratio of the linear power spectrum of the model M000n1 normalized in three different ways to the M000 power spectrum showing the suppression of the high k power spectrum relative to the low k par. The normalizations are chosen so that (black, solid) usual matter power spectrum, (blue, solid) cb power spectrum normalized so that the matter power spectrum has $\sigma_8 = 0.8,$ and (red, solid) cb power spectrum normalized so that the cb power spectrum has $\sigma_8 = 0.8$. The behavior of the three lines at low k is understood in terms of the amplitudes of fluctuations required for the normalizations. As seen in the ratio of different normalizations to the M000n1 matter power spectrum, the cb power spectrum normalized so that cb $\sigma_8 = 0.8 $ has a smaller amplitude and should be below the the cb power spectrum normalized so that cb power spectrum has a $\sigma_8 = 0.8$. Thus the red curve is always below the blue curve with a constant ratio. The low k asymoptotics is given by the relative amplitudes and the expected values are shown by the green dashed curves. The relative behavior of the black, blue, and red curves can be better understood from the previous curve showing the ratios. However, the supression of any of these curves with respect to their low k values cannot be understood in terms of such numbers which are post-processing/counting issues, while the supression is a real effect due to lack of sourcing of the Poisson equation. In comparison, the effect of sourcing due to the clustering of neutrinos is negligible as seen in the power spectrum paper.

Power Spectra at different redshifts:


In [41]:
#M000files
#transfer and power spectrum files at redshift 1, 2, 3
tkfile1 = dirn +"example_data/M000/m000n0_transfer_247_out.dat"
pkfile1 = dirn +"example_data/M000/m000n0_matterpower_247.dat"
pkfile2 = dirn +"example_data/M000/m000n0_matterpower_163.dat"
pkfile3 = dirn + "example_data/M000/m000n0_matterpower_121.dat"
tkfile3 = dirn + "example_data/M000/m000n0_transfer_121_out.dat"

#M000n1 files
#power spectrum file at redshift  1, 2, 3
ntkfile1 = dirn +"example_data/M000n1/m000n1_transfer_247_out.dat"
npkfile1 = dirn +"example_data/M000n1/m000n1_matterpower_247.dat"
npkfile2 = dirn +"example_data/M000n1/m000n1_matterpower_163.dat"
npkfile3 = dirn + "example_data/M000n1/m000n1_matterpower_121.dat"


psfromtk1 = psu.powerspectrum(koverh = None, asciifile = tkfile1, cosmo = M000)
psfrompk1 = psu.powerspectrum(koverh = None, asciifile = pkfile1, cosmo = M000)
psfromtk10 = psu.powerspectrum(koverh = None, asciifile = tkfile, cosmo = M000s,z =1.0,method= "CAMBoutgrowth")
psfromtk2  = psu.powerspectrum(koverh = None, asciifile = pkfile, cosmo = M000,z =2.0,method = "CAMBoutgrowth")
psfromtk2A = psu.powerspectrum(koverh = None, asciifile = tkfile, cosmo = M000,z =2.0,method = "CAMBoutgrowth")

M000psfrompk1  = psu.powerspectrum(koverh = None, asciifile = pkfile1)
M000psfrompk2  = psu.powerspectrum(koverh = None, asciifile = pkfile2)
M000psfrompk3  = psu.powerspectrum(koverh = None, asciifile = pkfile3)
#M000psfromtk10 = psu.powerspectrum(koverh = None, asciifile = tkfile, cosmo = M000s, z = 1.0, method  = "CAMBoutgrowth")
#M000psfromtk20 = psu.powerspectrum(koverh = None, asciifile = tkfile, cosmo = M000s, z = 2.0, method  = "CAMBoutgrowth")
#M000psfromtk30 = psu.powerspectrum(koverh = None, asciifile = tkfile, cosmo = M000s, z = 3.0, method  = "CAMBoutgrowth")
M000psfromtk10 = psu.powerspectrum(koverh = None, asciifile = tkfile, cosmo = M000, z = 1.006, method  = "CAMBoutgrowth")
M000psfromtk20 = psu.powerspectrum(koverh = None, asciifile = tkfile, cosmo = M000, z = 2.02, method  = "CAMBoutgrowth")
M000psfromtk30 = psu.powerspectrum(koverh = None, asciifile = tkfile, cosmo = M000, z = 3.04, method  = "CAMBoutgrowth")


M000n1psfrompk1  = psu.powerspectrum(koverh = None, asciifile = npkfile1)
M000n1psfrompk2  = psu.powerspectrum(koverh = None, asciifile = npkfile2)
M000n1psfrompk3  = psu.powerspectrum(koverh = None, asciifile = npkfile3)
M000n1psfromtk10 = psu.powerspectrum(koverh = None, asciifile = ntkfile, cosmo = M000n1, z = 1.006, method  = "CAMBoutgrowth")
M000n1psfromtk20 = psu.powerspectrum(koverh = None, asciifile = ntkfile, cosmo = M000n1, z = 2.02, method  = "CAMBoutgrowth")
M000n1psfromtk30 = psu.powerspectrum(koverh = None, asciifile = ntkfile, cosmo = M000n1, z = 3.04, method  = "CAMBoutgrowth")

mykoverh = M000n1psfrompk1[0] 
iM000psfrompk1  = psu.powerspectrum(koverh = mykoverh, asciifile = npkfile1)
iM000psfrompk2  = psu.powerspectrum(koverh = mykoverh, asciifile = npkfile2)
iM000psfrompk3  = psu.powerspectrum(koverh = mykoverh, asciifile = npkfile3)
iM000psfromtk10 = psu.powerspectrum(koverh = mykoverh, asciifile = tkfile, cosmo = M000s, z = 1.0, method  = "CAMBoutgrowth")
iM000psfromtk20 = psu.powerspectrum(koverh = mykoverh, asciifile = tkfile, cosmo = M000s, z = 2.0, method  = "CAMBoutgrowth")
iM000psfromtk30 = psu.powerspectrum(koverh = mykoverh, asciifile = tkfile, cosmo = M000s, z = 3.0, method  = "CAMBoutgrowth")

Growth Functions:

The growth function for the models are different. M000 has a scale independent growth function, while M000n1 has a scale-dependent growth function whose low k limit approximately (Note that the M000 model still has 3.04 massless neutrino species, while the massive neutrino model does not ... it was run with 0.046 massless neutrinos, which would better have been 0. with a higher temperature) coincide with M000 (for small redshifts, where the massive neutrinos behave as CDM). In the high k limit, the neutrinos do not contribute to the clustering.


In [14]:
z = np.arange(0.00,3.1,0.1) 
gfig, gax0, gax1 = pu.settwopanel(setdifflimits =[0.9,1.1])
gax0.plot ( z, M000.growth(z) [0], 'k-', lw =2.0, label = "M000")
gax0.plot (z, M000n1.growth(z)[0], 'r--',label = "M000n1", lw =2 )
gax0.plot(z , 1./(1+z), 'b.',label = "a")
gax0.set_ylabel("D(z)")
gax0xticks  = gax0.axes.get_xticklabels()
gax0.set_xticklabels(gax0xticks[::2], visible= False)
gax0.legend(loc="best")
gax1.plot(z, (1./M000.growth(z)[0]*M000n1.growth(z)[0]),'r--',lw =2.0)
gax1.plot(z, 1./M000.growth(z)[0]/(1.+z), 'b.')
gax1ticks = gax1.axes.get_yticklabels()
gax1.set_xlabel("z")
gax1.yaxis.set_ticks([0.9,0.95,1.0, 1.05, 1.1])
gax1.yaxis.set_ticklabels([0.9,"",1.0,"",1.1])
#gax1.set_yticklabels([0.95,1.0], visible= True)
#print len(gax1ticks)
#plt.setp(gax1ticks[::,2],visible=True)
#


Out[14]:
[<matplotlib.text.Text at 0x112acdf10>,
 <matplotlib.text.Text at 0x112acbf50>,
 <matplotlib.text.Text at 0x112a9f290>,
 <matplotlib.text.Text at 0x112acd550>,
 <matplotlib.text.Text at 0x112acb210>]

In [40]:
zadrian = np.arange(0.,11.,1.0)
print zadrian
print M000.growth(z = zadrian) [0]
print M000n1.growth(z = zadrian)[0]


[  0.   1.   2.   3.   4.   5.   6.   7.   8.   9.  10.]
[ 1.          0.6249225   0.43284744  0.32803081  0.26346592  0.21997139
  0.1887494   0.16527179  0.14698305  0.13233882  0.12034967]
[ 1.          0.63866592  0.44987636  0.34528185  0.28009808  0.23577741
  0.20371769  0.17945041  0.16043804  0.14513584  0.13255041]

Caption : The growth function of the model with neutrinos is higher than the growth function in a similar model without neutrinos: in order to obtain the same amplitude of a fluctuation at z = 1.0, the fluctuation should have been larger at a higher z for the model with neutrinos, because with the absence of some dark matter the fluctuation will grow less. The ratio reaches about 5 percent at z = 3.0. We will see this effect on the high redshift power spectra obtained from the z = 0, power spectra scaled by the growth function. The blue dotted curve shows the growth function that would have been obtained in a massless neutrino EdS model with no dark energy, showing that dark energy slows the clustering of matter appreciably.


In [51]:
fig = plt.figure(figsize = (16,6))

ax = plt.subplot(1,2,1)
ay = plt.subplot(1,2,2)


ax.yaxis.set_major_formatter(majorformatter)
ax.plot(M000psfrompk1[0],M000psfromtk10[1]/M000psfrompk1[1])
ax.plot(M000psfrompk2[0],M000psfromtk20[1]/M000psfrompk2[1])
ax.plot(M000psfrompk2[0],M000psfromtk30[1]/M000psfrompk3[1])
#ax.plot(M000psfrompk1[0],M000psfromtk10[1]/M000psfrompk1[1])
#ax.plot(M000psfrompk2[0],M000psfromtk20[1]/M000psfrompk2[1])
#ax.plot(M000psfrompk2[0],M000psfromtk30[1]/M000psfrompk3[1])
ax.grid(True)
ax.set_xlabel("k (h/Mpc)")
ax.set_ylabel("PK(z=0)$D^2$(z)/PK(z=z)")
ax.set_ylim(0.995,1.005)
ax.set_title("M000")
ax.set_xscale('log')
#
ay.plot(M000n1psfrompk1[0],M000n1psfromtk10[1]/M000n1psfrompk1[1],label = "z = 1")
ay.plot(M000n1psfrompk2[0],M000n1psfromtk20[1]/M000n1psfrompk2[1],label = "z = 2")
ay.plot(M000n1psfrompk2[0],M000n1psfromtk30[1]/M000n1psfrompk3[1],label = "z = 3")
#ay.plot(M000n1psfrompk1[0],M000n1psfromtk10[1]/M000n1psfrompk1[1])
#ay.plot(M000n1psfrompk2[0],M000n1psfromtk20[1]/M000n1psfrompk2[1])
#ay.plot(M000n1psfrompk2[0],M000n1psfromtk30[1]/M000n1psfrompk3[1])
ay.set_xscale('log')
ay.legend(loc="best")
ay.set_title("M000n1")
ay.grid(True)
ay.set_ylim(ymin = 0.99)

#az.plot(M000n1psfrompk1[0],M000n1psfromtk10[1]/M000n1psfrompk1[1])
#az.plot(M000n1psfrompk2[0],M000n1psfromtk20[1]/M000n1psfrompk2[1])
#az.plot(M000n1psfrompk2[0],M000n1psfromtk30[1]/M000n1psfrompk3[1])
#az.plot(M000n1psfrompk1[0],M000n1psfromtk10[1]/iM000psfrompk1[1])
#az.plot(M000n1psfrompk2[0],M000n1psfromtk20[1]/iM000psfrompk2[1])
#az.plot(M000n1psfrompk2[0],M000n1psfromtk30[1]/iM000psfrompk3[1])
#az.set_xscale('log')
#az.set_ylim(1.,1.0.3)
#az.grid(True)


Out[51]:
(0.99, 1.1200000000000001)

Caption : The ratio of power spectra for each model at different redshifts calculated in two different ways. The numerator is obtained from the z = 0 power spectrum by multiplying by the square of the scale independent growth function plotted above. The denominator is obtained from CAMB by evolving a set of coupled equations obtaining scale dependent growth with neutrinos sourcing the growth at low k, but hardly at high k. The little difference (about 2 percent for z = 3 in the LCDM model) at high k and about 0.5 percent


In [49]:
fig = plt.figure(figsize = (16,6))
ay = plt.subplot(1,2,2)
az = plt.subplot(1,2,1)
ay.plot(M000n1psfrompk1[0],M000n1psfromtk10[1]/M000n1psfrompk1[1])
ay.plot(M000n1psfrompk2[0],M000n1psfromtk20[1]/M000n1psfrompk2[1])
ay.plot(M000n1psfrompk2[0],M000n1psfromtk30[1]/M000n1psfrompk3[1])
#ay.plot(M000n1psfrompk1[0],M000n1psfromtk10[1]/M000n1psfrompk1[1])
#ay.plot(M000n1psfrompk2[0],M000n1psfromtk20[1]/M000n1psfrompk2[1])
#ay.plot(M000n1psfrompk2[0],M000n1psfromtk30[1]/M000n1psfrompk3[1])
ay.set_xscale('log')
ay.set_ylim(0.995,1.005)
ay.set_xlim(xmin = 0.01)
ay.grid(True)

def gsratio(redshift):
        redshift = np.asarray(redshift)
        gg = M000n1.growth(z= redshift)[0]
        gs = M000.growth(z = redshift) [0]
        
        return (gg/gs)**2 
az.plot(M000n1psfrompk1[0],M000n1psfromtk10[1]/M000n1psfrompk1[1]/gsratio(1.006))
az.plot(M000n1psfrompk2[0],M000n1psfromtk20[1]/M000n1psfrompk2[1]/gsratio(2.02))
az.plot(M000n1psfrompk2[0],M000n1psfromtk30[1]/M000n1psfrompk3[1]/gsratio(3.04))
#az.plot(M000n1psfrompk1[0],M000n1psfromtk10[1]/M000n1psfrompk1[1])
#az.plot(M000n1psfrompk2[0],M000n1psfromtk20[1]/M000n1psfrompk2[1])
#az.plot(M000n1psfrompk2[0],M000n1psfromtk30[1]/M000n1psfrompk3[1])
az.set_xscale('log')
az.set_xlim(xmax =0.01)
az.set_ylim(0.995,1.005)
az.grid(True)
az.yaxis.set_major_formatter(majorformatter)
ay.yaxis.set_major_formatter(majorformatter)


Caption : The ratio of power spectra for M000n1 at different redshifts calculated in two different ways. The large ratio (> 10 percent) at low k is multiplied by the ratio of the square of growth factors of M000 and M000n1 (which is around 5 percent, shown in a figure), to effectively use the growth function of M000 instead of M000n1. At these high k, this is appropriate as the growth is sourced by both neutrinos and matter.

Mass Function Results from Fitting Formulae


In [26]:
Masses = np.logspace(8,15,100)
psfrompkn = psu.powerspectrum(koverh = None, asciifile  = npkfile)#power spectrum output 
psfromtknm = psu.powerspectrum(koverh = None, asciifile  = ntkfile, pstype ="matter", cosmo = M000n1)#transfer output, tot
psfromtkncbs = psu.powerspectrum(koverh = None, asciifile  = ntkfile, pstype ="cb",sigma8type="matter", cosmo = M000n1) #transfer output cb, notmakized to sigma8
psfromtkncbscb = psu.powerspectrum(koverh = None, asciifile  = ntkfile, pstype ="cb",sigma8type="cb", cosmo = M000n1) #transfer output cb, notmakized to sigma8
sigM = psu.sigmaM(M =Masses,ps= psfrompk ,cosmo = M000)
dndlnM_M000 = psu.dndlnM(M =Masses,ps= psfrompk ,cosmo = M000)
dndlnM_M000nm = psu.dndlnM(M =Masses,ps= psfromtknm ,cosmo = M000n1)
dndlnM_M000ncbs = psu.dndlnM(M =Masses,ps= psfromtkncbs ,cosmo = M000n1)
dndlnM_M000ncbscb = psu.dndlnM(M =Masses,ps= psfromtkncbscb ,cosmo = M000n1)
#plt.plot(Masses, dndlnM_M000,'k-',lw=2.0)
plt.plot(Masses, dndlnM_M000/dndlnM_M000nm,'k-',lw=2.0, label = "M000")
plt.plot(Masses, dndlnM_M000ncbs/dndlnM_M000nm,'r-',lw=2.0,label = "cb, $\sigma_8(m)$")
plt.plot(Masses, dndlnM_M000ncbscb/dndlnM_M000nm,'b--',lw=2.0,label = "cb, $\sigma_8(cb)$")
plt.xscale('log')
plt.ylabel("dndlnM (pstype)/dndlnM(cbn)")
plt.ylim(0.8,1.2)
plt.grid(True)
plt.legend(loc="best")


0.333 0.788 0.807 1.795
0.333 0.788 0.807 1.795
0.333 0.788 0.807 1.795
0.333 0.788 0.807 1.795
Out[26]:
<matplotlib.legend.Legend at 0x114e1d1d0>

In [27]:
plt.plot(Masses, dndlnM_M000/dndlnM_M000ncbs,'k-',lw=2.0, label = "M000")
plt.plot(Masses, dndlnM_M000ncbs/dndlnM_M000ncbs,'r-',lw=2.0,label = "cb, $\sigma_8(m)$")
plt.plot(Masses, dndlnM_M000ncbscb/dndlnM_M000ncbs,'b--',lw=2.0,label = "cb, $\sigma_8(cb)$")
plt.xscale('log')
plt.ylabel("dndlnM (pstype)/dndlnM(cbn)")
plt.ylim(0.8,1.2)
plt.grid(True)
plt.legend(loc="best")


Out[27]:
<matplotlib.legend.Legend at 0x114e372d0>

In [ ]:
sigMn = psu.sigmaM(M =Masses,ps= psfrompkn ,cosmo = M000n1)
plt.plot(Masses, sigM,'b-',label = "M000")
plt.plot(Masses, sigMn,'r--',label = "M000n1")
plt.grid(True)
plt.xscale('log')
plt.xlabel("Mass $M_\odot$")
plt.ylabel("$\sigma (M)$")
plt.legend(loc="best")

In [21]:
mf_fig, mf_ax0, mf_ax1 = pu.settwopanel(setdifflimits=[0.85,1.15])
mf_ax0.loglog(Masses,dndlnM_M000, 'k-',lw=2.0, label = "M000")
mf_ax0.loglog(Masses,dndlnM_M000n1, 'r--',lw =2.0,label = "M000n1")
mf_ax0.xaxis.set_ticklabels("",visible=False)
mf_ax1.plot(Masses, dndlnM_M000n1/dndlnM_M000)
#mf_ax1.yaxis.set_ticks(np.arange(0.9,1.15,0.05))
mf_ax1.set_yticks([0.9,1.0,1.10])
#mf_ax1.yaxis.set_ticklabels([0.9,"",1.0,"",1.10])
mf_ax1.set_ylim(0.9,1.1)
mf_ax1.set_xscale('log')
mf_ax1.set_xlabel("M $M_\odot$")
mf_ax0.set_ylabel("dn/dln M $(h Mpc)^{-3}$")
mf_ax0.set_title("Check y units")


Out[21]:
<matplotlib.text.Text at 0x115f16f10>

In [22]:
plt.plot(Masses, dndlnM_M000n1/dndlnM_M000)
plt.xscale('log')
plt.grid(True)
plt.title("Mass Function Suppression Due to Neutrinos using matter power spectrum with neutrinos")


Out[22]:
<matplotlib.text.Text at 0x111c6d4d0>

Mass functions for different neutrino power spectra


In [23]:
plt.loglog(Masses, psu.dndlnM(Masses, ps = psfrompkn , cosmo = M000n1))
plt.loglog(Masses, psu.dndlnM(Masses, ps = psfromtknm, cosmo = M000n1))
plt.loglog(Masses, psu.dndlnM(Masses, ps = psfromtkncbs, cosmo = M000n1))
plt.loglog(Masses, psu.dndlnM(Masses, ps = psfromtkncbscb, cosmo = M000n1))


0.333 0.788 0.807 1.795
0.333 0.788 0.807 1.795
0.333 0.788 0.807 1.795
0.333 0.788 0.807 1.795
Out[23]:
[<matplotlib.lines.Line2D at 0x115e3e690>]

Mass Functions for Different Normaliztions


In [24]:
#plt.plot(Masses, psu.dndlnM(Masses, ps = psfrompkn , cosmo = M000n1)/psu.dndlnM(Masses, ps = psfrompk, cosmo= M000), lw = 6,label = "Pk neutrino")
plt.plot(Masses, psu.dndlnM(Masses, ps = psfromtknm, cosmo = M000n1)/psu.dndlnM(Masses, ps = psfrompkn, cosmo= M000n1), label = "TK/PKn")
plt.plot(Masses, psu.dndlnM(Masses, ps = psfromtkncbs, cosmo = M000n1)/psu.dndlnM(Masses, ps = psfrompkn, cosmo= M000n1),label = "TKcb/PKn", lw =2, color = 'b', ls = 'dotted')
plt.plot(Masses, psu.dndlnM(Masses, ps = psfromtkncbscb, cosmo = M000n1)/psu.dndlnM(Masses, ps = psfrompkn, cosmo= M000n1),color = "r", lw =2,ls = 'dashed',label = "TKcb $\sigma_8$ cb/PKn")
plt.xscale('log')
plt.grid(True)
plt.legend(loc= "best")
plt.ylim(0.9,1.1)
plt.xlim(1.0e8,1.0e15)
plt.xlabel(r'Mass $M_\odot$')
plt.ylabel('ratio')


0.333 0.788 0.807 1.795
0.333 0.788 0.807 1.795
0.333 0.788 0.807 1.795
0.333 0.788 0.807 1.795
0.333 0.788 0.807 1.795
0.333 0.788 0.807 1.795
Out[24]:
<matplotlib.text.Text at 0x1124f3d10>

Caption: Ratio of Mass function of M000n1 at z = 0, using power spectrum in different ways. The denominator always uses the CAMB matter power spectrum output of M000n1 and is called Pkn. The black line shows the ratio of mass function when the power spectrum is calculated from the transfer function and the sigma8 value. This should be $1.0$ and any deviations are simply due to the sigma8 calculation differing from CAMB. The blue dotted curve uses the cb part of the power spectrum, with the same sigma8 normalization, so each component has the same primordial amplitude as those used in the black curve. The red dashed curve shows the cb power spectrum, normalized so that the cb power spectrum is normalized to sigma8 =0.8. The power spectrum is the same as the matter power spectrum for this model at low k and higher by a factor of (rhom/rhocb)^2 at high k.

How Does the Mass function of these different types evolve with redshift


In [27]:
#plt.plot(Masses, psu.dndlnM(Masses, ps = psfrompkn , cosmo = M000n1)/psu.dndlnM(Masses, ps = psfrompk, cosmo= M000), lw = 6,label = "Pk neutrino")
#psfromtkn2  = psu.powerspectrum(koverh = None, asciifile = ntkfile, cosmo = M000n1,z =2.0,method = "CAMBoutgrowth", pstype ='cb', sigma8type= "cb")
import halomass as hmf

psfromtk2  = psu.powerspectrum(koverh = None, asciifile = tkfile, cosmo = M000,z =2.0,method = "CAMBoutgrowth")
psfromtkn2 = psu.powerspectrum(koverh = None,z = 2.0, asciifile  = ntkfile, pstype ="cb",sigma8type="cb", cosmo = M000n1,method="CAMBoutgrowth")#transfer function output cb, normalized to isgma8 for cb
plt.plot(Masses, psu.dndlnM(Masses, ps = psfromtk2, cosmo = M000)/psu.dndlnM(Masses, ps = psfrompk, cosmo= M000), label = "TK 2/PK 0")
plt.plot(Masses, psu.dndlnM(Masses, ps = psfromtkn2, cosmo = M000n1)/psu.dndlnM(Masses, ps = psfromtkncbscb, cosmo= M000n1),label = "TKcb 2/PKn 0", lw =2, color = 'r', ls = 'dashed')
plt.plot(Masses, hmf.dndlnM0(Masses,ps = psfromtkms,cosmo = M000s, z= 2.0)/hmf.dndlnM0(Masses, ps = psfromtkms, cosmo = M000s, z =0.0) ,label="Suman")

#plt.plot(Masses, psu.dndlnM(Masses, ps = psfromtkncbscb, cosmo = M000n1)/psu.dndlnM(Masses, ps = psfromtkn2, cosmo= M000n1),color = "r", lw =2,ls = 'dashed',label = "TKcb $\sigma_8$ cb/PKn")
plt.xscale('log')
plt.grid(True)
plt.legend(loc= "best")
#plt.ylim(0.8,1.2)
plt.xlim(1.0e11,1.0e15)
plt.xlabel(r'Mass $M_\odot$')
plt.ylabel('ratio')


0.333 0.788 0.807 1.795
0.333 0.788 0.807 1.795
0.333 0.788 0.807 1.795
0.333 0.788 0.807 1.795
0.295094350742 0.779390315288 0.807 1.795
0.333 0.788 0.807 1.795
Out[27]:
<matplotlib.text.Text at 0x113ba1310>

Caption: Evolution of the mass function with neutrinos is not tremendously different from the evolution of mass function without neutrinos. The curves show the ratio mass at z = 2.0 (for an exaggerated effect which is still tiny), to the mass function at z = 0.0 for both the M000 model and M000n1 model. For the M000n1 model, we use the cb power spectrum normalized so that the cb sigma8 = 0.8


In [ ]: