cosmology.sm

__init_cosmology
	# set global variables for the cosmology 
	define Omegam 0.30711		# matter
	define Omegal 0.69289		# dark energy
	define Omegab 0.04825		# baryons
	define sigma8 0.8288		# mass variance
	define hubble0 0.6777		# reduced Hubble constant 
	define Omegan 0				# neutrinos
	define mstar 4.944e+12		# typical non-linear mass [M_sun]
	define deltac 1.675748522	# critical density for collapse
	## neutrino physics
	#define f_nu ($Omegan/$Omegam)
	#define m_nu ($Omegan*93.14*$hubble0**2)	# neutrino masses [eV]
	#define z_nr (1890*$m_nu-1)					# non-relativistic transition
	#define k_nr (0.018*($Omegam*$m_nu)**(1/2))	# non-relativistic scale [h/Mpc]

rhocz 1
	# return the critical density as a function of redshift
	local set z $1
	set $0 = $rho_c *($Omegam*(1+z)**3+$Omegal+(1-$Omegam-$Omegal)*(1+z)**2)

deltavz 1
	# return the virial overdensity as a function of redshift
	# Bryan & Norman (1998) approximation
	local set z $1
	local set omx = $Omegam*(1+z)**3/($Omegam*(1+z)**3+$Omegal+(1-$Omegam-$Omegal)*(1+z)**2) - 1
	set $0 = (18*PI**2+82*omx-39*omx**2)

delta_c 12
	# : delta_c z [inverse=0]
	# interpolate linear critical overdensity from table
	# needs variable dczfile to be set with the path to a file with 
	# an overdensity - redshift table
	local define reverse 0
	if ($?2) {define reverse $2}
	set delta local
	set redshift local
	data "$!dczfile"
	read {delta 1 redshift 2}
	if (!$reverse) {
	  interp2 redshift delta $1 $0
	} else {
	  interp2 delta redshift $1 $0
	}

lmspk 12
	# : lmspk M [inverse=0]
	# interpolate mass variance from table
	# needs variable lmsfile to be set with the path to a file with 
	# a log(mass) - power spectrum table
	set lms local
	set spk local
	data "$!lmsfile"
	read {lms 1 spk 2}
	if (!$?2) {
		local set lm = lg($1)
		interp2 lms spk lm $0
	} else {
		sort {spk lms}
		interp2 spk lms $1 $0
	}

omz 34
	# : omz z mass f [z0=0]
	# return scaled variable omega from redshift, mass and mass fraction
	local define f $3
	local set z0 = 0
	if ($?4) {set z0 = $4}
	local set df = delta_c($1)
	local define d0 $(delta_c(z0))
	local set Sf = lmspk($f*$2)
	local set S0 = lmspk($2)
	set $0 = (df - $d0)/sqrt(Sf - S0)

growthf 1
	# growth factor as a function of redshift
	# only works for a flat ΛCDM cosmology
	local set z = $1
	local set a = 1/(1+z)
	local set Bottom = a + $Omegam*(1-a) + $Omegal*(a**3 - a)
	local set Omega = $Omegam/Bottom		# Omega(a)
	local set El = (a**3)*$Omegal/Bottom
	local set gb = Omega**(4/7) - El		# Lahav & Suto 2004
	local set gbottom = gb + (1 + 1/2*Omega)*(1 + El/70)
	local set ga = 2.5*Omega/gbottom		# relative growth factor
	set $0 = a*ga							# D(a) not normalized!

gnfw 1
	# : gnfw r/rvir
	# integrated nwf profile
	set $0=ln(1+$1)-$1/(1+$1)

rhom_nfw 2
	# : rhom_nfw r/rvir c
	# encloded density inside an nfw profile
	local define z0 0
	local define Dv (deltavz($z0))
	set $0 = 1/$1**3 * gnfw($2*$1)/gnfw($2)

massr_nwf 2
	# : massr_nwf r/rvir c
	# enclosed mass inside an nfw profile
	set $0 = gnfw($2*$1)/gnfw($2)