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

**z**
- Redshift, where z must be > -1

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

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

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

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

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

**w[0]**
- 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[0] and w'. The dark energy EoS evolution is parameterised such that P=w.rho.c^2 and w=w[0]+2.w'.(1-1/(1+z)). This means rhoDE=rhoDE[0].e^{-6.w'.(1-1/(1 + z))}/(1 + z)^{-(3 + 3.w_0 + 6.w')}, as per Wright (2006).

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[0]' into the OmegaL[0] field to set OmegaL[0] to
1-OmegaM[0]
for all calculations.
The
**Reference Set**
box lets you set the variables (H[0], OmegaM[0], OmegaL[0], OmegaR[0] and Sigma8[0]) 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.

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[0]' into the OmegaL[0] field to set OmegaL[0] to
1-OmegaM[0]
for all calculations.
The
**Reference Set**
box lets you set the variables (H[0], OmegaM[0], OmegaL[0], OmegaR[0] and Sigma8[0]) 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.**

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[0], OmegaM[0], OmegaL[0] and OmegaR[0]) 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!

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!

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

Driver S.P. & Robotham A.S.G., 2010, MNRAS, 407, 2131

Hamilton A.J.S., 2001, MNRAS, 322, 419

Hogg D.W., et al., 1999 (arXiv 9905116)

Lahav O., et al., 1991, MNRAS, 251, 136

Liske J., 2000, MNRAS, 319, 557L

Peacock J.A., 1999, Cosmological Physics, Cambridge University Press

Wright E.L., 2006, PASP, 118, 1711

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).

function(z = 1, H0 = 100, OmegaM = 0.3, OmegaL = 1 - OmegaM - OmegaR, OmegaR = 0, age = FALSE, ref, error = FALSE) { 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 Einv = function(z, OmegaM, OmegaL, OmegaR, OmegaK) { 1/sqrt(OmegaR * (1 + z)^4 + OmegaM * (1 + z)^3 + OmegaK * (1 + z)^2 + OmegaL) } if (age) { Einvz = function(z, OmegaM, OmegaL, OmegaR, OmegaK) { 1/(sqrt(OmegaR * (1 + z)^4 + OmegaM * (1 + z)^3 + OmegaK * (1 + z)^2 + OmegaL) * (1 + z)) } } temp = function(z, H0, OmegaM, OmegaL, OmegaR, OmegaK) { HubDist = (299792.458/H0) temp = suppressWarnings(integrate(Einv, 0, z, OmegaM = OmegaM, OmegaL = OmegaL, OmegaR = OmegaR, OmegaK = OmegaK, subdivisions = 1000L)) CoDist = HubDist * temp$value if (error) { if (z > 0) { RelError = abs(temp$abs.error/temp$value) } 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 } AngSize = AngDist * (pi/(180 * 60 * 60)) * 1000 if (age) { HT = (3.08568025e+19/(H0 * 31556926))/1e+09 UniAge = HT * integrate(Einvz, 0, Inf, OmegaM = OmegaM, OmegaL = OmegaL, OmegaR = OmegaR, OmegaK = OmegaK, subdivisions = 1000L)$value zAge = HT * integrate(Einvz, 0, z, OmegaM = OmegaM, OmegaL = OmegaL, OmegaR = OmegaR, OmegaK = OmegaK, subdivisions = 1000L)$value } if (error) { if (age) { return = c(z = z, a = a, CoDist = CoDist, LumDist = LumDist, AngDist = AngDist, CoDistTran = CoDistTran, DistMod = DistMod, AngSize = AngSize, 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, AngSize = AngSize, CoVol = CoVol, RelError = RelError) } } else { if (age) { return = c(z = z, a = a, CoDist = CoDist, LumDist = LumDist, AngDist = AngDist, CoDistTran = CoDistTran, DistMod = DistMod, AngSize = AngSize, 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, AngSize = AngSize, CoVol = CoVol) } } } return(as.data.frame(t(Vectorize(temp)(z = z, H0 = H0, OmegaM = OmegaM, OmegaL = OmegaL, OmegaR = OmegaR, OmegaK = OmegaK)))) } }

function(z = 1, H0 = 100, OmegaM = 0.3, OmegaL = 1 - OmegaM - OmegaR, OmegaR = 0, Sigma8 = 0.8, fSigma8 = FALSE, Dist = "Co", Mass = "Msun", ref) { 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) { # For parameters that cant use OmegaR OmegaK = 1 - OmegaM - OmegaL OmegaSum = OmegaM * (1 + z)^3 + OmegaK * (1 + z)^2 + OmegaL Hz = H0 * sqrt(OmegaSum) OmegaMAtz = (OmegaM * (1 + z)^3)/OmegaSum OmegaLAtz = OmegaL/OmegaSum OmegaKAtz = (OmegaK * (1 + z)^2)/OmegaSum Einva3 = function(a, OmegaM, OmegaL, OmegaK) { 1/(a^3 * (sqrt(OmegaM * a^(-3) + OmegaK * a^(-2) + OmegaL))^3) } Factor = (5 * OmegaM/2) * (Hz/H0) * (1 + z) * integrate(Einva3, 0, 1/(1 + z), OmegaM = OmegaM, OmegaL = OmegaL, OmegaK = OmegaK, subdivisions = 1000L)$value Factor0 = (5 * OmegaM/2) * integrate(Einva3, 0, 1, OmegaM = OmegaM, OmegaL = OmegaL, OmegaK = OmegaK, subdivisions = 1000L)$value 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)) } # For parameters that can use OmegaR OmegaK = 1 - OmegaM - OmegaL - OmegaR OmegaSum = OmegaR * (1 + z)^4 + OmegaM * (1 + z)^3 + OmegaK * (1 + z)^2 + OmegaL Hz = H0 * sqrt(OmegaSum) OmegaMAtz = (OmegaM * (1 + z)^3)/OmegaSum OmegaLAtz = OmegaL/OmegaSum OmegaRAtz = (OmegaR * (1 + z)^4)/OmegaSum OmegaKAtz = (OmegaK * (1 + z)^2)/OmegaSum G = 6.67384e-11 # m^3 kg^-1 s^-2 km2m = 1000 Mpc2m = 3.08567758e+22 Msol2kg = 1.9891e+30 # kg RhoCrit = (3 * Hz^2)/(8 * pi * G) * (km2m^2)/Mpc2m^2 #kg/m^3 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 return = c(z = z, a = 1/(1 + z), H = Hz, OmegaM = OmegaMAtz, OmegaL = OmegaLAtz, OmegaR = OmegaRAtz, OmegaK = OmegaKAtz, 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)))) } }

ICRAR
2015, written by Aaron Robotham and Joseph Dunne