hombit / mesa2py

Yet another Python binding to MESA
0 stars 1 forks source link

Zero abundances is chemical composition. #5

Closed AndreyTavleev closed 2 years ago

AndreyTavleev commented 2 years ago

Wrong processing of zero abundances, when setting up the chemical composition:

from opacity import Opac
op = Opac({'h1': 1.0, 'he4': 0.0})

returns

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "opacity.pyx", line 247, in opacity.Opac.__new__
  File "opacity.pyx", line 243, in genexpr
  File "opacity.pyx", line 243, in genexpr
  File "opacity.pyx", line 241, in genexpr
  File "opacity.pyx", line 216, in opacity.signif
ValueError: math domain error

The error is in opacity.pyx in signif function (calculation of log10(0)). The solution is to change two lines:

def signif(x, p):
    """Inspired by https://stackoverflow.com/a/59888924/5437597"""
    x_positive = np.where(x != 0, np.abs(x), 10**p)  # new line
#    mags = 10 ** (p - m.floor(m.log10(x)))
    mags = 10 ** (p - m.floor(m.log10(x_positive)))  # change line
    return round(x * mags) / mags

Also why do you use math instead of numpy? math is used only here in two places.

hombit commented 2 years ago

Do we need to check for infinites, NaNs and subnormal numbers (not sure the last one is possible in python)?

hombit commented 2 years ago

Why do we need np.abs?

AndreyTavleev commented 2 years ago

I realised, that we check whether composition values are non-negative before the signif() call. So we don't need np.abs. We can additionally check composition for infty and NaN before the normalisation calculation.