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 [email protected]
For clarity this website uses the convention of a prefix c, p or l for comoving, physical (or proper) and luminosity distances respectively.
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).
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))))
}
}