I A glitch in Jorrie’s Cosmo-Calculator?

[pbuk]

This has an ODE solver: https://ccc-js.github.io/numeric2/index.html. I think that could be slightly abused to do integration.

I don't think that is a very good idea: (i) the project is even more stale than Lightcone7 (no commits for 10 years); (ii) robust methods for quadrature are very different from methods for initial value problems.

Or we could just use the approach in the existing module.

I think that would be better.

I don't think it's a bad approach and seems pretty responsive (although I'd move some stuff around to make it a bit easier to see what's going on, and you can get rid of one pass of the integrator, I think). It's just grown messy and needs a dead code pruning round.

I agree with all of that. An approach I have found successful in the past is to extract the existing code into a stand-alone module and implement continuous integration with some high level functional tests to check you don't break anything before refactoring the existing code so that you can figure out how it works/why it doesn't.

Once you have done this you can prioritize what needs to be done next (quick fix to address bugs? rewrite algorithm for better stability/accuracy/speed? Etc...).

 

[PAllen]

(ii) robust methods for quadrature are very different from methods for initial value problems.
.

While I have no comment on the referenced project, far more development of adaptive and variable step size has occurred for ODEs than quadrature. In particular, for spikey functions poorly approximated by polynomials, a given accuracy is often achieved with much less computation by casting the quadrature as an ODE, and using modern ODE methods.

 

[pbuk]

Here it is https://lightcone7.github.io/LightCone7.html with the calculations split out into a separate module (source here) with the bare minimum of changes to get it to work independently.

Needs some tidying up and documentation: PM me with github user details if you want to use/take over this.

 

[pbuk]

In particular, for spikey functions poorly approximated by polynomials, a given accuracy is often achieved with much less computation by casting the quadrature as an ODE, and using modern ODE methods.

I didn't realize that, although these functions are generally smooth.

Far more development of adaptive and variable step size has occurred for ODEs than quadrature.

Yes this is true in general, but there is not much implemented in Javascript for ODEs either.

 

[PAllen]

I didn't realize that, although these functions are generally smooth.

Agreed.

Yes this is true in general, but there is not much implemented in Javascript for ODEs either.

Probably true, but in other languages, publicly available implementations of most any method should be available. The SciPy quadrature interface @Ibix mentioned looks good. It uses adaptive step size. Scipy also has ODE methods.

[edit: I see this was all discussed before - browser availability is a requirement. I have no idea what is available open source in this context.]

 

[Jorrie]

I have some experience of collaborative Javascript projects, @Jorrie do you mind if I set something up as a illustration of how to transition the calculation module to a robust modern approach?

By all means, do that.

 

[JimJCW]

Install python, man! . . .

That sounds like a good idea. Thanks for calculating the D_now values! A comparison is shown below:

z	Gnedin	Jorrie, 2022	ibix
0.1	1.41004	1.4084	1.41
0.09	1.27214	1.270609	1.27
0.08	1.13354	1.132106	1.13
0.07	0.99426	0.992931	0.994
0.06	0.85427	0.85305	0.854
0.05	0.7136	0.712477	0.713
0.04	0.57225	0.571225	0.572
0.03	0.43023	0.429278	0.43
0.02	0.28748	0.286651	0.28
0.01	0.14408	0.14336	0.144
0	0	0.000621	0
-0.01	0.14474	0.14526	0.145
-0.02	0.29015	0.290567	0.29
-0.03	0.43623	0.436536	0.436
-0.04	0.58295	0.583161	0.583
-0.05	0.73032	0.730433	0.73
-0.06	0.87834	0.878358	0.878
-0.07	1.02699	1.026926	1.03
-0.08	1.17631	1.176124	1.18
-0.09	1.32623	1.325968	1.33
-0.1	1.47679	1.476416	1.48

If we plot the three D_now vs z curves, they are just about indistinguishable. If we add the (D_Gnedin – D_ibix) result to the plot shown in Post #36, we get the following:

In general, (D_Gnedin - D_ibix) fluctuates around the y = 0 line. The greater deviations from the line near the two ends are probably due to round-off effect of the calculated D_ibix. I am still puzzled by the slightly different behaviors of the two data sets.

 

[pbuk]

I am still puzzled by the slightly different behaviors of the two data sets.

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.

 

[Ibix]

in short it needs rewriting.

It looks to me that Calculate() can be split up into three steps

1. Compute a list of the $z$ values requested, making sure $z=0$ is there
2. Do the integrals $\int_0^S\frac{1}{H(s)}ds$ and $\int_0^S\frac{1}{sH(s)}ds$ for all $S$ corresponding to those $z$ and also for $S=\infty$.
3. Compute all output statistics for each $z$ (some of which depend on the integral up to $S=1$ or $S=\infty$, hence doing this after the integration run).

I would expect that this is just reordering existing code, so ought to make little to no difference to the output. Then we could work on doing the integration in a better way.

Any thoughts? I haven't looked at your github page yet @pbuk, but will do so when I'm on something with a slightly larger screen.

 

[Jorrie]

Sounds good. The present Calculate() also has a separate T_age module, simply because there is a legacy requirement to immediate see the effect on present age when any relevant input parameter is changed. But this could be integrated with the $\int_0^S\frac{1}{H(s)}ds$ module.

 

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

 

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