Commit c6a16465 authored by Hannah Middleton's avatar Hannah Middleton

updating cosomlogy for plotting script

parent 97afc357
import numpy as np import numpy as np
from scipy.integrate import quad from scipy.integrate import quad
import astropy.units as u
from astropy.cosmology import Planck18 as cosmo from astropy.cosmology import Planck18 as cosmo
import readPosterior
def zInt(beta,gamma,zlow,zhigh): def zInt(beta,gamma,zlow,zhigh):
OmegaM = 0.3 OmegaM = cosmo.Om0
OmegaL = 0.7 OmegaL = cosmo.Ode0
Ho = 1.0/13.8
# Ho from Planck18, which we convert to units of 1/Gyrs
Gyr_in_s = 1.e9 * 1.*u.year.to(u.s)
Mpc_in_m = 1.e6 * 1.*u.parsec.to(u.m)
H0_in_km_per_s_per_Mpc = cosmo.H(0).value
H0_in_per_Gyr = H0_in_km_per_s_per_Mpc * Gyr_in_s / (Mpc_in_m / 1000)
Ho = H0_in_per_Gyr
zfn = lambda zi: ( (1.0 + zi)**(beta-1.0) * np.exp(zi/gamma) ) / \ zfn = lambda zi: ( (1.0 + zi)**(beta-1.0) * np.exp(zi/gamma) ) / \
( Ho * ( OmegaM*(1.0+zi)**3. + OmegaL )**(1./2.) ) ( Ho * ( OmegaM*(1.0+zi)**3. + OmegaL )**(1./2.) )
...@@ -46,18 +53,19 @@ def dndlog10M(theta,log10Ms): ...@@ -46,18 +53,19 @@ def dndlog10M(theta,log10Ms):
def main(): def main():
pathToRuns = '../../runs/simpleModel/logNormLikeMstar6to10/' pathToRuns = '../../runs/simpleModel/logNormLikeMstar6to10/'
pathToRuns = '../models/agnostic_model/posterior.dat'
nRuns=5 nRuns=5
simpleModelData = readPosterior.readPosterior(pathToRuns, nRuns=nRuns) simpleModelData = np.genfromtxt(pathToRuns,names=True)#readPosterior.readPosterior(pathToRuns, nRuns=nRuns)
# set up M range and write to file # set up M range and write to file
log10MRange = np.linspace(6.0, 11.0, 40) log10MRange = np.linspace(6.0, 11.0, 40)
writeMs = open('{}combined/ms.dat'.format(pathToRuns),'w') writeMs = open('tmp2_ms.dat'.format(pathToRuns),'w')
for logM in log10MRange: for logM in log10MRange:
writeMs.write('{}\n'.format(10.0**logM)) writeMs.write('{}\n'.format(10.0**logM))
writeMs.close() writeMs.close()
# place to write results # place to write results
writeLines = open('{}combined/mLines.dat'.format(pathToRuns),'w') writeLines = open('tmp2_mLines.dat'.format(pathToRuns),'w')
# compute lines # compute lines
for i in range(len(simpleModelData)): for i in range(len(simpleModelData)):
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment