# Cosmology Calculator

#### Set Variables:

Welcome to ICRAR's Cosmology Calculator! This calculator was written by Aaron Robotham and Joseph Dunne in the programming language R, and uses the library Shiny to provide the interface. The functions are available as part of Aaron Robotham's R celestial package.

If this website is used to assist published work, please could you reference it by including a link to the front page. Also, please do get in touch with Aaron Robotham if you have any problems, feedback or ideas at aaron.robotham@uwa.edu.au

For clarity this website uses the convention of a prefix c, p or l for comoving, physical (or proper) and luminosity distances respectively.

#### Description of key cosmological parameters used in calculator

z - Redshift, where z must be > -1

H Hubble constant as defined at z=0 / (km/s)/Mpc

OmegaM - Omega matter, the relative component of energy in mass, as defined at z=0

OmegaL - Omega lambda, the relative component of energy in dark energy, as defined at z=0

OmegaR - Omega radiation, the relative component of energy in radiation (including neutrinos), as defined at z=0

Sigma8 - Amplitude of the linear power spectrum on the scale of 8 Mpc/h, as defined at z=0

w - The value of dark energy equation of state at z=0, as defined at z=0

w' - The evolution term that governs how the dark energy equation of state evolves with redshift

Re the dark energy equation of state evolution parameters w and w'. The dark energy EoS evolution is parameterised such that P=w.rho.c^2 and w=w+2.w'.(1-1/(1+z)). This means rhoDE=rhoDE.e^{-6.w'.(1-1/(1 + z))}/(1 + z)^{-(3 + 3.w_0 + 6.w')}, as per Wright (2006).

#### Calculation

This tab is used to calculate various distance parameters.

To use this tab, fill in the variables under Set Variables and click the button to calculate variables at a certain redshift. Type '1-OmegaM' into the OmegaL field to set OmegaL to 1-OmegaM for all calculations. The Reference Set box lets you set the variables (H, OmegaM, OmegaL, OmegaR and Sigma8) to a reference set (e.g. Planck, WMAP). The fSigma8 flag varies whether the exact growth rate (f, when fSigma8 is unchecked) or the growth rate times Sigma8 (fSigma8, when fSigma8 ic checked) is computed.

Under Custom Calc, the calculation may be done using a chosen variable from the Variable menu. When the Value box contains a value, the custom calculation will be used next time the button is clicked. For ambiguous mappings (e.g. Angular Size has two possible redshift solutions for many values) you can specify whether the lower redshift or higher redshift solution is found by selecting the respective option from the 'Solution to find' option box.

#### Plot

This tab is used to plot various parameters.

To use this tab, fill in the variables under Set Variables and Plot Options, and click to produce plots across a redshift range. The first plot is a distance plot, and the second plot is a custom plot which may be modified using the options under Custom Plot. Type '1-OmegaM' into the OmegaL field to set OmegaL to 1-OmegaM for all calculations. The Reference Set box lets you set the variables (H, OmegaM, OmegaL, OmegaR and Sigma8) to a reference set (e.g. Planck, WMAP). The fSigma8 flag varies whether the exact growth rate (f, when fSigma8 is unchecked) or the growth rate times Sigma8 (fSigma8, when fSigma8 ic checked) is computed. Some of the plot options are as follows:

z Start - The starting redshift for the plots.

z End - The finishing redshift for the plots.

z Resolution - The number of points to plot between z Start and z End. The spacing is made in equal separations of log10 expansion factor (where a=1/1+z).

The resulting data can be saved using the options under Save Data.

#### Survey Design

This tab is used to produce multiple diagnostics of an area in the sky.

To use this tab, fill in the variables under Set Variables and click the button to calculate the Comoving Volume. The Reference Set box lets you set the variables (H, OmegaM, OmegaL and OmegaR) to a reference set (e.g. Planck, WMAP). If the area is unknown, an extra option under Find Area (optional) can be used the find the area of the sky given the latitude and longitude. The Area field is then updated by clicking the button. The Cosmic Variance options allows you to control whether the variance is determined from the area stated in the Area box, or whether it ignores this and uses the longitude and latitude geometry in the Find Area section of the web form.

The Regions button affects a number of quantities on the results page, but this should be clear from the parenthetic for all regions text. Quantities like volume and the number of halos etc are multiplied by the number of regions. Quantities regarding distances are only computed within a single region, since the geometry between regions cannot be specified, merely the number. Remember to hit , or your changes will not affect the page output!

#### Sky Separation

This tab is used to find various types of separation between two points at (potentially) different locations on the sky (separation given by Sep / deg ) and at different redshifts ( z1 and z2 ). It is only in the scenario of both of these being different that the calculation becomes non-trivial. To calculate these outputs I relied heavily on Peacock (1999) and to a lesser extent Liske (2000), but any mistakes are of course mine.

Whilst this code has been tested in a few regimes, I have access to very few explicit calculations to compare against. Please check that results are sensible and contact Aaron Robotham if you find any clear problems (please give me the exact inputs and explain quantitatively what the answer should be, please do not simply say *this doesn't look right*). Intuitively the comoving distance calculation makes the most sense, the others can be less obvious because of the non-standard observational frame we are considering. Depending on your precise problem you may want to construct your own effective redshift, remembering that LumDist = CoDist.(1+z_eff) and AngDist = CoDist/(1+z_eff) no matter how funky your cosmology is. Isn't it handy that we directly onserve the relative explansion factors of the Universe and not something trivial like distance!

The object at z1 should be thought of as the observer, and the object at z2 the observed. For comoving distance this makes no difference, but for zeff it tells us the redshift the object z1 is observing the object at z2 for the age of the Universe we observe z1 to be at (e.g. for a 737 cosmology and z1=1 we would be making this measurement for a Universe age of 5.75 Gyrs). This carries on into the angular size and luminosity distances, which are calculated using zeff. This is useful if you (for instance) wanted to now how much irradiance a source at z1 is experiencing due to (e.g.) a quasar at z2 at the moment we are observing z1 to be in. E.g. if z1 is observed to be quenched then this is a quantity we might want to calculate.

This!

#### R Code

Extracts of the R code run on the server to make some of the simpler calculations.

#### Reference Sets

 Planck 18 Planck Collaboration, 2018, arXiv, 1807.06209 Planck 15 Planck Collaboration, 2015, arXiv, 1502.01589 Planck 13 Planck Collaboration, 2014, A&A, 571, 16 (arXiv:1303.5076) WMAP9 Hinshaw G., et al., 2013, ApJS, 208, 19 (arXiv:1212.5226) WMAP7 Komatsu E., et al., 2010, ApJS, 192, 18 (arXiv:1001.4538) WMAP5 Komatsu E., et al., 2009, ApJS, 180, 306 (arXiv:0803.0547) WMAP3 Spergel D. N., et al., 2007, ApJS, 170, 377 (arXiv:astro-ph/0603449) WMAP1 Spergel D. N., et al., 2003, ApJS, 148, 175 (arXiv:astro-ph/0302209) Millennium Springel V., et al., 2005, Nature, 435, 629 (arXiv:astro-ph/0504097) GiggleZ Poole G. B., et al., 201?, MNRAS, ?, ? (arXiv:1407.0390)

Basic cosmological distance calculator R code used server-side to generate outputs.

Written by Aaron Robotham as part of the R celestial package (see Info tab for references).

cosdist
```function (z = 1, H0 = 100, OmegaM = 0.3, OmegaL = 1 - OmegaM -
OmegaR, OmegaR = 0, w0 = -1, wprime = 0, age = FALSE, ref,
error = FALSE)
{
.Einv=function(z, OmegaM, OmegaL, OmegaR, OmegaK, w0, wprime){
return(1/sqrt(OmegaR*(1+z)^4 + OmegaM * (1 + z)^3 + OmegaK * (1 + z)^2 + OmegaL*cosgrowRhoDE(z=z, w0=w0, wprime=wprime, rhoDE=1)))
}
.Einvz=function(z, OmegaM, OmegaL, OmegaR, OmegaK, w0, wprime){
return(1/(sqrt(OmegaR*(1+z)^4 + OmegaM * (1 + z)^3 + OmegaK * (1 + z)^2 + OmegaL*cosgrowRhoDE(z=z, w0=w0, wprime=wprime, rhoDE=1)) * (1 + z)))
}
z = as.numeric(z)
if (!all(is.finite(z))) {
stop("All z must be finite and numeric")
}
if (!all(z > -1)) {
stop("All z must be > -1")
}
if (!missing(ref)) {
params = .getcos(ref)
H0 = as.numeric(params["H0"])
OmegaM = as.numeric(params["OmegaM"])
OmegaL = as.numeric(params["OmegaL"])
if (!is.na(params["OmegaR"])) {
OmegaR = as.numeric(params["OmegaR"])
}
}
OmegaK = 1 - OmegaM - OmegaL - OmegaR
temp = function(z, H0, OmegaM, OmegaL, OmegaR, OmegaK, w0 = w0,
wprime = wprime) {
HubDist = (299792.458/H0)
temp = suppressWarnings(integral(.Einv, 0, z, OmegaM = OmegaM,
OmegaL = OmegaL, OmegaR = OmegaR, OmegaK = OmegaK,
w0 = w0, wprime = wprime))
CoDist = HubDist * temp
if (error) {
if (z > 0) {
RelError = 1e-08
}
else {
RelError = 0
}
}
if (OmegaK == 0) {
CoDistTran = CoDist
CoVol = ((4/3) * pi * CoDist^3)/1e+09
}
else {
if (OmegaK > 0) {
CoDistTran = HubDist * (1/sqrt(OmegaK)) * sinh(sqrt(OmegaK) *
CoDist/HubDist)
CoVol = ((4 * pi * HubDist^3/(2 * OmegaK)) *
((CoDistTran/HubDist) * sqrt(1 + OmegaK * (CoDistTran/HubDist)^2) -
(1/sqrt(abs(OmegaK))) * asinh(sqrt(abs(OmegaK)) *
(CoDistTran/HubDist))))/1e+09
}
if (OmegaK < 0) {
CoDistTran = HubDist * (1/sqrt(abs(OmegaK))) *
sin(sqrt(abs(OmegaK)) * CoDist/HubDist)
CoVol = ((4 * pi * HubDist^3/(2 * OmegaK)) *
((CoDistTran/HubDist) * sqrt(1 + OmegaK * (CoDistTran/HubDist)^2) -
(1/sqrt(abs(OmegaK))) * asin(sqrt(abs(OmegaK)) *
(CoDistTran/HubDist))))/1e+09
}
}
a = 1/(1 + z)
LumDist = (1 + z) * CoDistTran
AngDist = CoDistTran/(1 + z)
if (z >= 0) {
DistMod = 5 * log10(LumDist) + 25
}
else {
DistMod = NA
}
AngScale = AngDist * (pi/(180 * 60 * 60)) * 1000
if (age) {
HT = (3.08568025e+19/(H0 * 31556926))/1e+09
UniAge = HT * integral(.Einvz, 0, Inf, OmegaM = OmegaM,
OmegaL = OmegaL, OmegaR = OmegaR, OmegaK = OmegaK,
w0 = w0, wprime = wprime)
zAge = HT * integral(.Einvz, 0, z, OmegaM = OmegaM,
OmegaL = OmegaL, OmegaR = OmegaR, OmegaK = OmegaK,
w0 = w0, wprime = wprime)
}
if (error) {
if (age) {
return = c(z = z, a = a, CoDist = CoDist, LumDist = LumDist,
AngDist = AngDist, CoDistTran = CoDistTran,
DistMod = DistMod, AngScale = AngScale, CoVol = CoVol,
HubTime = HT, UniAgeNow = UniAge, UniAgeAtz = UniAge -
zAge, TravelTime = zAge, RelError = RelError)
}
else {
return = c(z = z, a = a, CoDist = CoDist, LumDist = LumDist,
AngDist = AngDist, CoDistTran = CoDistTran,
DistMod = DistMod, AngScale = AngScale, CoVol = CoVol,
RelError = RelError)
}
}
else {
if (age) {
return = c(z = z, a = a, CoDist = CoDist, LumDist = LumDist,
AngDist = AngDist, CoDistTran = CoDistTran,
DistMod = DistMod, AngScale = AngScale, CoVol = CoVol,
HubTime = HT, UniAgeNow = UniAge, UniAgeAtz = UniAge -
zAge, TravelTime = zAge)
}
else {
return = c(z = z, a = a, CoDist = CoDist, LumDist = LumDist,
AngDist = AngDist, CoDistTran = CoDistTran,
DistMod = DistMod, AngScale = AngScale, CoVol = CoVol)
}
}
}
return(as.data.frame(t(Vectorize(temp)(z = z, H0 = H0, OmegaM = OmegaM,
OmegaL = OmegaL, OmegaR = OmegaR, OmegaK = OmegaK, w0 = w0,
wprime = wprime))))
}
}```
cosgrow
```function (z = 1, H0 = 100, OmegaM = 0.3, OmegaL = 1 - OmegaM -
OmegaR, OmegaR = 0, w0 = -1, wprime = 0, Sigma8 = 0.8, fSigma8 = FALSE,
Dist = "Co", Mass = "Msun", ref)
{
.Einv=function(z, OmegaM, OmegaL, OmegaR, OmegaK, w0, wprime){
1/sqrt(OmegaR*(1+z)^4 + OmegaM * (1 + z)^3 + OmegaK * (1 + z)^2 + OmegaL*cosgrowRhoDE(z=z, w0=w0, wprime=wprime, rhoDE=1))
}
.Einva3=function(a, OmegaM, OmegaL, OmegaR, OmegaK, w0, wprime){
1/(a^3*(sqrt(OmegaR*a^(-4) + OmegaM*a^(-3) + OmegaK*a^(-2) + OmegaL*cosgrowRhoDE(z=1/a-1, w0=w0, wprime=wprime, rhoDE=1)))^3)
}
z = as.numeric(z)
if (!Dist %in% c("Co", "Ang", "m")) {
stop("Dist must be one of Co, Ang or m")
}
if (!Mass %in% c("Msun", "kg")) {
stop("Mass must be one of Msun or kg")
}
if (!all(is.finite(z))) {
stop("All z must be finite and numeric")
}
if (!all(z > -1)) {
stop("All z must be > -1")
}
if (!missing(ref)) {
params = .getcos(ref)
H0 = as.numeric(params["H0"])
OmegaM = as.numeric(params["OmegaM"])
OmegaL = as.numeric(params["OmegaL"])
if (!is.na(params["OmegaR"])) {
OmegaR = as.numeric(params["OmegaR"])
}
if (!is.na(params["Sigma8"])) {
Sigma8 = as.numeric(params["Sigma8"])
}
}
temp = function(z, H0, OmegaM, OmegaL, OmegaR, w0, wprime) {
OmegaK = 1 - OmegaM - OmegaL - OmegaR
OmegaSum = OmegaR * (1 + z)^4 + OmegaM * (1 + z)^3 +
OmegaK * (1 + z)^2 + OmegaL * cosgrowRhoDE(z = z,
w0 = w0, wprime = wprime, rhoDE = 1)
Hz = H0 * sqrt(OmegaSum)
CoVel = 299792.458 * integral(.Einv, 0, z, OmegaM = OmegaM,
OmegaL = OmegaL, OmegaR = OmegaR, OmegaK = OmegaK,
w0 = w0, wprime = wprime)
OmegaRatz = (OmegaR * (1 + z)^4)/OmegaSum
OmegaMatz = (OmegaM * (1 + z)^3)/OmegaSum
OmegaLatz = OmegaL * cosgrowRhoDE(z = z, w0 = w0, wprime = wprime,
rhoDE = 1)/OmegaSum
OmegaKatz = (OmegaK * (1 + z)^2)/OmegaSum
Factor = (5 * OmegaM/2) * (Hz/H0) * (1 + z) * integral(.Einva3,
0, 1/(1 + z), OmegaM = OmegaM, OmegaL = OmegaL, OmegaR = OmegaR,
OmegaK = OmegaK, w0 = w0, wprime = wprime)
Factor0 = (5 * OmegaM/2) * integral(.Einva3, 0, 1, OmegaM = OmegaM,
OmegaL = OmegaL, OmegaR = OmegaR, OmegaK = OmegaK,
w0 = w0, wprime = wprime)
Sigma8atz = Sigma8 * (Factor/Factor0)/(1 + z)
if (fSigma8 == FALSE) {
Rate = -1 - OmegaMatz/2 + OmegaLatz + (5 * OmegaMatz)/(2 *
Factor)
}
else {
Rate = Sigma8atz * (-1 - OmegaMatz/2 + OmegaLatz +
(5 * OmegaMatz)/(2 * Factor))
}
G = 6.67384e-11
km2m = 1000
Mpc2m = 3.08567758e+22
Msol2kg = 1.9891e+30
RhoCrit = (3 * Hz^2)/(8 * pi * G) * (km2m^2)/Mpc2m^2
if (Mass == "Msun") {
RhoCrit = RhoCrit/Msol2kg
}
if (Dist == "Ang") {
RhoCrit = RhoCrit * Mpc2m^3
}
if (Dist == "Co") {
RhoCrit = RhoCrit * Mpc2m^3/(1 + z)^3
}
RhoMean = RhoCrit * OmegaMatz
Decelq = OmegaMatz/2 + OmegaRatz - OmegaLatz
return = c(z = z, a = 1/(1 + z), H = Hz, CoVel = CoVel,
OmegaM = OmegaMatz, OmegaL = OmegaLatz, OmegaR = OmegaRatz,
OmegaK = OmegaKatz, Decelq = Decelq, Factor = Factor,
Rate = Rate, Sigma8 = Sigma8atz, RhoCrit = RhoCrit,
RhoMean = RhoMean)
}
return(as.data.frame(t(Vectorize(temp)(z = z, H0 = H0, OmegaM = OmegaM,
OmegaL = OmegaL, OmegaR = OmegaR, w0 = w0, wprime = wprime))))
}
}```

ICRAR 2015, written by Aaron Robotham and Joseph Dunne