Skip to content

Commit

Permalink
Merge pull request #357 from SNEWS2/sybenzvi/suppress-div-by-zero
Browse files Browse the repository at this point in the history
Suppress div-by-zero warning in spectrum calc.
  • Loading branch information
JostMigenda authored Aug 14, 2024
2 parents c46df13 + 9148335 commit 400d951
Showing 1 changed file with 9 additions and 3 deletions.
12 changes: 9 additions & 3 deletions python/snewpy/models/base.py
Original file line number Diff line number Diff line change
Expand Up @@ -297,14 +297,20 @@ def get_initial_spectra(self, t, E, flavors=Flavor):
L = np.expand_dims(L, axis=1)
Ea = np.expand_dims(Ea,axis=1)
a = np.expand_dims(a, axis=1)

# For numerical stability, evaluate log PDF and then exponentiate.
result = \
np.exp(np.log(L) - (2+a)*np.log(Ea) + (1+a)*np.log(1+a)
- loggamma(1+a) + a*np.log(E) - (1+a)*(E/Ea)) / (u.erg * u.s)
# Suppress div-by-zero and other warnings that do not matter here.
with np.errstate(divide='ignore', invalid='ignore'):
result = \
np.exp(np.log(L) - (2+a)*np.log(Ea) + (1+a)*np.log(1+a)
- loggamma(1+a) + a*np.log(E) - (1+a)*(E/Ea)) / (u.erg * u.s)

#remove bad values
result[np.isnan(result)] = 0
result[:, E[0]==0] = 0

#remove unnecessary dimensions, if E or t was scalar:
result = np.squeeze(result)
initialspectra[flavor] = result

return initialspectra

0 comments on commit 400d951

Please sign in to comment.