I A glitch in Jorrie’s Cosmo-Calculator?

[pbuk]

Any thoughts?

Yes that's exactly what I was thinking (plus step 0 being sanitizing inputs and other initialization).

 

[Jorrie]

I'm impressed with the apparent enthusiasm for the community project, which should finally become a PF Cosmo tool, decoupled from anyone member and managed by an elected committee on a neutral platform.

 

[pbuk]

Well kind of self-elected! PM GitHub user ids to be added as project admins.

 

[PAllen]

If there is no good open source quadrature package available for JavaScript, it seems to me that the methods described in section 4.5 of the third edition of “Numerical Recipes” are more than adequate for the required integrals, and are much less code than trying to reproduce something like QUADPACK.

 

[robphy]

I'm quite rusty on my cosmology... so I can't vouch for the equations.
When this thread started, I was interested in trying to reproduce the calculations in Desmos.

So, here is what I have so far:
https://www.desmos.com/calculator/dorewco5ms
See the table at the end to check the numerical calculations.

You'll have to check the equations ( based on calculations by Jorrie and Ibix ; see also work by pbuk ) all documented in the Desmos file.

I don't have any more time to play around with it.
Maybe someone might find it useful (once the equations are verified).

(A note on Desmos... the ordering of the entries doesn't matter.
Unlike geogebra or trinket/Glowscript,
every Desmos "save" generates a new unique URL... which the saver has to share for others to see.)

 

[PAllen]

@robphy , looking at Desmos worksheet (really nice, by the way, never heard of this tool), I have one question. You have a formula using $f_T$, but I don’t see a definition for it.

 

[robphy]

@robphy , looking at Desmos worksheet (really nice, by the way, never heard of this tool), I have one question. You have a formula using $f_T$, but I don’t see a definition for it.

It’s in the “constants” folder.
It’s a conversion unit.

 

[JimJCW]

I'm not. Jorrie's calculation uses a fixed step size, the integration steps don't match the result steps excactly, it doesn't treat the integration to infinity properly and it only manages to get as close as it does by using far more integration steps than it needs to - in short it needs rewriting.

I have been calculating D_now using various calculators. The results are often slightly different from one another. The results calculated with cosmic and Toth calculators are added to the table shown in Post #57. They are consistent with that calculated by @Ibix but with more digits kept:

z​	Gnedin​	Jorrie, 2022​	ibix​	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.287	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

For the discussions blow, let’s assume that the cosmic result is correct and use it as a reference for comparison.

For two properly working calculators, the difference in their calculated D_now’s as a function of z should be a slightly fluctuating plot around the y = 0 line, as exemplified by the D_cosmic - D_Toth plot shown in the following figure:

The D_cosmic - D_Jorrie plot, however, shows a peculiar behavior: there is an abrupt break around z = 0 and D_Jorrie(z) deviates systematically from D_cosmic(z). It will be nice if we can get a clue of why this is happening. If you or @Jorrie have some ideas about the location(s) in the Calculation.js file that is responsible for this behavior, please let us know.

 

[pbuk]

If you or @Jorrie have some ideas about the location(s) in the Calculation.js file that is responsible for this behavior, please let us know.

It's more complicated than that, see post #55 above. Anyway we have moved on from there, the Dnow calculations seem to be much improved over the range near z = 0, but other calculations need to be checked and a few other things sorted. The results in the table below are taken from the current development version (1.2.0) running at https://lightcone7.github.io/LightCone7-develop.html.

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

 

[JimJCW]

It's more complicated than that, see post #55 above. Anyway we have moved on from there, the Dnow calculations seem to be much improved over the range near z = 0, but other calculations need to be checked and a few other things sorted. The results in the table below are taken from the current development version (1.2.0) running at https://lightcone7.github.io/LightCone7-develop.html.

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. Some of the D_now values probably were not correctly entered in your table. For example, D_now = 0 for z = 0.03.

If we add (D_cosmic - D_v1.2.0) vs z curve to the figure shown in Post #68, we get the following result:

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). A better indication of this is shown below:

 

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