astropy: BlackBody bolometric flux is wrong if scale has units of dimensionless_unscaled

The astropy.modeling.models.BlackBody class has the wrong bolometric flux if scale argument is passed as a Quantity with dimensionless_unscaled units, but the correct bolometric flux if scale is simply a float.

Description

Expected behavior

Expected output from sample code:

4.823870774433646e-16 erg / (cm2 s)
4.823870774433646e-16 erg / (cm2 s)

Actual behavior

Actual output from sample code:

4.5930032795393893e+33 erg / (cm2 s)
4.823870774433646e-16 erg / (cm2 s)

Steps to Reproduce

Sample code:

from astropy.modeling.models import BlackBody
from astropy import units as u
import numpy as np

T = 3000 * u.K
r = 1e14 * u.cm
DL = 100 * u.Mpc
scale = np.pi * (r / DL)**2

print(BlackBody(temperature=T, scale=scale).bolometric_flux)
print(BlackBody(temperature=T, scale=scale.to_value(u.dimensionless_unscaled)).bolometric_flux)

System Details

>>> import numpy; print("Numpy", numpy.__version__)
Numpy 1.20.2
>>> import astropy; print("astropy", astropy.__version__)
astropy 4.3.dev758+g1ed1d945a
>>> import scipy; print("Scipy", scipy.__version__)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ModuleNotFoundError: No module named 'scipy'
>>> import matplotlib; print("Matplotlib", matplotlib.__version__)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ModuleNotFoundError: No module named 'matplotlib'

About this issue

  • Original URL
  • State: closed
  • Created 3 years ago
  • Comments: 34 (34 by maintainers)

Commits related to this issue

Most upvoted comments

Would we ruffle too many feathers

Can’t be worse than the episode when we deprecated clobber in io.fits… 😅

synphot keeps its own BlackBody1D class (perhaps renamed to BlackBodyFlux1D to mirror ConstFlux1D)

synphot never used the new blackbody stuff here, so I think it can be safely left out of the changes here. If you feel strongly about its model names, feel free to open issue at https://github.com/spacetelescope/synphot_refactor/issues but I don’t think it will affect anything at astropy or vice versa. 😅

@eteq @pllim - it might be possible to achieve this same use-case (not having to worry about thinking about solid angle if you don’t intend to make calls to bolometric_flux) in a single class by allowing solid_angle = None for the flux case and absorbing the steradians into the scale factor. That case would then need to raise an informative error for calls to bolometric_flux to avoid the ambiguity issue. The tradeoff I see is more complex argument validation logic and extended documentation in a single class rather than multiple classes for different use-cases.

If no one thinks of any major drawbacks/concerns, I will take a stab at that implementation and come up with examples for each of the use-cases discussed so far and we can then reconsider if splitting into separate classes is warranted.

Thanks for all the good ideas!

@kecnry, I’m concerned that overloading the scale to handle either a unitless value or a value with units of steradians is a footgun, because depending on the units you pass, it may or may not add a factor of pi. This is a footgun because people often think of steradians as being dimensionless.

No, I think @astrofrog introduced scale for fitting. The functional, uh, functions that we have removed did not have scaling.