BlackBody model

[karllark]

BlackBody model that is a straightforward implementation of Planck's Law. This implementation is more flexible than the existing BlackBody1D model.

Added surface brightness equivalencies to the spectral_density equivalencies. Useful for BlackBody and in general.

Docs and tests are included in this PR.

The BlackBody is the 1st model located in physical_models.py. Other physics motivated models should go in this file (distinguished from functional_models.py that are those that are more mathematically motivated). A later PR can/could move some models from functional to physical as appropriate.

This new function replaces the existing blackbody_lambda and blackbody_nu functions and BlackBody1D model. See #9066 for discussion of why this is desirable. These existing functions have been deprecated as part of this PR. Documentation showing how to reproduce the results from all three with BlackBody has been included.

Closes #9066

[pllim]

Please do not merge until I have a chance to review this fully and test synphot against it. Thanks!

[bsipocz]

@pllim - please remove the label once you're happy with it downstream.

todo for the bot: stop the merge button to appear while a PR has the needs-downstream-testing label.

[pllim]

I'll start with nitpicking the change log. 😉

More to come...

[pllim]

Partial review on user documentation

During the hack day, we also talked about maybe examples of more advanced usage like redshifted and modified blackbodies would be nice.

[pllim]

Partial review on tests

I don't see any tests that apply non-default scales with different units.

Might also want to run test_blackbody.py with warning suppression off (see setup.cfg) and actively catch and test for AstropyDeprecationWarning (a grep for this in tests in other submodules should give you an idea).

EDIT: On warning suppression, also see this:

astropy/astropy/conftest.py

Line 22 in 09b7796

include_astropy_deprecations=False,

[pllim]

As discussed offline, the overflow test needs to be adapted for the refactored module as well.

Also, I wonder if this test needs to be ported over as well. It was originally adapted from #1480.

astropy/astropy/modeling/tests/test_blackbody.py

Lines 56 to 73 in 09b7796

	@pytest.mark.skipif('not HAS_SCIPY')
	def test_blackbody_scipy():
	"""Test Planck function.
	.. note:: Needs ``scipy`` to work.
	"""
	flux_unit = u.Watt / (u.m ** 2 * u.um)
	wave = np.logspace(0, 8, 100000) * u.AA
	temp = 100. * u.K
	with np.errstate(all='ignore'):
	bb_nu = blackbody_nu(wave, temp) * u.sr
	flux = bb_nu.to(flux_unit, u.spectral_density(wave)) / u.sr
	lum = wave.to(u.um)
	intflux = integrate.trapz(flux.value, x=lum.value)
	ans = const.sigma_sb * temp ** 4 / np.pi
	np.testing.assert_allclose(intflux, ans.value, rtol=0.01) # 1% accuracy

[pllim]

Note for the future: This refactoring gets rid of the ability to fit by bolometric flux (#4855) and also the lambda_max property. I am neutral on this. Just mainly a note for my future self.

[pllim]

I have done a full review and given (probably even too much?) feedback here. I have a companion synphot PR at spacetelescope/synphot_refactor#224 . I'll review again when all comments have been addressed.

[karllark]

@pllim : review (all of it) much appreciated. I'm doing a bit of work on the PR right now. You may be happy to learn that our in room conversations did convince me to put in lambda_max (and nu_max) as properties. Sometimes, I just need time. ;-)

[nden]

@astrofrog Since you wrote the original one - do you have time to review this?

[karllark]

I don't understand the issue with building the docs. They work fine when I create them on my laptop. Maybe a version issue? And I can't figure out the units issue from the error log.

[pllim]

Any situation where your example would generate something in 1 / (Hz sr)?

astropy.units.core.UnitConversionError: '1 / (Hz sr)' and 'erg / (cm2 Hz s)' (spectral flux density) are not convertible

Did you clean your doc build before rebuilding? Sometimes you could be skipping that example because an older (successful) cache already exists.

Anyway, looks like you need a rebase, as there is a conflict in change log.

[karllark]

Anyway, looks like you need a rebase, as there is a conflict in change log.

Just did this earlier today for the same reason. Then another modeling PR was merged making a conflict reappear. I'll wait a little bit and see if there are any other changes needed and then rebase. :-)

[nden]

@karllark Sorry for the multiple conflicts caused. That's it from me today. Please squash your commits when you rebase. This looks good and as far as I can tell it's ready to merge but I'm leaving the honor to @pllim .

[karllark]

Any situation where your example would generate something in 1 / (Hz sr)?

astropy.units.core.UnitConversionError: '1 / (Hz sr)' and 'erg / (cm2 Hz s)' (spectral flux density) are not convertible

Did you clean your doc build before rebuilding? Sometimes you could be skipping that example because an older (successful) cache already exists.

Cleaned, rebuild fine. Checked plotting code manually by running it in file. No idea where this error is from. I'm going to rebase, squash down to a reasonable set of commits, and then push. Maybe the error will go away (unlikely). If it doesn't, I don't know what to try next. Maybe someone else can clone this branch and investigate?

[bsipocz]

I haven't yet tested the failing stuff with this branch, but couldn't not notice that you add docs for the deprecated functionality? Is that necessary?

I mean we shouldn't really advertize and thus document now deprecated stuff when there is another suggested way to get to the same result.

[bsipocz]

OK, my local docs build is not working (due to other reasons), but the example plots works as expected on python3.7 and numpy dev and mpl dev.

[bsipocz]

Based on the tracebacks and generated artefacts, I'm fairly certain that the failure is triggered by the first example, specifically this rendering of it:
https://47429-2081289-gh.circle-artifacts.com/0/home/circleci/project/docs/_build/html/modeling/blackbody_deprecated-1.py

I cannot make it fail locally though with any of the version combinations I tried so far.

[pllim]

Ah, I get it! Somewhere else in the documentation build (remember that namespace is shared for all the documentation pages), this happened:

>>> from astropy.visualization import quantity_support
>>> quantity_support()  

As a result, it tried to convert bb(wavelengths)/bb.bolometric_flux (this has unit of 1 / (Hz sr)) to the unit of bb1d(wavelengths) (this has unit of erg / (cm2 Hz s)) and failed. So, it emits a warning with the traceback but plots the stuff anyway. In our CI, we convert all Sphinx warnings into failures. And there you have it.

That said, if bb(wavelengths)/bb.bolometric_flux and bb1d(wavelengths) are supposed to be equivalent, how do we explain that when they don't have the same units? @karllark , I know you tried to explain it to me but it slipped my grasp again. 😬

[bsipocz]

Oh, great. Once we find some of those infinite free time we might want to fix that docs pages are independent, etc.

However, my question about advertizing the deprecated behaviour still stands. Shouldn't it be removed fully?

[pllim]

To answer @bsipocz , "However, my question about advertizing the deprecated behaviour still stands. Shouldn't it be removed fully?" My answer is, maybe not this round, as it is @astrofrog 's request for examples to illustrate how the new BlackBody model can replace deprecated class and functions. Perhaps we can open a follow up issue to remove the deprecated examples in 4.1 or 4.2?

I agree with @nden that this is pretty much ready, except that the broken HTML job needs to be fixed first.

[pllim]

Interesting... I wouldn't have thought of using scale in that way. I would have done something more involved, like:

 bb_nu(wavelengths).to(u.erg / (u.cm ** 2 * u.AA * u.s * u.sr), u.spectral_density(wavelengths))

So, thanks for including this example!

[karllark]

Using scale this way make it "easier" for the user and (after my latest commit) also means that the scale units can be used to set the output units even during fitting. This can't be done using the return_units functionality of modeling as spectral_density equivalencies is needed and this isn't supported by return_units at this point.

[karllark]

Latest commit should fix the html docs bug (if I understood it correctly). Minor text cleanup as well.

This commit also means that the results are returned in scale units even during fitting (been trying to figure out how to make this happen for sometime and a solution finally dropped into my brain).

Hopefully this is the last commit needed (fingers-crossed).

[nden]

Everything looks in order and approved, so I'm merging it. Thanks @karllark !