search

villano-lab / nrCascadeSim

calculating the NR spectrum resulting from neutron-capture cascades.
MIT License
0 stars 1 forks source link

[JOSS Review]: Possible loss of precision in cross-section computation #72

Open altavir opened 2 years ago

altavir commented 2 years ago

There are multiple places in the code where the loss of precision is possible like this one:

https://github.com/villano-lab/nrCascadeSim/blob/4f349d1de796e099db4af9f381b8ece16c894f33/src/cascadeProd.c#L237

It is not clear from the method documentation what is the mass unit. If it is is in GeVs, it seems to be within error boundary for double precision, but could still cause a loss of numerical precision if done in a cycle.

villaa commented 2 years ago

@altavir what does it mean "if done in a cycle?" --thanks

villaa commented 2 years ago

Probably the best thing to do is start with identifying the places this loss of precision can happen:

....analogous functions for other isotopes...

....and maybe the Weisskopf functions....

altavir commented 2 years ago

Consider a simple operation: a = a + b/c. If b and c are very small or very large, you will add a small numerical uncertainty to a. If it is done once - uncertainty is small. Say, 0.001. Then you do that in a cycle:

for (int i = 1; i <= 1000; ++i) {
       a = a + b/c;
}

The absolute uncertainties are summed and you get 50% uncerainty in the result.

For double-precision floating point, you usually need to think about number of order like 10e12 or 10e-12. For numbers like 1e9 it probably should be fine, but you need to check it. In the case above you have (Egam*1e6)^2. If you have so for keV photon, you will have like 10^18 which is large enough to be checked.

You need to remember that for floating point, the precision is relative, not absolute.

I usually either compare results with different precision, or write corner-case tests.

villaa commented 2 years ago

Ok, the first obvious thing is that many of these computations are just copied multiple times. That means we should make functions for them instead of copying lines of code. Useful functions:

Rmax(M, E, Ntilde...) g(escale,xvar) tmax R NS2 zfac Ecm ...various return values un-named....

Also consider making absolute constants like "pi" #define's