SourceSpectrum integration not consistent if wavelength units vary

[astronomyk]

Dear pysynphot team,

firstly thanks for the package! Super useful.

Now the possible bug, either with pysynphot, or my understanding of things: The integrated flux over a specific wavelength range should remain constant in energy regardless of the wavelength units being used, right? (or is my education lacking in some aspect)

When I convert between wavelength units (ang to um), the integrated flux varies by the scaling factor between the two wavelength units (i.e. 1E4)

psp.setref(waveset=(0.1E4, 20E4, 1E4)) ## in Angstrom
bb = psp.BlackBody(3000)

bb.convert("Angstrom")
flam = bb.integrate("Flam")
print("Flam_integrated [erg s-1 cm-2]", flam)

resulting in

Flam_integrated [erg s-1 cm-2] 2.065223378057641

Where as

bb.convert("um")
flam = bb.integrate("Flam")
print("Flam_integrated [erg s-1 cm-2]", flam)

gives me

Flam_integrated [erg s-1 cm-2] 0.00020652233780576413

I'm assuming this comes from the fact that the integrator uses the units seen by the user, and not the internal units (i.e. Angstrom and Photlam).

I think the offending line of code is 571 from spectrum.py : SourceSpectrum.integrate()

wave, flux = self.getArrays()

because according to the doc string from getArrays

"""Return wavelength and flux arrays in user units.

instead of using internal units.

However, this functionality might be intentional and it is my understanding of Flux that has the issue. If so, I would be very grateful if you could point out what I missed / forgot during my Astronomy 101 classes.

Thanks in advance! Kieran

[pllim]

Hello. Thanks for reporting this. I agree with you that this is a bug, as changing the wavelength unit should not change flux integration result. The workaround is to only integrate in Angstrom space. That one is tested sufficiently and should be correct. If you have to convert wavelength to micron for other operations, please convert it back to Angstrom prior to integration.

In your example though, I think you are missing E-12 in your answer. I get 2.334850283236709e-12 from PySynphot (using the Angstrom case). Also the number, besides the exponent, does not quite agree between us. So either you did not provide a fully standalone example or you might be using a very old version where something else had changed since.

Alternately, there is also https://synphot.readthedocs.io , which is an Astropy affiliated package (www.astropy.org). Here is how to use it to accomplish similar calculations as you wrote above:

>>> import numpy as np
>>> from astropy import units as u
>>> from synphot import SourceSpectrum
>>> from synphot.models import BlackBodyNorm1D
>>> bb = SourceSpectrum(BlackBodyNorm1D, temperature=3000)
>>> w = np.logspace(3, np.log10(2e5), num=10000) * u.AA   # Based on your setref
>>> bb.integrate(wavelengths=w, flux_unit='flam')
<Quantity 2.33323452e-12 erg / (cm2 s)>
>>> w_um = w.to(u.micron)  # If you really need micron
>>> bb.integrate(wavelengths=w_um, flux_unit='flam')
<Quantity 2.33323452e-12 erg / (cm2 s)>

Just be advised of spacetelescope/synphot_refactor#166 that is still being worked on if you use effstim in synphot.

[astronomyk]

And again, I'm super impressed at how quick you are at replying to issues! Good stuff!

Indeed you're right about the exponent and the number. Not quite sure what I was doing there - in the meantime I've deleted that cell from the python notebook.

Thanks for the heads-up regarding Synphot. Would you recommend switching to that package, given that it's astropy affiliated, or sticking to Pysynphot? Which one will enjoy more attention from the developers in the future?

Thanks again! Kieran

[pllim]

If you are just starting out (i.e., don't have a lot of "legacy" code) and use Astropy a lot (i.e., units, models), I would recommend synphot over PySynphot. Thank you very much for your interest!