BUG: parse_spec for "em" might be wrong if you use synphot>=0.1.2

[pllim]

synphot changed unit handling behavior for total flux input for GuassianFlux1D in spacetelescope/synphot_refactor#154 to fix another bug (spacetelescope/synphot_refactor#153). It was an oversight that parse_spec here was not updated to match.

For example: 'em(5000, 25, 1, flam)' should produce Gaussian with a total flux of 1 erg / s / cm^2 but because of how total_flux is passed in here (as a plain float), the assumption that synphot would interpret it as ph / s / cm^2 is now wrong.

I don't remember why I did not pass it in as a Quantity. I must gave a good reason a few years ago? :woman_shrugging: :grimacing:

[pllim]

I'll roll out a bug fix release this week. @vglaidler and @ibusko , I wonder if this bug caused some discrepancies you are seeing. I am not sure what versions you are using for your acceptance test.

[ibusko]

We are using

stsynphot 0.2.dev224 dev_0 <develop>

synphot 0.2.dev413 dev_0 <develop>

[pllim]

Then you might see discrepancy when you get your Gaussian spectrum using parse_spec('em(..., flam)'). If you use the regular Python API of SourceSpectrum(GaussianFlux1D, ...), then you are not affected. parse_spec('em(..., photlam)') should also be okay in theory.

[pllim]

To illustrate, let's pick a test case from #103 below:

from astropy import units as u
from astropy.tests.helper import assert_quantity_allclose  # For later

from synphot import SourceSpectrum, units
from synphot.models import GaussianFlux1D

from stsynphot import spparser

sp1 = spparser.parse_spec('em(5000, 25, 1, flam)')
w = sp1.waveset  # For later

# Pythonic equivalent
f = 1 * (units.FLAM * u.AA)
sp2 = SourceSpectrum(GaussianFlux1D, mean=5000 * u.AA, fwhm=25 * u.AA, total_flux=f)

With synphot 0.2.1 and stsynphot 0.2.1 (without the patch in #103):

>>> sp1.meta['expr']
'em(5000, 25, 2.51706e+11, FLAM)'
>>> assert_quantity_allclose(sp1(w), sp2(w))
...
AssertionError: 
Not equal to tolerance rtol=1e-07, atol=0

Mismatched elements: 100 / 100 (100%)
Max absolute difference: 2.38075292e+21
Max relative difference: 2.51705828e+11
 x: array([8.872240e+15, 1.455489e+16, 2.363970e+16, 3.801297e+16,
       6.051718e+16, 9.538557e+16, 1.488483e+17, 2.299651e+17,
       3.517523e+17, 5.326835e+17, 7.986536e+17, 1.185509e+18,...
 y: array([3.524845e+04, 5.782502e+04, 9.391795e+04, 1.510214e+05,
       2.404282e+05, 3.789566e+05, 5.913580e+05, 9.136262e+05,
       1.397474e+06, 2.116294e+06, 3.172964e+06, 4.709898e+06,...
>>> assert_quantity_allclose(sp1.integrate(wavelengths=w), sp2.integrate(wavelengths=w))
...
AssertionError: 
Not equal to tolerance rtol=1e-07, atol=0

Mismatched elements: 1 / 1 (100%)
Max absolute difference: 6.33557745e+22
Max relative difference: 2.51705828e+11
 x: array(6.335577e+22)
 y: array(2.517056e+11)

With the patch in #103 for stsynphot but same synphot version as above:

>>> sp1.meta['expr']
'em(5000, 25, 1, FLAM)'
>>> assert_quantity_allclose(sp1(w), sp2(w))  # No error
>>> assert_quantity_allclose(sp1.integrate(wavelengths=w), sp2.integrate(wavelengths=w))  # No error

[ibusko]

This is good to know. I will update my test environment with the new version. Thanks for fixing this!

[pllim]

Sorry for any inconvenience caused!