Planck Function

[cmccully]

I believe the planck function bb_func in planck.py appears to be missing a factor of wavelength. At line 99, x = factor * x * x * x and then there is one more wavelength factor divided out at the end, resulting in a total of wavelength ** -4, not -5 which is correct scaling. Should that line have a *= instead? I haven't checked if the normalization is correct or not.

[pllim]

Pretty sure it is correct. The result has been tested against IRAF synphot and similar functionality in Astropy. If you think the results are wrong, please re-open and provide a complete example. Thanks.

[cmccully]

I do not have permissions to reopen an issue.

If the code is correct, then it needs more documentation because the bbfunc and llam_SI show different behavior (the peaks are different because they are different by a factor of the wavelength, lambda). Example is below (plot is attached):

from matplotlib import pyplot
from pysynphot import planck
import numpy as np
lam = np.arange(2000, 12001)
b = planck.bbfunc(lam, 10000)
b2 = planck.llam_SI(lam / 1e10, 10000)
pyplot.plot(l, b / np.mean(b))
pyplot.plot(l, b2 / np.mean(b2))
pyplot.show()

[pllim]

It's weird that you can't re-open since you created it in the first place, but I re-opened it for you anyway.

The functions in planck are not meant to be used directly, but rather used indirectly by pysynphot via its spectrum objects (e.g., blackbody source spectrum or background thermal calculations). Those two do not produce the same output units, hence the difference you see.

What exactly are you looking for here?

[cmccully]

I was trying to validate a simple black body fitting code that I wrote, but the peak of my blackbody code did not match the one produced by pysynphot. So based on the readthedocs, I tried the following:
In [32]: from pysynphot import spectrum

In [33]: bb = spectrum.BlackBody(5000)

In [34]: bb.wave[np.argmax(bb.flux)]
Out[34]: 7339.7400297101867

Thus the peak of the blackbody produced by pysynphot is 7339.74 angstroms. However, my code produces a peak of 579.6 nm. I checked this with the following hyperphysics link:
http://hyperphysics.phy-astr.gsu.edu/hbase/wien.html

This matches well. I ended up tracking it down to be how the code calls the bbfunc and I noticed that it is missing a factor of 1/lambda. If I plot lambda F_lambda, my curve matches the one produced by pysynphot exactly.

[pllim]

Ah, I see. I agree that your answer of 579.6 nm is correct. In fact, it agrees with Astropy:

>>> import numpy as np
>>> from astropy import units as u
>>> from astropy.analytic_functions import blackbody_lambda
>>> w = np.arange(1000, 10000) * u.AA
>>> flux = blackbody_lambda(w, 5000 * u.K)
>>> w[np.argmax(flux)]
<Quantity 5796.0 Angstrom>

With pysynphot, you need to actually convert the blackbody spectrum to flam first. In its default photlam unit, you will get the wrong answer (as you have found):

>>> import numpy as np
>>> import pysynphot as S
>>> bb = S.BlackBody(5000)
>>> bb.convert('flam')
>>> bb.wave[np.argmax(bb.flux)]
5795.1366044833067

Not quite exactly 5796 but much closer.

[pllim]

This is more of a note for myself... In the refactored synphot for Astropy, it can be done like this:

>>> from synphot.models import BlackBodyNorm1D
>>> bb = BlackBodyNorm1D(temperature=5000)  # Same as bb(5000) in IRAF
>>> bb.lambda_max
5795.544199999998

[pllim]

@cmccully , if this sufficiently addresses your concerns, feel free to close the issue.

[pllim]

@cmccully , I haven't heard from you lately, so closing this issue.

[cmccully]

@pllim Sorry for the delay. I have been traveling. I think this basically answers my question. I would say that it would be useful to have a little more documentation in the planck.py to explain what that function returns. It would also be good to do the conversion to flam in the readthedocs as that shows the lambda^-5 function and then plots a different function. For that matter, why is the default units the blackbody object returns not flam? That seems counterintuitive to me. Anyway, thanks for the explanation!

[pllim]

@cmccully , is there a reason why you prefer the planck.py here to the one in astropy.analytic_functions? I would recommend the latter as the one here wasn't meant for public consumption.

[cmccully]

@pllim I don't really have a preference. I was just trying to do synthetic photometry and was trying to diagnose some differences that I was seeing. I didn't realize that I would need to farm out to another package to make the blackbody function as there is an example in the readthedocs.

[pllim]

Ah, okay. A lot of this stuff was inherited from the synphot package in IRAF, hence why unit is not automatically in FLAM, etc. We're already working on a replacement of Pysynphot to use Astropy models (and naturally astropy.analytic_functions). I will re-open this issue as a reminder for updating documentation (but it won't be that soon). Of course, you are also welcome to submit a PR as well. Thank you for bringing this to our attention.