Merge d0ea38e into 10f3843

CHANGES.rst

@@ -1,6 +1,11 @@
 0.2.0 (unreleased)
 ==================
 
+0.1.2 (unreleased)
+==================
+
+Bug fix for GaussianFlux1D ``total_flux`` unit handling. [#154]
+
 0.1.1 (2018-06-05)
 ==================
 

docs/index.rst

@@ -176,9 +176,10 @@ work, so they can be used directly. When in doubt, see if a model is in
     # 4000 Angstrom with a sigma of 20 Angstrom
     g_abs = BaseUnitlessSpectrum(GaussianAbsorption1D, amplitude=0.8,
                                  mean=4000, stddev=20)
-    # Create a Gaussian emission line with a total flux of 0.05 PHOTLAM
+    # Create a Gaussian emission line with the given total flux
     # centered at 3000 Angstrom with FWHM of 100 Angstrom
-    g_em = SourceSpectrum(GaussianFlux1D, total_flux=0.05*units.PHOTLAM,
+    g_em = SourceSpectrum(GaussianFlux1D,
+                          total_flux=3.5e-13*u.erg/(u.cm**2 * u.s),
                           mean=3000, fwhm=100)
     # Create a blackbody source spectrum with a temperature of 6000 K
     bb = SourceSpectrum(BlackBodyNorm1D, temperature=6000)
@@ -205,17 +206,15 @@ Models that built the spectrum::
     Model set size: 1
     Expression: ([0] + ([1] | [2])) * [3]
     Components:
-        [0]: <GaussianFlux1D(amplitude=0.0004697186393498257, mean=3000.0, ...>
+        [0]: <GaussianFlux1D(amplitude=...>
 
         [1]: <BlackBodyNorm1D(temperature=6000.0)>
 
         [2]: <Scale(factor=2.0)>
 
         [3]: <GaussianAbsorption1D(amplitude=0.8, mean=4000.0, stddev=20.0)>
     Parameters:
-          amplitude_0    mean_0    stddev_0   ... amplitude_3 mean_3 stddev_3
-        ---------------- ------ ------------- ... ----------- ------ --------
-        0.00046971863935 3000.0 42.4660900144 ...         0.8 4000.0     20.0
+        ...
 
 Redshift the source spectrum by :math:`z = 0.2`::
 

docs/synphot/from_pysyn_iraf.rst

@@ -95,7 +95,7 @@ Source Spectrum
 |alpha=expon)                          |                                      |
 +--------------------------------------+--------------------------------------+
 |SourceSpectrum(GaussianFlux1D, mean=mu|GaussianSource(flux, mu, fwhm,        |
-|, fwhm=fwhm, total_flux=flux*form)    |fluxunits=form)                       |
+|, fwhm=fwhm, total_flux=flux)         |fluxunits=form)                       |
 +--------------------------------------+--------------------------------------+
 |sp.taper()                            |sp.taper()                            |
 +--------------------------------------+--------------------------------------+
@@ -219,7 +219,7 @@ Source Spectrum
 |alpha=expon)                          |                                      |
 +--------------------------------------+--------------------------------------+
 |SourceSpectrum(GaussianFlux1D, mean=mu|em(mu, fwhm, flux, form)              |
-|, fwhm=fwhm, total_flux=flux*form)    |                                      |
+|, fwhm=fwhm, total_flux=flux)         |                                      |
 +--------------------------------------+--------------------------------------+
 |sp.normalize(val*form, band=bp)       |rn(sp, bp, val, form)                 |
 +--------------------------------------+--------------------------------------+

docs/synphot/spectrum.rst

@@ -194,11 +194,12 @@ Astropy). By default, integration is done in internal units::
 
     >>> from synphot import SourceSpectrum, units
     >>> from synphot.models import GaussianFlux1D
-    >>> sp = SourceSpectrum(GaussianFlux1D, mean=6000, fwhm=10, total_flux=1)
+    >>> sp = SourceSpectrum(GaussianFlux1D, mean=6000, fwhm=10,
+    ...                     total_flux=1*u.erg/(u.cm**2 * u.s))
     >>> sp.integrate()
-    <Quantity 0.99999922 ph / (cm2 s)>
+    <Quantity 3.02046763e+11 ph / (cm2 s)
     >>> sp.integrate(flux_unit=units.FLAM)
-    <Quantity 3.31074211e-12 erg / (cm2 s)>
+    <Quantity 0.99999972 erg / (cm2 s)>
 
 
 .. _synphot-empirical-source:
@@ -350,7 +351,8 @@ a continuum (e.g., :ref:`synphot-flat-spec`):
 
     from synphot import SourceSpectrum
     from synphot.models import GaussianFlux1D, ConstFlux1D
-    em = SourceSpectrum(GaussianFlux1D, mean=6000, fwhm=10, total_flux=1)
+    em = SourceSpectrum(GaussianFlux1D, mean=6000, fwhm=10,
+                        total_flux=3.3e-12*u.erg/(u.cm**2 * u.s))
     bg = SourceSpectrum(ConstFlux1D, amplitude=2)
     sp = bg - em
     sp.plot()
@@ -392,7 +394,7 @@ using input values (central wavelength, FWHM, and total flux) that are somewhat
 consistent with ASTROLIB PYSYNPHOT.
 
 The example below creates and plots a Gaussian source centered at 1.8 micron
-with FWHM of 200 nm and total flux of 18.3 ABMAG. As stated in
+with FWHM of 200 nm and the given total flux. As stated in
 :ref:`synphot_overview`, conversion to internal units happen behind the scenes:
 
 .. plot::
@@ -402,7 +404,7 @@ with FWHM of 200 nm and total flux of 18.3 ABMAG. As stated in
     from synphot import SourceSpectrum
     from synphot.models import GaussianFlux1D
     sp = SourceSpectrum(GaussianFlux1D, mean=1.8*u.micron, fwhm=200*u.nm,
-                        total_flux=18.3*u.ABmag)
+                        total_flux=1e-26*u.W/(u.m**2))
     sp.plot(title=sp.meta['expr'])
 
 

docs/synphot/tutorials.rst

@@ -31,7 +31,8 @@ FWHM of 100 Angstrom, and centered at 7000 Angstrom::
     >>> from synphot import units
     >>> from synphot.models import GaussianFlux1D
     >>> em = SourceSpectrum(
-    ...     GaussianFlux1D, total_flux=8e-14*units.FLAM, fwhm=100, mean=7000)
+    ...     GaussianFlux1D, total_flux=8e-14*u.erg/(u.cm**2 * u.s),
+    ...     fwhm=100, mean=7000)
 
 Add emission line to continuum spectrum::
 
@@ -82,10 +83,11 @@ continuum-normalized spectrum.
     sp_nse = BaseUnitlessSpectrum(Empirical1D, points=wave, lookup_table=nse)
     bb_noisy = bb * sp_nse
     # Apply emission and absorption lines to the noisy continuum.
+    tf_unit = u.erg / (u.cm**2 * u.s)
     g_em = SourceSpectrum(
-        GaussianFlux1D, total_flux=0.02, mean=15000, fwhm=500)
+        GaussianFlux1D, total_flux=2.65e-14*tf_unit, mean=15000, fwhm=500)
     g_ab = SourceSpectrum(
-        GaussianFlux1D, total_flux=0.015, mean=4500, fwhm=100)
+        GaussianFlux1D, total_flux=6.62e-14*tf_unit, mean=4500, fwhm=100)
     sp = bb_noisy + g_em - g_ab
     # Divide the noisy spectrum with lines with the original smooth continuum
     # to obtain the continuum-normalized spectrum.

synphot/models.py

@@ -446,7 +446,7 @@ class Gaussian1D(BaseGaussian1D):
 
 
 class GaussianAbsorption1D(BaseGaussian1D):
-    """Same as `astropy.modeling.functional_models.GaussianAbsorption1D`,
+    """Same as ``astropy.modeling.functional_models.GaussianAbsorption1D``,
     except with ``sampleset`` defined.
 
     """
@@ -472,12 +472,23 @@ class GaussianFlux1D(Gaussian1D):
 
     Parameters
     ----------
+    amplitude : float
+        Amplitude of the Gaussian in PHOTLAM.
+        Also see ``total_flux``.
+
+    mean : float
+        Mean of the Gaussian in Angstrom.
+
+    stddev : float
+        Standard deviation of the Gaussian in Angstrom.
+        Also see ``fwhm``.
+
     fwhm : float
         Full width at half maximum of the Gaussian in Angstrom.
         If given, this overrides ``stddev``.
 
     total_flux : float
-        Total flux under the Gaussian in PHOTLAM.
+        Total flux under the Gaussian in ``erg/s/cm^2``.
         If given, this overrides ``amplitude``.
 
     """
@@ -495,12 +506,21 @@ def __init__(self, *args, **kwargs):
         gaussian_amp_to_totflux = np.sqrt(2.0 * np.pi) * self.stddev
 
         if total_flux is None:
+            u_str = 'PHOTLAM'
             total_flux = self.amplitude * gaussian_amp_to_totflux
         else:
-            self.amplitude = total_flux / gaussian_amp_to_totflux
+            u_str = 'FLAM'
+            # total_flux is passed in unaltered, any conversion error would
+            # happen here.
+            tf_unit = u.erg / (u.cm * u.cm * u.s)
+            if isinstance(total_flux, u.Quantity):
+                total_flux = total_flux.to(tf_unit)
+            else:
+                total_flux *= tf_unit
+            self.amplitude = (total_flux / (gaussian_amp_to_totflux * u.AA)).to(units.PHOTLAM, u.spectral_density(self.mean.value * u.AA)).value  # noqa
 
-        self.meta['expr'] = 'em({0:g}, {1:g}, {2:g}, PHOTLAM)'.format(
-            self.mean.value, fwhm, total_flux)
+        self.meta['expr'] = 'em({0:g}, {1:g}, {2:g}, {3})'.format(
+            self.mean.value, fwhm, total_flux.value, u_str)
 
 
 class Lorentz1D(_models.Lorentz1D):

synphot/spectrum.py

@@ -88,7 +88,7 @@ class BaseSpectrum(object):
         'Gaussian1D': {'amplitude': 'flux', 'mean': 'wave', 'stddev': 'wave'},
         'GaussianAbsorption1D': {
             'amplitude': 'flux', 'mean': 'wave', 'stddev': 'wave'},
-        'GaussianFlux1D': {'total_flux': 'flux', 'amplitude': 'flux',
+        'GaussianFlux1D': {'total_flux': 'noconv', 'amplitude': 'flux',
                            'mean': 'wave', 'stddev': 'wave', 'fwhm': 'wave'},
         'LogParabola1D': {
             'amplitude': 'flux', 'x_0': 'wave',

synphot/tests/test_observation.py

@@ -141,8 +141,9 @@ def test_default_binset_from_bandpass(self):
         np.testing.assert_array_equal(obs2.binset, self.obs.bandpass.waveset)
 
     def test_default_binset_from_spectrum(self):
+        tf_unit = u.erg / (u.cm * u.cm * u.s)
         sp = SourceSpectrum(
-            GaussianFlux1D, mean=5000, total_flux=(1 * units.FLAM), fwhm=10)
+            GaussianFlux1D, mean=5000, total_flux=(1 * tf_unit), fwhm=10)
         bp = SpectralElement(Const1D, amplitude=1)
         obs2 = Observation(sp, bp, force='extrap')
         np.testing.assert_array_equal(obs2.binset, sp.waveset)

synphot/tests/test_spectrum.py

@@ -381,8 +381,10 @@ def test_eval(self):
 class TestGaussianSource(object):
     """Test source spectrum with GaussianFlux1D model."""
     def setup_class(self):
+        self._tf_unit = u.erg / (u.cm * u.cm * u.s)
+        tf = 4.96611456e-12 * self._tf_unit
         self.sp = SourceSpectrum(
-            GaussianFlux1D, total_flux=1, mean=4000, fwhm=100)
+            GaussianFlux1D, total_flux=tf, mean=4000, fwhm=100)
 
     def test_eval(self):
         y = self.sp([3900, 4000, 4060])
@@ -397,7 +399,7 @@ def test_totalflux(self):
 
         # FLAM
         x0 = (400 * u.nm).to(u.AA).value
-        totflux = units.convert_flux(x0, 1 * units.FLAM, units.PHOTLAM).value
+        totflux = 1 * self._tf_unit
         fwhm = (10 * u.nm).to(u.AA)
         sp2 = SourceSpectrum(
             GaussianFlux1D, total_flux=totflux, mean=x0, fwhm=fwhm)
@@ -415,6 +417,17 @@ def test_fwhm(self):
         np.testing.assert_allclose(bp.fwhm().value, 100, rtol=1e-3)  # 0.1%
 
 
+def test_gaussian_source_watts():
+    """https://github.com/spacetelescope/synphot_refactor/issues/153"""
+    mu = 1 * u.um
+    fwhm = (0.01 / 0.42466) * u.um
+    flux = 1 * u.W / u.m**2
+
+    sp = SourceSpectrum(GaussianFlux1D, mean=mu, fwhm=fwhm, total_flux=flux)
+    tf = sp.integrate(flux_unit=units.FLAM)
+    np.testing.assert_allclose(tf.to(flux.unit), flux, rtol=1e-4)
+
+
 class TestPowerLawSource(object):
     """Test source spectrum with PowerLawFlux1D model."""
     def setup_class(self):
@@ -597,8 +610,9 @@ def setup_class(self):
         self.bb = SourceSpectrum(BlackBodyNorm1D, temperature=5000)
 
         # Gaussian emission line: em(5500, 250, 1e-13, flam)
+        tf_unit = u.erg / (u.cm * u.cm * u.s)
         self.em = SourceSpectrum(GaussianFlux1D, mean=5500,
-                                 total_flux=(1e-13 * units.FLAM), fwhm=250)
+                                 total_flux=(1e-13 * tf_unit), fwhm=250)
 
         # ACS bandpass: band(acs,hrc,f555w)
         bandfile = get_pkg_data_filename(
@@ -723,8 +737,9 @@ def test_none(self):
         assert sp.waveset is None
 
     def test_sampleset(self):
+        tf_unit = u.erg / (u.cm * u.cm * u.s)
         sp = SourceSpectrum(
-            GaussianFlux1D, total_flux=(1 * units.FLAM), mean=5000, fwhm=10)
+            GaussianFlux1D, total_flux=(1 * tf_unit), mean=5000, fwhm=10)
         np.testing.assert_array_equal(sp.waveset.value, sp.model.sampleset())
 
     def test_box1d(self):
@@ -742,7 +757,7 @@ def test_composite_none(self):
         np.testing.assert_array_equal(bp.waveset, bp1.waveset)
 
     def test_composite(self):
-        totflux = 1 * units.FLAM
+        totflux = 1 * (u.erg / (u.cm * u.cm * u.s))
         g1 = SourceSpectrum(
             GaussianFlux1D, total_flux=totflux, mean=5000, fwhm=10)
         g2 = SourceSpectrum(
@@ -755,8 +770,9 @@ def test_composite(self):
             [999.49, 1000.49, 5019.95906231, 6699.59062307, 7509.7672007])
 
     def test_redshift(self):
+        tf_unit = u.erg / (u.cm * u.cm * u.s)
         sp = SourceSpectrum(
-            GaussianFlux1D, total_flux=(1 * units.FLAM), mean=5000, fwhm=10)
+            GaussianFlux1D, total_flux=(1 * tf_unit), mean=5000, fwhm=10)
         sp.z = 1.3
         m = RedshiftScaleFactor(z=1.3)
         w_step25_z0 = [4978.76695499, 4989.3834775, 5000, 5010.6165225]
@@ -779,7 +795,7 @@ class TestRedShift(object):
     """
     def setup_class(self):
         x0 = 5000
-        totflux = 1 * u.Jy
+        totflux = 1e-23 * (u.erg / (u.cm * u.cm * u.s))  # 1 Jy * Hz
         fwhm = 100
         self.sp_z0 = SourceSpectrum(
             GaussianFlux1D, total_flux=totflux, mean=x0, fwhm=fwhm)