Merge pull request #151 from pllim/fix-integrate-unit

Fix units

docs/index.rst

@@ -193,13 +193,13 @@ Sample the spectrum at 0.3 micron::
 
     >>> from astropy import units as u
     >>> sp(0.3 * u.micron)
-    <Quantity 0.001268063356632316 PHOTLAM>
+    <Quantity 0.0012685 PHOTLAM>
 
 Models that built the spectrum::
 
     >>> print(sp)
     SourceSpectrum at z=0.0
-    Model: CompoundModel5
+    Model: CompoundModel2
     Inputs: ('x',)
     Outputs: ('y',)
     Model set size: 1
@@ -233,7 +233,7 @@ integrate it::
 
     >>> sp_rn = sp.normalize(1 * u.Jy, band=bp)
     >>> sp_rn.integrate()
-    <Quantity 12900.947986966807 PHOTLAM>
+    <Quantity 12901.13466265 ph / (cm2 s)>
 
 Create an observation by passing the redshifted and normalized source spectrum
 through the box bandpass::
@@ -246,9 +246,9 @@ Calculate the count rate of the observation above for an 2-meter telescope:
     >>> import numpy as np
     >>> area = np.pi * (1 * u.m) ** 2
     >>> area
-    <Quantity 3.141592653589793 m2>
+    <Quantity 3.14159265 m2>
     >>> obs.countrate(area=area)
-    <Quantity 24219652.902789623 ct / s>
+    <Quantity 24219651.1854991 ct / s>
 
 
 .. _synphot_using:

docs/synphot/formulae.rst

@@ -108,7 +108,7 @@ Bandpass Equivalent Width
 -------------------------
 
 For a bandpass, :meth:`~synphot.spectrum.SpectralElement.equivwidth` implements
-the equivalent width. It gives the same value (but not unit) as
+the equivalent width. It gives the same value as
 :meth:`~synphot.spectrum.BaseSpectrum.integrate` and is equivalent to
 IRAF SYNPHOT ``bandpar`` result for ``equvw``.
 
@@ -121,7 +121,7 @@ Example::
     >>> bp.equivwidth()
     <Quantity 272.0108162945954 Angstrom>
     >>> bp.integrate()
-    <Quantity 272.0108162945954>
+    <Quantity 272.0108162945954 Angstrom>
 
 
 .. _synphot-formula-rectw:
@@ -236,7 +236,7 @@ It is equivalent to IRAF SYNPHOT ``bandpar`` result for ``uresp``.
 
 .. math::
 
-    \text{uresp} = \frac{hc}{a} \int P_{\lambda}\; \lambda\; d\lambda
+    \text{uresp} = \frac{hc}{a \int P_{\lambda}\; \lambda\; d\lambda}
 
 Example::
 

docs/synphot/spectrum.rst

@@ -196,9 +196,9 @@ Astropy). By default, integration is done in internal units::
     >>> from synphot.models import GaussianFlux1D
     >>> sp = SourceSpectrum(GaussianFlux1D, mean=6000, fwhm=10, total_flux=1)
     >>> sp.integrate()
-    <Quantity 0.9999992180687505 PHOTLAM>
+    <Quantity 0.99999922 ph / (cm2 s)>
     >>> sp.integrate(flux_unit=units.FLAM)
-    <Quantity 3.3107418773847306e-12 FLAM>
+    <Quantity 3.31074211e-12 erg / (cm2 s)>
 
 
 .. _synphot-empirical-source:

synphot/spectrum.py

@@ -445,9 +445,16 @@ def integrate(self, wavelengths=None, **kwargs):
             Invalid default integrator.
 
         synphot.exceptions.SynphotError
-            `waveset` is needed but undefined.
+            `waveset` is needed but undefined or cannot integrate
+            natively in the given ``flux_unit``.
 
         """
+        # Cannot integrate per Hz units natively across wavelength
+        # without converting them to per Angstrom unit first, so
+        # less misleading to just disallow that option for now.
+        if 'flux_unit' in kwargs:
+            self._validate_flux_unit(kwargs['flux_unit'], wav_only=True)
+
         x = self._validate_wavelengths(wavelengths)
 
         # TODO: When astropy.modeling.models supports this, need to
@@ -459,17 +466,33 @@ def integrate(self, wavelengths=None, **kwargs):
         except (AttributeError, NotImplementedError):
             if conf.default_integrator == 'trapezoid':
                 y = self(x, **kwargs)
-                result = abs(np.trapz(y.value, x=x.value)) * y.unit
+                result = abs(np.trapz(y.value, x=x.value))
+                result_unit = y.unit
             else:  # pragma: no cover
                 raise NotImplementedError(
                     'Analytic integral not available and default integrator '
                     '{0} is not supported'.format(conf.default_integrator))
         else:  # pragma: no cover
             start = x[0].value
             stop = x[-1].value
-            result = (m(stop) - m(start)) * self._internal_flux_unit
+            result = (m(stop) - m(start))
+            result_unit = self._internal_flux_unit
+
+        # Ensure final unit takes account of integration across wavelength
+        if result_unit != units.THROUGHPUT:
+            if result_unit == units.PHOTLAM:
+                result_unit = u.photon / (u.cm**2 * u.s)
+            elif result_unit == units.FLAM:
+                result_unit = u.erg / (u.cm**2 * u.s)
+            else:  # pragma: no cover
+                raise NotImplementedError(
+                    'Integration of {0} is not supported'.format(result_unit))
+        else:
+            # Ideally flux can use this too but unfortunately this
+            # operation results in confusing output unit for flux.
+            result_unit *= self._internal_wave_unit
 
-        return result
+        return result * result_unit
 
     def avgwave(self, wavelengths=None):
         """Calculate the :ref:`average wavelength <synphot-formula-avgwv>`.
@@ -761,17 +784,22 @@ class BaseSourceSpectrum(BaseSpectrum):
     _internal_flux_unit = units.PHOTLAM
 
     @staticmethod
-    def _validate_flux_unit(new_unit):
+    def _validate_flux_unit(new_unit, wav_only=False):
         """Make sure flux unit is valid."""
         new_unit = units.validate_unit(new_unit)
-        acceptable_types = [
-            'spectral flux density', 'spectral flux density wav',
-            'photon flux density', 'photon flux density wav']
+        acceptable_types = ['spectral flux density wav',
+                            'photon flux density wav']
+        acceptable_names = ['PHOTLAM', 'FLAM']
+
+        if not wav_only:  # Include per Hz units
+            acceptable_types += ['spectral flux density',
+                                 'photon flux density']
+            acceptable_names += ['PHOTNU', 'FNU', 'Jy']
 
         if new_unit.physical_type not in acceptable_types:
             raise exceptions.SynphotError(
-                'Source spectrum cannot operate in {0}. Acceptable units are '
-                'PHOTLAM, PHOTNU, FLAM, FNU, or Jy.'.format(new_unit))
+                'Source spectrum cannot operate in {0}. Acceptable units: '
+                '{1}'.format(new_unit, ','.join(acceptable_names)))
 
         return new_unit
 
@@ -921,10 +949,10 @@ def normalize(self, renorm_val, band=None, wavelengths=None, force=False,
                 renorm_unit_name == units.OBMAG.to_string()):
             # Special handling for non-density units
             flux_tmp = sp(w, flux_unit=u.count, area=area)
-            totalflux = flux_tmp.sum()
+            totalflux = flux_tmp.sum().value
             stdflux = 1.0
         else:
-            totalflux = sp.integrate(wavelengths=wavelengths)
+            totalflux = sp.integrate(wavelengths=wavelengths).value
 
             # VEGAMAG
             if renorm_unit_name == units.VEGAMAG.to_string():
@@ -952,17 +980,17 @@ def normalize(self, renorm_val, band=None, wavelengths=None, force=False,
                 up = stdspec * band
                 stdflux = up.integrate(wavelengths=wavelengths).value
 
-        utils.validate_totalflux(totalflux.value)
+        utils.validate_totalflux(totalflux)
 
         # Renormalize in magnitudes
         if (renorm_val.unit.decompose() == u.mag or
                 isinstance(renorm_val.unit, u.LogUnit)):
             const = renorm_val.value + (2.5 *
-                                        np.log10(totalflux.value / stdflux))
+                                        np.log10(totalflux / stdflux))
             newsp = self.__mul__(10**(-0.4 * const))
         # Renormalize in linear flux units
         else:
-            const = renorm_val.value * (stdflux / totalflux.value)
+            const = renorm_val.value * (stdflux / totalflux)
             newsp = self.__mul__(const)
 
         newsp.warnings = warndict
@@ -1662,8 +1690,7 @@ def equivwidth(self, wavelengths=None):
             Bandpass equivalent width.
 
         """
-        return (self.integrate(wavelengths=wavelengths).value *
-                self._internal_wave_unit)
+        return self.integrate(wavelengths=wavelengths)
 
     def rectwidth(self, wavelengths=None):
         """Calculate :ref:`bandpass rectangular width <synphot-formula-rectw>`.
@@ -1724,7 +1751,7 @@ def emflx(self, area, wavelengths=None):
         Returns
         -------
         em_flux : `~astropy.units.quantity.Quantity`
-            Equivalent monochromatic flux in FLAM.
+            Equivalent monochromatic flux.
 
         """
         t_lambda = self.tlambda(wavelengths=wavelengths)
@@ -1733,9 +1760,8 @@ def emflx(self, area, wavelengths=None):
             em_flux = 0.0 * units.FLAM
         else:
             uresp = self.unit_response(area, wavelengths=wavelengths)
-            rectw = self.rectwidth(wavelengths=wavelengths).value
-            fac = self.tpeak(wavelengths=wavelengths) / t_lambda
-            em_flux = uresp * rectw * fac
+            equvw = self.equivwidth(wavelengths=wavelengths).value
+            em_flux = uresp * equvw / t_lambda
 
         return em_flux
 

synphot/tests/test_spectrum.py

@@ -173,13 +173,20 @@ def test_conversion(self):
         np.testing.assert_allclose(y.value, 2.282950185743497e-26, rtol=1e-6)
 
     def test_integrate(self):
+        expected_unit = u.erg / (u.cm**2 * u.s)
         # Whole range
         f = self.sp.integrate(flux_unit=units.FLAM)
         np.testing.assert_allclose(f.value, 8.460125829057308e-12, rtol=1e-5)
+        assert f.unit == expected_unit
 
         # Given range
         f = self.sp.integrate(wavelengths=_wave, flux_unit=units.FLAM)
         np.testing.assert_allclose(f.value, 4.810058069909525e-14, rtol=1e-5)
+        assert f.unit == expected_unit
+
+        # Unsupported unit
+        with pytest.raises(exceptions.SynphotError):
+            self.sp.integrate(flux_unit=u.Jy)
 
     def test_taper(self):
         # Original spectrum already tapered -- nothing done
@@ -222,10 +229,12 @@ def test_integrate(self):
         # Whole range (same as EQUVW)
         f = self.bp.integrate()
         np.testing.assert_allclose(f.value, 272.01081629459344)
+        assert f.unit == u.AA
 
         # Given range
         f = self.bp.integrate(wavelengths=_wave)
         np.testing.assert_allclose(f.value, 1.2062975715374322, rtol=2.5e-6)
+        assert f.unit == u.AA
 
     def test_avgwave(self):
         """Compare AVGWAVE with old SYNPHOT result."""
@@ -300,6 +309,7 @@ def test_equivwidth(self):
         """Compare EQUVW with ASTROLIB PYSYNPHOT result."""
         w = self.bp.equivwidth()
         np.testing.assert_allclose(w.value, 272.01081629459344, rtol=1e-6)
+        assert w.unit == u.AA
 
     def test_rectw(self):
         """Compare RECTW with old SYNPHOT result."""
@@ -381,8 +391,9 @@ def test_eval(self):
 
     def test_totalflux(self):
         # PHOTLAM
-        val = self.sp.integrate().value
-        np.testing.assert_allclose(val, 1, rtol=1e-5)
+        f = self.sp.integrate()
+        np.testing.assert_allclose(f.value, 1, rtol=1e-5)
+        assert f.unit == u.photon / (u.cm**2 * u.s)
 
         # FLAM
         x0 = (400 * u.nm).to(u.AA).value
@@ -660,16 +671,6 @@ def test_renorm_vegamag(self, sp_type, ans_countrate):
                              vegaspec=_vspec)
         self._compare_countrate(rn_sp, ans_countrate)
 
-    def test_renorm_noband_count(self):
-        """No bandpass. This option is not offered by ASTROLIB PYSYNPHOT
-        but can be indirectly calculated using a very large box as bandpass.
-
-        """
-        rn_sp = self.em.normalize(1 * u.ct, area=_area)
-        x = rn_sp.integrate(flux_unit=u.ct, area=_area)
-        ans = 10.615454634451927
-        np.testing.assert_allclose(x.value, ans, rtol=1e-3)
-
     def test_renorm_noband_jy(self):
         """Replace this with real test when it is implemented."""
         with pytest.raises(NotImplementedError):