A glitch in Jorrie’s Cosmo-Calculator?

[pbuk]

I noticed that using version (1.2.0), I can’t duplicate your D_now(z) in one run of the calculator; I had to calculate one or two values at a time.

Yes, it doesn't currently do linear intervals; this should be easy to add.

The D_v1.2.0(z) result shows a big improvement over the Jorrie 2022 result, but it still deviates systematically from D_cosmic(z).

I'm not really surprised at this deviation, steps to improve it include:

• Check the precision of the physical constants and model values
• Try midpoint instead of trapezium rule which always overestimates for a concave function.
• Add error estimation/adaptive step size.

But these are all enhancements for 1.4.0.

 

[pbuk]

I'm not really surprised at this deviation, steps to improve it include:

These calculators use different values for the constants* (and Cosmic doesn't take $\Omega_R$ into account at all). It's not surprising they give different results with a regular pattern. Note that the "new Jorrie" site does not yet allow overriding the default values (again this is something that is easy to add but pointless until we are happy with ALL the calculations, not just $D_{now}$.

*e.g. $H_0$ Jorrie, pbuk 67.74 km/s/Mpc, cosmic 67.04, Gnedin 68.14.

 

[JimJCW]

These calculators use different values for the constants* (and Cosmic doesn't take $\Omega_R$ into account at all). It's not surprising they give different results with a regular pattern. Note that the "new Jorrie" site does not yet allow overriding the default values (again this is something that is easy to add but pointless until we are happy with ALL the calculations, not just $D_{now}$.

*e.g. $H_0$ Jorrie, pbuk 67.74 km/s/Mpc, cosmic 67.04, Gnedin 68.14.

You are right; all other calculators I have tried do not include Ω_R. That’s why I like and use Jorrie’s calculator. The effect of Ω_R is significant only during early cosmological times (high z values) and it probably does not have any noticeable effects in our D_now calculations.

All calculators I tried use H₀, Ω_m and Ω_Λ as input parameters and I use H₀ = 67.74 km/s/Mpc, Ω_m = 0.309 and Ω_Λ = 0.691 for all of them.

 

[pbuk]

You are right; all other calculators I have tried do not include Ω_R. That’s why I like and use Jorrie’s calculator. The effect of Ω_R is significant only during early cosmological times (high z values) and it probably does not have any noticeable effects in our D_now calculations.

It would be easy enough to exclude it for better comparisons: then with the same constants the only thing that should be different is the integration algorithm.

All calculators I tried use H₀, Ω_m and Ω_Λ as input parameters and I use H₀ = 67.74 km/s/Mpc, Ω_m = 0.309 and Ω_Λ = 0.691 for all of them.

Ah OK. Is Ω_m = 0.309 correct or is it Ω_m = 0.3089 as per Jorrie's defaults?

 

[Jorrie]

Ah OK. Is Ω_m = 0.309 correct or is it Ω_m = 0.3089 as per Jorrie's defaults?

Keep in mind that the accuracy of our knowledge of the Omega's is far worse than the rounding contributions. It is good to chase mathematical integrity in algorithms, but is it worth it to get far ahead of the quality of the observations that we work with?

 

[pbuk]

is it worth it to get far ahead of the quality of the observations that we work with?

Not in actual use, but in testing it is important to determine whether differences in results arise from different models or poor integration methods.

 

[JimJCW]

Ah OK. Is Ω_m = 0.309 correct or is it Ω_m = 0.3089 as per Jorrie's defaults?

See discussion in StackExchange:

If we change z_eq from 3370 to 999998 in the current version of the calculator (i.e., neglect the radiation energy density), Ω_m,0 will change from 0.3089 to 0.309.

In the current version, changing z_eq in the calculator does not change the calculated D_now. When it is fixed, we can compare the calculated D_now using large z_eq values with that calculated with the cosmic calculator (no radiation energy term).

 

[JimJCW]

Hi @Ibix,

In the figures included in Post #70, D_v1.2.0(z) by @pbuk shows a small, systematic deviation from D_cosmic(z) near z = 0. It is interesting to see whether this is caused by the omission of the radiation term Ω_R,0 in the cosmic model. Your python version includes the radiation term and may provide information concerning this. Would you please post your result with more digits? (See Post #57)

z​	Gnedin​	Jorrie, 2022​	ibix​	ibix_2​	cosmic​	Toth​
0.1	1.41004	1.4084	1.41		1.409592	1.409608
0.09	1.27214	1.270609	1.27		1.27174	1.271756
0.08	1.13354	1.132106	1.13		1.133189	1.133219
0.07	0.99426	0.992931	0.994		0.993944	0.993964
0.06	0.85427	0.85305	0.854		0.854005	0.854024
0.05	0.7136	0.712477	0.713		0.713377	0.713399
0.04	0.57225	0.571225	0.572		0.572063	0.57209
0.03	0.43023	0.429278	0.43		0.430065	0.430062
0.02	0.28748	0.286651	0.28		0.287385	0.287382
0.01	0.14408	0.14336	0.144		0.144029	0.144017
0	0	0.000621	0		0	0
-0.01	0.14474	0.14526	0.145		0.144698	0.144702
-0.02	0.29015	0.290567	0.29		0.290061	0.290057
-0.03	0.43623	0.436536	0.436		0.436084	0.436097
-0.04	0.58295	0.583161	0.583		0.582763	0.582789
-0.05	0.73032	0.730433	0.73		0.730094	0.730101
-0.06	0.87834	0.878358	0.878		0.878072	0.878098
-0.07	1.02699	1.026926	1.03		1.026688	1.026715
-0.08	1.17631	1.176124	1.18		1.175941	1.175951
-0.09	1.32623	1.325968	1.33		1.325824	1.32584
-0.1	1.47679	1.476416	1.48		1.476332	1.476349

 

[pbuk]

Hi @JimJCW,

Yes, this is part of my plan for generating test data. I am away on a sailing trip at the moment so can't work on it for a while though.

 

[JimJCW]

Hi @Jorrie, @pbuk, @lbix,

I think the cosmology calculator by The International Centre for Radio Astronomy Research, ICRAR, can be useful for testing in our calculator development. It is written in R language and can have Ω_R as an input parameter. The inputs to this calculator for PLANCK 2015 Data are:

H₀ = 67.74​

Ω_Λ = 0.691​

Ω_R = 0.0000916 [obtained from Ω_R = Ω_m / (z_eq + 1)]​

Ω_m = 1 - 0.691 - 0.0000916 = 0.3089​

I have some questions about some of the terms:

What is the difference between ‘The Comoving Radial Distance LoS to z’ and ‘The Comoving Radial Distance Tran to z’?​

​

What is Sigma8[0]?​

 

[Jorrie]

Hi @Jorrie, @pbuk, @lbix,

I think the cosmology calculator by The International Centre for Radio Astronomy Research, ICRAR, can be useful for testing in our calculator development. It is written in R language and can have Ω_R as an input parameter. ...

I have some questions about some of the terms:

What is the difference between ‘The Comoving Radial Distance LoS to z’ and ‘The Comoving Radial Distance Tran to z’?​

​

What is Sigma8[0]?​

You will find the definitions of comoving distances in Hogg
(arXiv:astro-ph/9905116v4 16 Dec 2000) section 5.

I have no clue about Sigma8[0]

 

[pbuk]

What is Sigma8[0]?

Apparently it "measures the amplitude of the (linear) power spectrum on the scale of 8 h-1 Mpc". I have no idea what any of this means (I'm a mathematician and computer scientist, not a cosmologist).

 

[pbuk]

In the figures included in Post #70, D_v1.2.0(z) by @pbuk shows a small, systematic deviation from D_cosmic(z) near z = 0. It is interesting to see whether this is caused by the omission of the radiation term Ω_R,0 in the cosmic model. Your python version includes the radiation term and may provide information concerning this. Would you please post your result with more digits? (See Post #57)

I have posted results from Ibix's Python program below with a LOT more digits. I have also written a much better integrator for my Javascript code* which I am pleased to say agrees almost to machine precision with Python's scipy.integration.quad here. I have posted these results (with differing digits shown in bold) as "pbuk g7k15" below. I may be able to get this into GitHub and the demo site over the weekend.

I have also run Ibix's code with the radiation parameters $\Omega_{R0}$ and $\Omega_R$ removed. You can see from the table below that it doesn't make much difference so there must be some other reason for the differences from the other calculators.

z​	Gnedin​	Jorrie, 2022​	ibix​	pbuk g7k15​	ibix no $\Omega_R$​	cosmic​	Toth​
0.1	1.41004	1.4084	1.409699537167435	1.409699537167435	1.409707	1.409592	1.409608
0.09	1.27214	1.270609	1.2718364599617440	1.2718364599617438	1.271842	1.27174	1.271756
0.08	1.13354	1.132106	1.133275788110779	1.133275788110779	1.133280	1.133189	1.133219
0.07	0.99426	0.992931	0.9940198604891496	0.9940198604891496	0.994023	0.993944	0.993964
0.06	0.85427	0.85305	0.8540712154740808	0.8540712154740807	0.854074	0.854005	0.854024
0.05	0.7136	0.712477	0.7134325928253945	0.7134325928253945	0.713434	0.713377	0.713399
0.04	0.57225	0.571225	0.5721069353785664	0.5721069353785664	0.572108	0.572063	0.57209
0.03	0.43023	0.429278	0.4300973905418238	0.4300973905418237	0.430098	0.430065	0.430062
0.02	0.28748	0.286651	0.2874073115884052	0.2874073115884052	0.287408	0.287385	0.287382
0.01	0.14408	0.14336	0.1440402587352972	0.1440402587352972	0.144040	0.144029	0.144017
0	0	0.000621	0	0	0	0	0
-0.01	0.14474	0.14526	0.1447094881728547	0.1447094881728547	0.144709	0.144698	0.144702
-0.02	0.29015	0.290567	0.2900840205228997	0.2900840205228997	0.290084	0.290061	0.290057
-0.03	0.43623	0.436536	0.43611920283095362	0.43611920283095357	0.436119	0.436084	0.436097
-0.04	0.58295	0.583161	0.5828104320285896	0.5828104320285897	0.582809	0.582763	0.582789
-0.05	0.73032	0.730433	0.7301528965670417	0.7301528965670419	0.730151	0.730094	0.730101
-0.06	0.87834	0.878358	0.8781415770516691	0.8781415770516691	0.878139	0.878072	0.878098
-0.07	1.02699	1.026926	1.0267712471476862	1.0267712471476862	1.026768	1.026688	1.026715
-0.08	1.17631	1.176124	1.1760364747623280	1.1760364747623278	1.176033	1.175941	1.175951
-0.09	1.32623	1.325968	1.3259316235080327	1.3259316235080325	1.325927	1.325824	1.32584
-0.1	1.47679	1.476416	1.4764508544505859	1.4764508544505859	1.476445	1.476332	1.476349

* Given that scipy.integrate.quad simply links the venerable Fortran library QUADPACK in, I am quite pleased that my Gauss-Kronrod [7, 15] implementation seems to perform as well, on this problem at least. I need to make sure it handles the discontinuity at $s = 0$ well and add transformation for the integration to $\infty$ before I put it into the 'new' LightCone7.

 

[JimJCW]

Hi @pbuk,

Thanks for the results (Post #83) obtained with the python program by @Ibix. It is nice to know that when the Ω_R term is included, pbuk g7k15 = version (1.2.0) = ibix. Please confirm that if Ω_R = 0, pbuk g7k15 = ibix is also true.

When comparing the various results, what do you think about using the following assumptions?

The cosmic calculator works properly when Ω_R = 0.​

​

The ICRAR calculator works properly when the Ω_R term is included.​

 

[Jorrie]

I need to make sure it handles the discontinuity at $s = 0$ well and add transformation for the integration to $\infty$ before I put it into the 'new' LightCone7.

I presume you meant "at $z = 0$" ?
There is no real discontinuity at z = 0. One could just as well have had distances (and some other outputs) go negative for future times (negative z). The reason the original Lightcone team have chosen to keep things positive was for better resolution on the vertical graph axis, while keeping the graphs vertically compact.

Maybe we should start talking about Lightcone8 (or X) due to the major rework and so avoid confusion...?

 

[pbuk]

Please confirm that if Ω_R = 0, pbuk g7k15 = ibix is also true.

Yes this is true: the calculations in pbuk g7k15, ibix and jorrie are in fact identical except for the numerical integration method summarised below

• jorrie uses a simple fixed step method with the 'glitch' in the title of this thread.
• ibix uses Python's scipy.integrate.quad which links the ancient, well-tested FORTRAN 77 code QUADPACK.
• pbuk g7k15 (note this is not yet the version used on the website) uses a new JavaScript routine written from scratch. I believe that for infinite ranges this uses the same transformation and underlying algorithm as QUADPACK (15-point Gauss-Kronrod), although there are differences in the error estimation and hence adaptive step calculations, however for finite ranges QUADPACK uses a higher order (21 point) method: this does not seem to make any material difference.

When comparing the various results, what do you think about using the following assumptions?

The cosmic calculator works properly when Ω_R = 0.​

​

The ICRAR calculator works properly when the Ω_R term is included.​

I'm not sure, other calculators seem to include other parameters, for instance Ned Wright's includes this:

// by Ned Wright
// 25 Jul 1999 - revised 2 Dec 2005
// Copyright Edward L. Wright, all rights reserved.

// ...

// correction for annihilations of particles not present now like e+/e-
// added 13-Aug-03 based on T_vs_t.f
  var lpz = Math.log((1+1.0*z))/Math.log(10.0);
  var dzage = 0;
  if (lpz >  7.500) dzage = 0.002 * (lpz -  7.500);
  if (lpz >  8.000) dzage = 0.014 * (lpz -  8.000) +  0.001;
  if (lpz >  8.500) dzage = 0.040 * (lpz -  8.500) +  0.008;
  if (lpz >  9.000) dzage = 0.020 * (lpz -  9.000) +  0.028;
  if (lpz >  9.500) dzage = 0.019 * (lpz -  9.500) +  0.039;
  if (lpz > 10.000) dzage = 0.048;
  if (lpz > 10.775) dzage = 0.035 * (lpz - 10.775) +  0.048;
  if (lpz > 11.851) dzage = 0.069 * (lpz - 11.851) +  0.086;
  if (lpz > 12.258) dzage = 0.461 * (lpz - 12.258) +  0.114;
  if (lpz > 12.382) dzage = 0.024 * (lpz - 12.382) +  0.171;
  if (lpz > 13.055) dzage = 0.013 * (lpz - 13.055) +  0.188;
  if (lpz > 14.081) dzage = 0.013 * (lpz - 14.081) +  0.201;
  if (lpz > 15.107) dzage = 0.214;
  zage = zage*Math.pow(10.0,dzage);

There is clearly a lot more to this than I realized at first: I'm just a humble coder

 

[pbuk]

I presume you meant "at $z = 0$" ?

Well I meant at $s = 0$ because I was confusing the calculation of $T_{now} = \int_{s=S}^\infty \frac{1}{sH(s)}$ with $D_{hor} = \frac1S \int_{s=0}^S \frac{1}{H(s)}$ and thinking we had to compute $\int_{s=0}^S \frac{1}{sH(s)}$ somewhere! So there is no discontinuity: that's good, my integrator doesn't handle discontinuities as well as QUADPACK.

Maybe we should start talking about Lightcone8 (or X) due to the major rework and so avoid confusion...?

Yes, a bit of housekeeping to do here.

 

[JimJCW]

It may be interesting to compare D_now vs z curves from different calculators for larger z values. The result is shown below:

This suggests that for large z values, the results from all these calculators are almost indistinguishable and the glitch discussed in this thread is very small. I am working on the comparison for small z values.

Note that D_Jorrie,2021 is the modified version discussed in Posts #15 and #17.

 

[JimJCW]

In the following discussion, I assume that the ICRAR calculator works properly and use it as a reference for comparison.

These figures suggest that there exist real discrepancies among the three calculators. When we modify and/or rewrite Jorrie’s calculator, maybe we can use ICRAR calculator as a tool for testing.

 

[pbuk]

When we modify and/or rewrite Jorrie’s calculator, maybe we can use ICRAR calculator as a tool for testing.

It's already rewritten, now working here: https://light-cone-calc.github.io. There's not much point testing it against the ICRAR calculator though as that uses a more sophisticated model.

 

[JimJCW]

It's already rewritten, now working here: https://light-cone-calc.github.io. There's not much point testing it against the ICRAR calculator though as that uses a more sophisticated model.

Great! I can make more comparisons.

Please explain what you mean by “. . . the ICRAR calculator . . . uses a more sophisticated model”. Doesn’t it calculate the same D_now?

 

[pbuk]

Please explain what you mean by “. . . the ICRAR calculator . . . uses a more sophisticated model”. Doesn’t it calculate the same D_now?

Not exactly, there is a much more sophisticated treatment of dark energy. You can see the code (in R) if you go to https://cosmocalc.icrar.org/ and click the "R code" tab.

 

[JimJCW]

Not exactly, there is a much more sophisticated treatment of dark energy. You can see the code (in R) if you go to https://cosmocalc.icrar.org/ and click the "R code" tab.

How about comparing the cosmic model with your new version (setting Ω_R = 0)? Similar discrepancy exists.

 

[pbuk]

How about comparing the cosmic model with your new version (setting Ω_R = 0)? Similar discrepancy exists.

Rather than reverse engineer output I'd rather look at the implementations themselves. For instance the constants used in Jorrie's original model include these:

// https://github.com/light-cone-calc/light-cone-calc/blob/main/src/model.ts
const physicalConstants = {
  rhoConst: 1.7885e9, // 3 / (8 pi G)
  secInGy: 3.1536e16, // s / Gyr
  tempNow: 2.725, // CMB temperature now
  convertToGyr: 1 / 978, // Convert km/s/Mpc -> Gyr^-1
};

whereas cosmic uses these

// https://github.com/joshkempner/Cosmic.NET/blob/main/Cosmic.NET/cosmo/Cosmo.cs
        private const double c = 2.99792458e5; // speed of light
        private const double G = 6.67259e-8;   // gravitational constant
        private const double KmPerMpc = 3.08567758e19;
        private const double TropicalYear = 3.1556926e7; // in seconds

Also, as noted in the code, I'm not sure about the calculation of $\Omega_m$

//  s_eq: 1 + 3370, // Stretch when OmegaM=OmegaR
...

  //@TODO check this - should it be s_eq + 1 as the original, or as below
  // from Ibix?
  // const OmegaM = ((Omega - OmegaL) * s_eq) / (s_eq + 1); // Energy density of matter
  // const OmegaR = OmegaM / s_eq; // Energy density of radiation
  const OmegaM = ((Omega - OmegaL) * (s_eq + 1)) / (s_eq + 2); // Energy density of matter

A better place to discuss this is probably https://github.com/light-cone-calc/light-cone-calc/issues - in fact I've already noted the last point here https://github.com/light-cone-calc/light-cone-calc/issues/6.

 

[Jorrie]

Pbuk, Jim and Ibix, thanks for your great work!
We have a local 4-day weekend starting tomorrow, so I will take a look good at your Lightcone8 version and also at the issue that you raised on Github and discuss it there.

 

[Jorrie]

Yes, There is an error in my old Legacy js code, line 356
Ibix is correct, i.e.
// const OmegaM = ((Omega - OmegaL) * s_eq) / (s_eq + 1); // Energy density of matter

Tiny, tiny influence, considering the uncertainty on the z_eq value.

@pbuk, which one did you use in Lightcone8?

 

[pbuk]

Ibix is correct, i.e.
// const OmegaM = ((Omega - OmegaL) * s_eq) / (s_eq + 1); // Energy density of matter

Tiny, tiny influence, considering the uncertainty on the z_eq value.

Yes of course.

@pbuk, which one did you use in Lightcone8?

I used yours, but it's a quick fix.

 

[pbuk]

Tiny, tiny influence, considering the uncertainty on the z_eq value.

Tiny but important: this was the remaining part of the problem causing the offset at z = 0.

z: 0
Dnow: 3.552713678800501e-15
Dthen: 3.552713678800501e-15
Vnow: 2.460744627831758e-16
Vthen: 2.4607446278317573e-16

These values should all be 0 of course. Now fixed.

 

[Jorrie]

z: 0
Dnow: 3.552713678800501e-15
Dthen: 3.552713678800501e-15
Vnow: 2.460744627831758e-16
Vthen: 2.4607446278317573e-16

These values should all be 0 of course. Now fixed.

Great stuff! I will link https://light-cone-calc.github.io/ in my signature and you can do the same, if you so wish.
Maybe do a last edit to take some credit for the improvement?

 

[JimJCW]

Hi @Jorrie, @pbuk,

When trying LightCone8, I noticed that if z_lower = 0, a duplicated line is printed at z = 0. For example:

 

[pbuk]

Hi @Jorrie, @pbuk,

When trying LightCone8, I noticed that if z_lower = 0, a duplicated line is printed at z = 0.

Thanks for that, it is an unintended side effect of fixing the range so it includes z = 0 instead of z = 1e-6 for instance. I'll fix the fix.

 

[JimJCW]

Hi @Jorrie, @pbuk,

Let’s review the current situation of Jorrie’s calculator. First, LightCone8 updated by @pbuk has eliminated the glitch discussed in this thread (see Post #1):

ICRAR is an online calculator that includes Ω_R,0 as an input parameter, just like LightCone8 (z_eq = 3370). A comparison of the calculated D_now’s using the two models is shown below:

The two results are almost indistinguishable. However, when their difference is plotted against z, the result shows,

The fact that ‘ICRAR - LightCone8’ is not a horizontal line near y = 0 suggests that there exists discrepancy between the two calculators. To look into this further, we note that ICRAR (Ω_R=0), LightCone8 (Ω_R=0), and cosmic (no Ω_R) should give the same results, but they don’t:

The figure suggests that when Ω_R=0, ICRAR and cosmic are consistent with each other, but not with LightCone8. The source of this discrepancy is unknown. It would be nice, though, to be able to identify it.

To summarize:

LightCone8 has removed the glitch discussed in this thread. However, it is noticed that there exist discrepancies between LightCone8 and some other online calculators, such as ICRAR and cosmic.​

 

[pbuk]

The figure suggests that when Ω_R=0, ICRAR and cosmic are consistent with each other, but not with LightCone8. The source of this discrepancy is unknown. It would be nice, though, to be able to identify it.

My first suspicion is the values of certain constants, particularly the conversion factor from $km\\s^{-1}\\Mpsc^{-1}$ to $Gyr^{-1}$ which is a simple $\frac 1 {978}$ in LightCone8 instead of something like $\frac{365.2421897 \times 86400 \times 1000000 \pi}{96939420213600000}$.

I've raised an issue and will look at it: https://github.com/light-cone-calc/light-cone-calc/issues/11

 

[Jorrie]

I have commented on the Github issues page.

 

[JimJCW]

Hi @Jorrie and @pbuk,

I noticed some discrepancies in results from LightCone8. Let’s use the following input values:

Equations in Tutorial, part III,

Ω_m = (Ω - Ω_Λ) S_eq / (1 + S_eq), where S_eq = z_eq + 1​

Ω_R = Ω_m / S_eq​

give,

Ω_m = 0.308908363
Ω_R = 0.000091637

These numbers are compared with output of LightCone8 below:

​	z	OmegaM​	OmegaL​	OmegaR​	OmegaT​
LightCone8​	0​	0.309000000​	6.910000000	0.00009169139​	1.000091691​
From Equations​	0​	0.308908363​	6.91​	0.000091637​	1​

Note that under Conversions in the output, LightCone8 gives,

Conclusion: The values of Ω_m and Ω_T calculated with LightCone8 are questionable.