Repository: spacetelescope/synphot_refactor

Issue:
Implement analytic integration for 1D models
Unlikely #88 would happen anytime soon, so just go ahead an implement it here for the models that we use. No point to use trapezoid integration for simple things like `Box1D`.

Discussion:

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index 0e11d85..a373c48 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -14,6 +14,10 @@
 - Added option for ``synphot.utils.download_data()`` to download to the cache
   instead of a specific location. Please note that new option is not fully
   compatible with customization using ``synphot.cfg``. [#211]
+- Added option to use analytic integral for some models. However, for backward
+  compatibility, the default is still trapezoid integration. [#252]
+- Trapezoid integration now provides unsigned area for ``RickerWavelet1D``.
+  [#252]
 
 0.2.1 (2019-12-20)
 ==================
diff --git a/docs/index.rst b/docs/index.rst
index 86901ae..915158c 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -151,15 +151,13 @@ work, so they can be used directly. When in doubt, see if a model is in
 
     >>> from astropy import units as u
     >>> from synphot import units, SourceSpectrum
-    >>> from synphot.models import (BlackBodyNorm1D, GaussianAbsorption1D,
-    ...                             GaussianFlux1D)
-    >>> from synphot.spectrum import BaseUnitlessSpectrum
+    >>> from synphot.models import BlackBodyNorm1D, GaussianFlux1D
 
-Create a Gaussian absorption profile that absorps 80% of the flux at
+Create a Gaussian absorption line with the given amplitude centered at
 4000 Angstrom with a sigma of 20 Angstrom::
 
-    >>> g_abs = BaseUnitlessSpectrum(GaussianAbsorption1D, amplitude=0.8,
-    ...                              mean=4000, stddev=20)
+    >>> g_abs = SourceSpectrum(GaussianFlux1D, amplitude=1*u.mJy,
+    ...                        mean=4000, stddev=20)
 
 Create a Gaussian emission line with the given total flux
 centered at 3000 Angstrom with FWHM of 100 Angstrom::
@@ -175,7 +173,7 @@ Create a blackbody source spectrum with a temperature of 6000 K::
 Combine the above components to create a source spectrum that is twice
 the original blackbody flux with the Gaussian emission and absorption lines::
 
-    >>> sp = g_abs * (g_em + 2 * bb)
+    >>> sp = 2 * bb + g_em - g_abs
 
 Plot the spectrum, zooming in on the line features::
 
@@ -185,30 +183,19 @@ Plot the spectrum, zooming in on the line features::
 
     from astropy import units as u
     from synphot import units, SourceSpectrum
-    from synphot.models import (BlackBodyNorm1D, GaussianAbsorption1D,
-                                GaussianFlux1D)
-    from synphot.spectrum import BaseUnitlessSpectrum
-    # Create a Gaussian absorption profile that absorps 80% of the flux at
-    # 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 the given total flux
-    # centered at 3000 Angstrom with FWHM of 100 Angstrom
+    from synphot.models import BlackBodyNorm1D, GaussianFlux1D
+    g_abs = SourceSpectrum(GaussianFlux1D, amplitude=1*u.mJy,
+                           mean=4000, stddev=20)
     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)
-    # Combine the above components to create a source spectrum that is twice
-    # the original blackbody flux with the Gaussian emission and absorption
-    # lines
-    sp = g_abs * (g_em + 2 * bb)
-    # Plot the spectrum, zooming in on the line features
+    sp = 2 * bb + g_em - g_abs
     sp.plot(left=1, right=7000)
 
 Sample the spectrum at 0.3 micron::
 
-    >>> sp(0.3 * u.micron)  # doctest: +FLOAT_CMP +IGNORE_WARNINGS
+    >>> sp(0.3 * u.micron)  # doctest: +FLOAT_CMP
     <Quantity 0.00129536 PHOTLAM>
 
 Or sample the same thing but in a different flux unit::
@@ -220,7 +207,7 @@ Sample the spectrum at its native wavelength set::
 
     >>> sp(sp.waveset)  # doctest: +FLOAT_CMP +ELLIPSIS
     <Quantity [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
-               4.47483287e-05, 4.34264666e-05, 4.21423394e-05] PHOTLAM>
+               4.47483872e-05, 4.34265235e-05, 4.21423946e-05] PHOTLAM>
 
 Models that built the spectrum::
 
@@ -230,15 +217,15 @@ Models that built the spectrum::
     Inputs: ('x',)
     Outputs: ('y',)
     Model set size: 1
-    Expression: ([0] + ([1] | [2])) * [3]
+    Expression: ([0] | [1]) + [2] - [3]
     Components:
-        [0]: <GaussianFlux1D(...>
+        [0]: <BlackBodyNorm1D(...)>
     <BLANKLINE>
-        [1]: <BlackBodyNorm1D(...)>
+        [1]: <Scale(...)>
     <BLANKLINE>
-        [2]: <Scale(...)>
+        [2]: <GaussianFlux1D(...)>
     <BLANKLINE>
-        [3]: <GaussianAbsorption1D(...)>
+        [3]: <GaussianFlux1D(...)>
     Parameters:
         ...
 
@@ -258,7 +245,7 @@ integrate it::
 
     >>> sp_rn = sp.normalize(1 * u.Jy, band=bp)
     >>> sp_rn.integrate()  # doctest: +FLOAT_CMP
-    <Quantity 12883.09903646 ph / (cm2 s)>
+    <Quantity 12540.4613615 ph / (cm2 s)>
 
 Create an observation by passing the redshifted and normalized source spectrum
 through the box bandpass::
@@ -273,7 +260,7 @@ Calculate the count rate of the observation above for an 2-meter telescope:
     >>> area  # doctest: +FLOAT_CMP
     <Quantity 3.14159265 m2>
     >>> obs.countrate(area=area)  # doctest: +FLOAT_CMP
-    <Quantity 24219650.71772057 ct / s>
+    <Quantity 24219649.46880912 ct / s>
 
 
 .. _synphot_using:
diff --git a/docs/synphot/overview.rst b/docs/synphot/overview.rst
index 4bca019..68d9cc0 100644
--- a/docs/synphot/overview.rst
+++ b/docs/synphot/overview.rst
@@ -97,8 +97,7 @@ and special notes:
 |`~synphot.models.Gaussian1D`                             |**synphot**|Like `astropy.modeling.functional_models.Gaussian1D`          |
 |                                                         |           |but with ``sampleset``.                                       |
 +---------------------------------------------------------+-----------+--------------------------------------------------------------+
-|`~synphot.models.GaussianAbsorption1D`                   |**synphot**|Like `astropy.modeling.functional_models.GaussianAbsorption1D`|
-|                                                         |           |but with ``sampleset``.                                       |
+|`~synphot.models.GaussianAbsorption1D`                   |**synphot**|This is ``1 - Gaussian1D`` model.                             |
 +---------------------------------------------------------+-----------+--------------------------------------------------------------+
 |`~synphot.models.GaussianFlux1D`                         |**synphot**|Like `~synphot.models.Gaussian1D` but allows                  |
 |                                                         |           |backward-compatible parameters like total flux and            |
diff --git a/docs/synphot/spectrum.rst b/docs/synphot/spectrum.rst
index 36a12eb..de3e55f 100644
--- a/docs/synphot/spectrum.rst
+++ b/docs/synphot/spectrum.rst
@@ -188,20 +188,27 @@ necessary to achieve the desired normalization::
    :alt: Renormalize blackbody.
 
 **Integration** is done with the
-:meth:`~synphot.spectrum.BaseSpectrum.integrate`
-method. It uses trapezoid integration (but could be expanded to perform
-analytical calculations instead in the future when that is supported by
-Astropy). By default, integration is done in internal units::
+:meth:`~synphot.spectrum.BaseSpectrum.integrate` method. It can use trapezoid
+or analytical integration (if the latter is available). The type of integration
+being done can be controlled with software configuration or keyword.
+By default, trapezoid integration is done in internal units::
 
     >>> from astropy import units as u
-    >>> from synphot import SourceSpectrum, units
+    >>> from synphot import SourceSpectrum, units, conf
     >>> from synphot.models import GaussianFlux1D
-    >>> sp = SourceSpectrum(GaussianFlux1D, mean=6000, fwhm=10,
-    ...                     total_flux=1*u.erg/(u.cm**2 * u.s))
+    >>> sp = SourceSpectrum(GaussianFlux1D, mean=6000*u.AA, fwhm=10*u.AA,
+    ...                     total_flux=1*(u.erg/(u.cm**2 * u.s)))
     >>> sp.integrate()  # doctest: +FLOAT_CMP
-    <Quantity 3.02046763e+11 ph / (cm2 s)>
+    <Quantity 3.02046758e+11 ph / (cm2 s)>
+    >>> with conf.set_temp('default_integrator', 'analytical'):
+    ...     print(f'{repr(sp.integrate())}')  # doctest: +FLOAT_CMP
+    <Quantity 3.02046994e+11 ph / (cm2 s)>
+    >>> sp.integrate(integration_type='analytical')  # doctest: +FLOAT_CMP
+    <Quantity 3.02046994e+11 ph / (cm2 s)>
     >>> sp.integrate(flux_unit=units.FLAM)  # doctest: +FLOAT_CMP
     <Quantity 0.99999972 erg / (cm2 s)>
+    >>> sp.integrate(flux_unit=units.FLAM, integration_type='analytical')  # doctest: +FLOAT_CMP
+    <Quantity 1. erg / (cm2 s)>
 
 
 .. _synphot-empirical-source:
@@ -343,10 +350,9 @@ The example below creates and plots a flat source with the amplitude of
 Gaussian Absorption
 -------------------
 
-There are two ways to create a Gaussian absorption feature;
-You can choose whichever method that better suits your needs.
-One is to first create :ref:`synphot-gaussian` and then subtract it from
-a continuum (e.g., :ref:`synphot-flat-spec`):
+To create a Gaussian absorption feature, you can first create
+:ref:`synphot-gaussian` and then subtract it from a continuum
+(e.g., :ref:`synphot-flat-spec`):
 
 .. plot::
     :include-source:
@@ -360,22 +366,6 @@ a continuum (e.g., :ref:`synphot-flat-spec`):
     sp = bg - em
     sp.plot()
 
-The other way is to create a unitless absorption profile and then multiply it
-to a continuum (e.g., :ref:`synphot-flat-spec`):
-
-.. plot::
-    :include-source:
-
-    from astropy.stats.funcs import gaussian_fwhm_to_sigma
-    from synphot import SourceSpectrum, BaseUnitlessSpectrum
-    from synphot.models import GaussianAbsorption1D, ConstFlux1D
-    sig = 10 * gaussian_fwhm_to_sigma
-    ab = BaseUnitlessSpectrum(GaussianAbsorption1D, mean=6000, stddev=sig,
-                              amplitude=0.047)
-    bg = SourceSpectrum(ConstFlux1D, amplitude=2)
-    sp = bg * ab
-    sp.plot()
-
 
 .. _synphot-gaussian:
 
diff --git a/synphot/compat.py b/synphot/compat.py
index 4c75f8b..c0a23d1 100644
--- a/synphot/compat.py
+++ b/synphot/compat.py
@@ -19,6 +19,8 @@ else:
     HAS_DUST_EXTINCTION = True
 
 
-__all__ = ['ASTROPY_LT_4_0', 'HAS_SPECUTILS', 'HAS_DUST_EXTINCTION']
+__all__ = ['ASTROPY_LT_4_1', 'ASTROPY_LT_4_0', 'HAS_SPECUTILS',
+           'HAS_DUST_EXTINCTION']
 
+ASTROPY_LT_4_1 = not minversion(astropy, '4.1')
 ASTROPY_LT_4_0 = not minversion(astropy, '4.0')
diff --git a/synphot/config.py b/synphot/config.py
index 7b43416..b2345af 100644
--- a/synphot/config.py
+++ b/synphot/config.py
@@ -16,10 +16,9 @@ __all__ = ['conf']
 class Conf(ConfigNamespace):
     """Configuration parameters."""
 
-    # TODO: Add more options here in the future
+    # First option in list is the default.
     default_integrator = ConfigItem(
-        ['trapezoid'],
-        'Default integrator to use when analytical integration is unavailable')
+        ['trapezoid', 'analytical'], 'Default integrator to use')
 
     # STANDARD STARS
     vega_file = ConfigItem(
diff --git a/synphot/models.py b/synphot/models.py
index f124dd9..fb17386 100644
--- a/synphot/models.py
+++ b/synphot/models.py
@@ -2,6 +2,7 @@
 """Spectrum models not in `astropy.modeling`."""
 
 # STDLIB
+import math
 import warnings
 from copy import deepcopy
 from functools import partial
@@ -134,6 +135,12 @@ class BlackBody1D(Fittable1DModel):
 
         return bbflux.value
 
+    def integrate(self, *args):
+        with u.add_enabled_equivalencies(u.temperature()):
+            t = u.Quantity(self.temperature, u.K)
+
+        return (const.sigma_sb * t ** 4 / math.pi)  # per steradian
+
 
 class BlackBodyNorm1D(BlackBody1D):
     """Create a normalized :ref:`blackbody spectrum <synphot-planck-law>`
@@ -176,6 +183,9 @@ class BlackBodyNorm1D(BlackBody1D):
         bbflux = super(BlackBodyNorm1D, self).evaluate(x, temperature)
         return bbflux * self._omega
 
+    def integrate(self, *args):
+        return super().integrate(*args) * self._omega
+
 
 class Box1D(_models.Box1D):
     """Same as `astropy.modeling.functional_models.Box1D`, except with
@@ -215,6 +225,13 @@ class Box1D(_models.Box1D):
 
         return np.asarray(w)
 
+    def integrate(self, *args):
+        # TODO: Remove unit hardcoding when we use model with units natively.
+        with u.add_enabled_equivalencies(u.spectral()):
+            w = u.Quantity(self.width, u.AA)
+
+        return self.amplitude * w
+
 
 class ConstFlux1D(_models.Const1D):
     """One dimensional constant flux model.
@@ -269,6 +286,19 @@ class ConstFlux1D(_models.Const1D):
         y = units.convert_flux(x, a, units.PHOTLAM)
         return y.value
 
+    def integrate(self, x):
+        # TODO: Remove unit hardcoding when we use model with units natively.
+        # TODO: We do not handle wav_unit as wave number nor energy for now.
+        if 'wav' in self._flux_unit.physical_type:
+            wav_unit = u.AA
+        else:
+            wav_unit = u.Hz
+        with u.add_enabled_equivalencies(u.spectral()):
+            x = u.Quantity(x, wav_unit)
+            amp = u.Quantity(self.amplitude, self._flux_unit)
+
+        return (max(x) - min(x)) * amp
+
 
 class Empirical1D(Tabular1D):
     """Empirical (sampled) spectrum or bandpass model.
@@ -414,6 +444,8 @@ class BaseGaussian1D(_models.Gaussian1D):
     ``sampleset`` defined.
 
     """
+    _sqrt_2_pi = math.sqrt(2 * math.pi)
+
     def sampleset(self, factor_step=0.1, **kwargs):
         """Return ``x`` array that samples the feature.
 
@@ -448,9 +480,16 @@ class Gaussian1D(BaseGaussian1D):
     ``sampleset`` defined.
 
     """
-    pass
+    def integrate(self, *args):
+        # TODO: Remove unit hardcoding when we use model with units natively.
+        with u.add_enabled_equivalencies(u.spectral()):
+            stddev = u.Quantity(self.stddev, u.AA)
 
+        return self.amplitude * stddev * self._sqrt_2_pi
 
+
+# TODO: Deprecate this?
+# This is not really supported anymore but kept for backward compatibility.
 class GaussianAbsorption1D(BaseGaussian1D):
     """Same as ``astropy.modeling.functional_models.GaussianAbsorption1D``,
     except with ``sampleset`` defined.
@@ -509,7 +548,7 @@ class GaussianFlux1D(Gaussian1D):
         else:
             self.stddev = fwhm * gaussian_fwhm_to_sigma
 
-        gaussian_amp_to_totflux = np.sqrt(2.0 * np.pi) * self.stddev
+        gaussian_amp_to_totflux = self._sqrt_2_pi * self.stddev
 
         if total_flux is None:
             u_str = 'PHOTLAM'
@@ -529,6 +568,10 @@ class GaussianFlux1D(Gaussian1D):
         self.meta['expr'] = 'em({0:g}, {1:g}, {2:g}, {3})'.format(
             self.mean.value, fwhm, total_flux, u_str)
 
+    def integrate(self, *args):
+        # TODO: Remove unit hardcoding when we use model with units natively.
+        return super(GaussianFlux1D, self).integrate(*args) * units.PHOTLAM
+
 
 class Lorentz1D(_models.Lorentz1D):
     """Same as `astropy.modeling.functional_models.Lorentz1D`, except with
@@ -558,6 +601,18 @@ class Lorentz1D(_models.Lorentz1D):
 
         return np.asarray(w)
 
+    def integrate(self, x):
+        # TODO: Remove unit hardcoding when we use model with units natively.
+        with u.add_enabled_equivalencies(u.spectral()):
+            x = u.Quantity(x, u.AA)
+            x_0 = u.Quantity(self.x_0, u.AA)
+            gamma = u.Quantity(self.fwhm, u.AA) * 0.5
+
+        a1 = np.arctan((min(x) - x_0) / gamma)
+        a2 = np.arctan((max(x) - x_0) / gamma)
+        da = (a2 - a1).to(u.dimensionless_unscaled, u.dimensionless_angles())
+        return self.amplitude * gamma * da
+
 
 class RickerWavelet1D(_RickerWavelet1D):
     """Same as `astropy.modeling.functional_models.RickerWavelet1D`, except
@@ -587,6 +642,37 @@ class RickerWavelet1D(_RickerWavelet1D):
 
         return np.asarray(w)
 
+    def integrate(self, x):
+        # TODO: Remove unit hardcoding when we use model with units natively.
+        with u.add_enabled_equivalencies(u.spectral()):
+            x = u.Quantity(x, u.AA)
+            x_0 = u.Quantity(self.x_0, u.AA)
+            sig = u.Quantity(self.sigma, u.AA)
+
+        # Roots, where y=0
+        root_left = x_0 - sig
+        root_right = x_0 + sig
+
+        x_min = min(x)
+        x_max = max(x)
+        if x_min >= root_left or x_max <= root_right:
+            raise NotImplementedError(
+                'Partial analytic integration not supported')
+
+        sig2 = sig * sig
+
+        def _int_subregion(xx1, xx2):
+            dx_min = xx1 - x_0
+            dx_max = xx2 - x_0
+            a1 = dx_min * np.exp(-0.5 * dx_min * dx_min / sig2)
+            a2 = dx_max * np.exp(-0.5 * dx_max * dx_max / sig2)
+            return abs(a2 - a1)
+
+        # Unsigned area
+        return self.amplitude * (_int_subregion(x_min, root_left) +
+                                 _int_subregion(root_left, root_right) +
+                                 _int_subregion(root_right, x_max))
+
 
 # TODO: Emit proper deprecation warning.
 # https://github.com/spacetelescope/synphot_refactor/issues/249
@@ -642,6 +728,17 @@ class PowerLawFlux1D(_models.PowerLaw1D):
         flux = units.convert_flux(x, y, units.PHOTLAM)
         return flux.value
 
+    def integrate(self, x):
+        # TODO: Remove unit hardcoding when we use model with units natively.
+        with u.add_enabled_equivalencies(u.spectral()):
+            x = u.Quantity(x, u.AA)
+            x_0 = u.Quantity(self.x_0, u.AA)
+            amp = u.Quantity(self.amplitude, self._flux_unit)
+
+        fac = 1 - self.alpha
+        denom = x_0 ** -self.alpha * fac
+        return amp * (max(x) ** fac - min(x) ** fac) / denom
+
 
 class Trapezoid1D(_models.Trapezoid1D):
     """Same as `astropy.modeling.functional_models.Trapezoid1D`, except with
@@ -662,6 +759,14 @@ class Trapezoid1D(_models.Trapezoid1D):
 
         return np.asarray(w)
 
+    def integrate(self, *args):
+        # TODO: Remove unit hardcoding when we use model with units natively.
+        with u.add_enabled_equivalencies(u.spectral()):
+            width = u.Quantity(self.width, u.AA)
+            slope = u.Quantity(self.slope, 1 / u.AA)
+
+        return self.amplitude * (width + self.amplitude / slope)
+
 
 # Functions below are for sampleset magic.
 
diff --git a/synphot/observation.py b/synphot/observation.py
index 3cabd36..c108661 100644
--- a/synphot/observation.py
+++ b/synphot/observation.py
@@ -471,7 +471,8 @@ class Observation(BaseSourceSpectrum):
             if not isinstance(vegaspec, SourceSpectrum):
                 raise exceptions.SynphotError('Vega spectrum is missing.')
             num = self.integrate(wavelengths=wavelengths)
-            den = (vegaspec * self.bandpass).integrate()
+            den = (vegaspec * self.bandpass).integrate(
+                integration_type='trapezoid')
             utils.validate_totalflux(num)
             utils.validate_totalflux(den)
             return (2.5 * (math.log10(den.value) -
diff --git a/synphot/spectrum.py b/synphot/spectrum.py
index fcb464c..28d748c 100644
--- a/synphot/spectrum.py
+++ b/synphot/spectrum.py
@@ -416,30 +416,49 @@ class BaseSpectrum(object):
         """
         return self.__truediv__(other)
 
-    def integrate(self, wavelengths=None, **kwargs):
+    def integrate(self, wavelengths=None, integration_type=None, **kwargs):
         """Perform integration.
 
-        This uses any analytical integral that the
-        underlying model has (i.e., ``self.model.integral``).
-        If unavailable, it uses the default fall-back integrator
-        set in the ``default_integrator`` configuration item.
-
-        If wavelengths are provided, flux or throughput is first resampled.
+        When integration is not analytical and wavelengths are provided,
+        flux or throughput is first resampled.
         This is useful when user wants to integrate at specific end points
         or use custom spacing; In that case, user can pass in desired
         sampling array generated with :func:`numpy.linspace`,
         :func:`numpy.logspace`, etc.
         If not provided, then `waveset` is used.
 
+        When integration is analytical, wavelengths provided are only used
+        if applicable for the particular model. If the model does not
+        support analytical integration, it will fall back to simple
+        trapezoid integration.
+
         Parameters
         ----------
         wavelengths : array-like, `~astropy.units.quantity.Quantity`, or `None`
             Wavelength values for integration.
+            This is ignored by analytical integration if not applicable.
             If not a Quantity, assumed to be in Angstrom.
             If `None`, `waveset` is used.
 
+        integration_type : {None, 'trapezoid', 'analytical'}
+            Defines how the integration is done, either by simple trapezoid
+            integration or analytical formula. If `None`, the option is pulled
+            from ``synphot.config.conf.default_integrator``. If analytical
+            integration is requested but no possible, trapezoid integration
+            is done anyway.
+
+        flux_unit : str, `~astropy.units.core.Unit`, or `None`
+            **This option is only available for source spectrum.**
+            For trapezoid integration, flux is converted to this unit for
+            sampling before integration. For analytical integration, while
+            calculation is done differently, providing this option would
+            result in equivalent behavior as trapezoid integration, where
+            possible, for consistency. If not given, internal unit is used.
+
         kwargs : dict
-            Optional keywords to ``__call__`` for sampling.
+            **This option is only available for source spectrum.**
+            Other optional keywords besides ``flux_unit`` to ``__call__``
+            for sampling when integration type is not analytical.
 
         Returns
         -------
@@ -449,57 +468,77 @@ class BaseSpectrum(object):
         Raises
         ------
         NotImplementedError
-            Invalid default integrator.
+            Invalid integration type.
 
         synphot.exceptions.SynphotError
             `waveset` is needed but undefined or cannot integrate
             natively in the given ``flux_unit``.
 
         """
+        # For non-analytical integration:
         # 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)
+        # For analytical integration: Keep the same behavior for consistency.
+        flux_unit = kwargs.get('flux_unit')
+        is_unitless = self._internal_flux_unit == units.THROUGHPUT
+        if flux_unit is not None:
+            if is_unitless:
+                raise exceptions.SynphotError(
+                    'flux_unit cannot be used with unitless spectrum')
+            else:
+                self._validate_flux_unit(flux_unit, wav_only=True)
 
         x = self._validate_wavelengths(wavelengths)
 
-        # TODO: When astropy.modeling.models supports this, need to
-        #       make sure that this actually works, and gives correct unit.
-        # https://github.com/astropy/astropy/issues/5033
-        # https://github.com/astropy/astropy/pull/5108
-        try:
-            m = self.model.integral
-        except (AttributeError, NotImplementedError):
-            if conf.default_integrator == 'trapezoid':
-                y = self(x, **kwargs)
-                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))
-            result_unit = self._internal_flux_unit
+        if integration_type is None:
+            integration_type = conf.default_integrator
+
+        # NOTE: Emitting warning is not done because it overcomplicates
+        #       the logic here.
+        # Fallback to trapezoid if analytical not possible.
+        if (integration_type == 'analytical' and
+                not hasattr(self.model, 'integrate')):
+            integration_type = 'trapezoid'
+
+        if integration_type == 'trapezoid':
+            y = abs(self(x, **kwargs))  # Unsigned area
+            result = abs(np.trapz(y, x=x))
+
+        elif integration_type == 'analytical':
+            result = self.model.integrate(x)
+
+            # TODO: Remove unit hardcoding when we use model with units
+            #       natively.
+            if not is_unitless and result.unit.physical_type == 'length':
+                result = result * self._internal_flux_unit
+
+            # NOTE: flux_unit is flux density, not integrated.
+            # Use wavelength for unit conversion, if applicable.
+            if not is_unitless and flux_unit is not None:
+                modelname = self.model.__class__.__name__
+                if modelname in self._model_fconv_wav:
+                    pname_wav = self._model_fconv_wav[modelname]
+                    wav = getattr(self.model, pname_wav).value
+                    with u.add_enabled_equivalencies(u.spectral()):
+                        wav = u.Quantity(wav, self._internal_wave_unit)
+                    to_unit = flux_unit * self._internal_wave_unit
+                    result = units.convert_flux(wav, result, to_unit)
+        else:
+            raise NotImplementedError(
+                '{} is not a supported integration '
+                'type'.format(integration_type))
 
         # 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
+        # and make it pretty, where applicable.
+        if not is_unitless:
+            result_unit_str = result.unit.to_string()
+            if 'ph' in result_unit_str or 'PHOTLAM' in result_unit_str:
+                result = result.to(u.photon / (u.cm**2 * u.s))
+            else:  # FLAM
+                result = result.to(u.erg / (u.cm**2 * u.s))
 
-        return result * result_unit
+        return result
 
     def avgwave(self, wavelengths=None):
         """Calculate the :ref:`average wavelength <synphot-formula-avgwv>`.
@@ -824,7 +863,7 @@ class BaseSpectrum(object):
             If not a Quantity, assumed to be in Angstrom.
             If `None`, ``self.waveset`` is used.
 
-        flux_unit : str or `~astropy.units.core.Unit` or `None`
+        flux_unit : str, `~astropy.units.core.Unit`, or `None`
             This option is not applicable to unitless spectrum like bandpass.
             Flux is converted to this unit before written out.
             If not given, internal unit is used.
@@ -889,7 +928,7 @@ class BaseSourceSpectrum(BaseSpectrum):
             Wavelength values for sampling. If not a Quantity,
             assumed to be in Angstrom.
 
-        flux_unit : str or `~astropy.units.core.Unit` or `None`
+        flux_unit : str, `~astropy.units.core.Unit`, or `None`
             Flux is converted to this unit.
             If not given, internal unit is used.
 
@@ -1029,7 +1068,8 @@ class BaseSourceSpectrum(BaseSpectrum):
             totalflux = flux_tmp.sum().value
             stdflux = 1.0
         else:
-            totalflux = sp.integrate(wavelengths=wavelengths).value
+            totalflux = sp.integrate(wavelengths=wavelengths,
+                                     integration_type='trapezoid')
 
             # VEGAMAG
             if renorm_unit_name == units.VEGAMAG.to_string():
@@ -1055,7 +1095,8 @@ class BaseSourceSpectrum(BaseSpectrum):
                 raise NotImplementedError('Must provide a bandpass')
             else:
                 up = stdspec * band
-                stdflux = up.integrate(wavelengths=wavelengths).value
+                stdflux = up.integrate(wavelengths=wavelengths,
+                                       integration_type='trapezoid')
 
         utils.validate_totalflux(totalflux)
 
@@ -1230,7 +1271,7 @@ class SourceSpectrum(BaseSourceSpectrum):
             If not a Quantity, assumed to be in Angstrom.
             If `None`, ``self.waveset`` is used.
 
-        flux_unit : str or `~astropy.units.core.Unit` or `None`
+        flux_unit : str, `~astropy.units.core.Unit`, or `None`
             Flux is converted to this unit for plotting.
             If not given, internal unit is used.
 
@@ -1264,7 +1305,7 @@ class SourceSpectrum(BaseSourceSpectrum):
             If not a Quantity, assumed to be in Angstrom.
             If `None`, ``self.waveset`` is used.
 
-        flux_unit : str or `~astropy.units.core.Unit` or `None`
+        flux_unit : str, `~astropy.units.core.Unit`, or `None`
             Flux is converted to this unit before written out.
             If not given, internal unit is used.
 
@@ -1483,9 +1524,9 @@ class SpectralElement(BaseUnitlessSpectrum):
 
         x1 = self._validate_wavelengths(wavelengths)
         y1 = self(x1)
-        a = x1[y1 > 0].value
+        a = x1[y1 > 0]
 
-        b = other._validate_wavelengths(wavelengths).value
+        b = other._validate_wavelengths(wavelengths)
 
         result = utils.overlap_status(a, b)
 
@@ -1504,20 +1545,23 @@ class SpectralElement(BaseUnitlessSpectrum):
             # Check if the lack of overlap is significant.
             else:
                 # Get all the flux
-                totalflux = self.integrate(wavelengths=wavelengths).value
+                totalflux = self.integrate(wavelengths=wavelengths,
+                                           integration_type='trapezoid')
                 utils.validate_totalflux(totalflux)
 
                 a_min, a_max = a.min(), a.max()
                 b_min, b_max = b.min(), b.max()
 
                 # Now get the other two pieces
-                excluded = 0.0
+                excluded = 0.0 * totalflux.unit
                 if a_min < b_min:
                     excluded += self.integrate(
-                        wavelengths=np.array([a_min, b_min])).value
+                        wavelengths=u.Quantity([a_min, b_min]),
+                        integration_type='trapezoid')
                 if a_max > b_max:
                     excluded += self.integrate(
-                        wavelengths=np.array([b_max, a_max])).value
+                        wavelengths=u.Quantity([b_max, a_max]),
+                        integration_type='trapezoid')
 
                 if excluded / totalflux < threshold:
                     result = 'partial_most'
@@ -1551,13 +1595,13 @@ class SpectralElement(BaseUnitlessSpectrum):
         a = units.validate_quantity(area, units.AREA)
 
         # Only correct if wavelengths are in Angstrom.
-        x = self._validate_wavelengths(wavelengths).value
+        x = self._validate_wavelengths(wavelengths).to(u.AA)
 
-        y = self(x).value * x
+        y = self(x) * x
         int_val = abs(np.trapz(y, x=x))
         uresp = units.HC / (a.cgs * int_val)
 
-        return uresp.value * units.FLAM
+        return (uresp / u.s).to(units.FLAM)
 
     def rmswidth(self, wavelengths=None, threshold=None):
         """Calculate the :ref:`bandpass RMS width <synphot-formula-rmswidth>`.
@@ -1586,33 +1630,32 @@ class SpectralElement(BaseUnitlessSpectrum):
             Threshold is invalid.
 
         """
-        x = self._validate_wavelengths(wavelengths).value
-        y = self(x).value
+        x = self._validate_wavelengths(wavelengths)
+        y = self(x)
 
         if threshold is None:
             wave = x
             thru = y
         else:
-            if (isinstance(threshold, numbers.Real) or
-                (isinstance(threshold, u.Quantity) and
-                 threshold.unit == self._internal_flux_unit)):
+            try:
                 mask = y >= threshold
-            else:
+            except Exception as e:
                 raise exceptions.SynphotError(
-                    '{0} is not a valid threshold'.format(threshold))
+                    '{0} is not a valid threshold: '
+                    '{1}'.format(threshold, str(e)))
             wave = x[mask]
             thru = y[mask]
 
-        a = self.avgwave(wavelengths=wavelengths).value
+        a = self.avgwave(wavelengths=wavelengths)
         num = np.trapz((wave - a) ** 2 * thru, x=wave)
         den = np.trapz(thru, x=wave)
 
         if den == 0:  # pragma: no cover
-            rms_width = 0.0
+            rms_width = 0.0 * a.unit
         else:
             rms_width = np.sqrt(abs(num / den))
 
-        return rms_width * self._internal_wave_unit
+        return rms_width
 
     def photbw(self, wavelengths=None, threshold=None):
         """Calculate the
@@ -1644,37 +1687,36 @@ class SpectralElement(BaseUnitlessSpectrum):
             Threshold is invalid.
 
         """
-        x = self._validate_wavelengths(wavelengths).value
-        y = self(x).value
+        x = self._validate_wavelengths(wavelengths)
+        y = self(x)
 
         if threshold is None:
             wave = x
             thru = y
         else:
-            if (isinstance(threshold, numbers.Real) or
-                (isinstance(threshold, u.Quantity) and
-                 threshold.unit == self._internal_flux_unit)):
+            try:
                 mask = y >= threshold
-            else:
+            except Exception as e:
                 raise exceptions.SynphotError(
-                    '{0} is not a valid threshold'.format(threshold))
+                    '{0} is not a valid threshold: '
+                    '{1}'.format(threshold, str(e)))
             wave = x[mask]
             thru = y[mask]
 
-        a = self.barlam(wavelengths=wavelengths).value
+        a = self.barlam(wavelengths=wavelengths)
 
         if a == 0:
-            bandw = 0.0
+            bandw = 0.0 * a.unit
         else:
             num = np.trapz(thru * np.log(wave / a) ** 2 / wave, x=wave)
             den = np.trapz(thru / wave, x=wave)
 
             if den == 0:  # pragma: no cover
-                bandw = 0.0
+                bandw = 0.0 * a.unit
             else:
                 bandw = a * np.sqrt(abs(num / den))
 
-        return bandw * self._internal_wave_unit
+        return bandw
 
     def fwhm(self, **kwargs):
         """Calculate :ref:`synphot-formula-fwhm` of equivalent gaussian.
@@ -1725,7 +1767,7 @@ class SpectralElement(BaseUnitlessSpectrum):
             Peak bandpass throughput.
 
         """
-        x = self._validate_wavelengths(wavelengths).value
+        x = self._validate_wavelengths(wavelengths)
         return self(x).max()
 
     def wpeak(self, wavelengths=None):
@@ -1751,15 +1793,13 @@ class SpectralElement(BaseUnitlessSpectrum):
         x = self._validate_wavelengths(wavelengths)
         return x[self(x) == self.tpeak(wavelengths=wavelengths)][0]
 
-    def equivwidth(self, wavelengths=None):
+    def equivwidth(self, **kwargs):
         """Calculate :ref:`bandpass equivalent width <synphot-formula-equvw>`.
 
         Parameters
         ----------
-        wavelengths : array-like, `~astropy.units.quantity.Quantity`, or `None`
-            Wavelength values for sampling.
-            If not a Quantity, assumed to be in Angstrom.
-            If `None`, ``self.waveset`` is used.
+        kwargs : dict
+            See :meth:`integrate`.
 
         Returns
         -------
@@ -1767,7 +1807,7 @@ class SpectralElement(BaseUnitlessSpectrum):
             Bandpass equivalent width.
 
         """
-        return self.integrate(wavelengths=wavelengths)
+        return self.integrate(**kwargs)
 
     def rectwidth(self, wavelengths=None):
         """Calculate :ref:`bandpass rectangular width <synphot-formula-rectw>`.
@@ -1785,7 +1825,8 @@ class SpectralElement(BaseUnitlessSpectrum):
             Bandpass rectangular width.
 
         """
-        equvw = self.equivwidth(wavelengths=wavelengths)
+        equvw = self.equivwidth(wavelengths=wavelengths,
+                                integration_type='trapezoid')
         tpeak = self.tpeak(wavelengths=wavelengths)
 
         if tpeak.value == 0:  # pragma: no cover
@@ -1811,10 +1852,9 @@ class SpectralElement(BaseUnitlessSpectrum):
             Dimensionless efficiency.
 
         """
-        x = self._validate_wavelengths(wavelengths).value
-        y = self(x).value
-        qtlam = abs(np.trapz(y / x, x=x))
-        return qtlam * u.dimensionless_unscaled
+        x = self._validate_wavelengths(wavelengths)
+        y = self(x)
+        return abs(np.trapz(y / x, x=x))
 
     def emflx(self, area, wavelengths=None):
         """Calculate
@@ -1837,8 +1877,9 @@ class SpectralElement(BaseUnitlessSpectrum):
             em_flux = 0.0 * units.FLAM
         else:
             uresp = self.unit_response(area, wavelengths=wavelengths)
-            equvw = self.equivwidth(wavelengths=wavelengths).value
-            em_flux = uresp * equvw / t_lambda
+            equvw = self.equivwidth(wavelengths=wavelengths,
+                                    integration_type='trapezoid')
+            em_flux = uresp * equvw / (t_lambda * self._internal_wave_unit)
 
         return em_flux
 
diff --git a/synphot/units.py b/synphot/units.py
index 6ceca5c..555b7a5 100644
--- a/synphot/units.py
+++ b/synphot/units.py
@@ -10,6 +10,7 @@ from astropy import units as u
 
 # LOCAL
 from . import exceptions
+from .compat import ASTROPY_LT_4_1
 
 __all__ = ['H', 'C', 'HC', 'SR_PER_ARCSEC2', 'AREA', 'THROUGHPUT', 'PHOTLAM',
            'PHOTNU', 'FLAM', 'FNU', 'OBMAG', 'VEGAMAG',
@@ -61,6 +62,31 @@ u.add_enabled_units([PHOTLAM, PHOTNU, FLAM, FNU, OBMAG, VEGAMAG])
 # Flux conversion #
 # --------------- #
 
+
+# ASTROPY_LT_4_1: Remove this and just use spectral_density() when astropy
+#                 minversion is 4.1.
+def spectral_density_integrated(wav):
+    """Flux equivalencies for integrated flux.
+
+    .. note::
+
+        For ``astropy>=4.1``, just use
+        :func:`~astropy.units.equivalencies.spectral_density`.
+
+    """
+    la_f_la = u.erg / u.cm ** 2 / u.s
+    la_phot_f_la = u.photon / (u.cm ** 2 * u.s)
+
+    def converter_phot_f_la_to_f_la(x):
+        return HC * x / wav.to_value(u.AA, u.spectral())
+
+    def iconverter_phot_f_la_to_f_la(x):
+        return x * wav.to_value(u.AA, u.spectral()) / HC
+
+    return [(la_phot_f_la, la_f_la,
+             converter_phot_f_la_to_f_la, iconverter_phot_f_la_to_f_la)]
+
+
 def spectral_density_vega(wav, vegaflux):
     """Flux equivalencies between PHOTLAM and VEGAMAG.
 
@@ -196,7 +222,11 @@ def convert_flux(wavelengths, fluxes, out_flux_unit, **kwargs):
     if not isinstance(wavelengths, u.Quantity):
         wavelengths = wavelengths * u.AA
 
-    eqv = u.spectral_density(wavelengths)
+    if ASTROPY_LT_4_1:
+        eqv = (u.spectral_density(wavelengths) +
+               spectral_density_integrated(wavelengths))
+    else:
+        eqv = u.spectral_density(wavelengths)
 
     # Use built-in astropy equivalencies
     try:


Test patch:
diff --git a/synphot/tests/test_integrator.py b/synphot/tests/test_integrator.py
index 856f266..82f5f5e 100644
--- a/synphot/tests/test_integrator.py
+++ b/synphot/tests/test_integrator.py
@@ -1,41 +1,188 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
-"""Test different integration methods."""
+"""Test different integration methods.
+Things not covered here but somewhere else:
+
+* ``conf.set_temp('default_integrator', 'analytical')`` and its default
+  covered in ``spectrum.rst``.
+* The following integrations are tested in ``test_spectrum.py``:
+    * ``BlackBodyNorm1D`` source, and indirectly, ``BlackBody1D``.
+    * ``Box1D`` bandpass.
+    * ``Empirical1D`` source and bandpass.
+    * ``GaussianFlux1D`` source.
+    * ``PowerLawFlux1D`` source.
+    * Compound source and bandpass.
+
+"""
 
 # THIRD-PARTY
 import numpy as np
-from astropy.modeling.models import Polynomial1D
+import pytest
+from astropy import units as u
+from astropy.modeling.models import Const1D
+from astropy.tests.helper import assert_quantity_allclose
 from numpy.testing import assert_allclose
 
 # LOCAL
-from ..models import Box1D
-from ..spectrum import SourceSpectrum
+from .. import models, units
+from ..exceptions import SynphotError
+from ..spectrum import SourceSpectrum, SpectralElement
+
+
+class TestSourceConstFlux1D:
+    def setup_class(self):
+        self.sp = SourceSpectrum(models.ConstFlux1D, amplitude=1 * units.FLAM)
+        self.w = np.arange(5499, 5600) * u.AA
+        self.ans_flam = 100 * (u.erg / (u.cm * u.cm * u.s))
+        self.ans_photlam = 2.79343128e+13 * (u.ph / (u.cm * u.cm * u.s))
+
+    def test_integrate(self):
+        assert_quantity_allclose(
+            self.sp.integrate(wavelengths=self.w), self.ans_photlam)
+        assert_quantity_allclose(
+            self.sp.integrate(wavelengths=self.w, flux_unit='flam'),
+            self.ans_flam)
+        assert_quantity_allclose(
+            self.sp.integrate(wavelengths=self.w, flux_unit='flam',
+                              integration_type='analytical'),
+            self.ans_flam)
+
+        with pytest.raises(SynphotError, match='cannot operate in Jy'):
+            self.sp.integrate(wavelengths=self.w, flux_unit=u.Jy)
+
+    @pytest.mark.xfail(reason='Cannot convert to PHOTLAM')
+    def test_integrate_wontfix(self):
+        """Constant spectrum in FLAM is not constant in PHOTLAM, so cannot
+        apply simple unit conversion to integrated flux.
+        """
+        assert_quantity_allclose(
+            self.sp.integrate(wavelengths=self.w,
+                              integration_type='analytical'),
+            self.ans_photlam)
+
+    def test_integrate_freq_space(self):
+        """When flux is constant in frequency space."""
+        sp = SourceSpectrum(models.ConstFlux1D, amplitude=1 * u.Jy)
+        nu = self.w.to(u.Hz, u.spectral())
+
+        # Answer in this setup is basically
+        # ((max(nu) - min(nu)) * u.Jy).to(u.erg / (u.cm * u.cm * u.s))
+        ans_flam = 9.73703017e-11 * (u.erg / (u.cm * u.cm * u.s))
+
+        assert_quantity_allclose(
+            sp.integrate(wavelengths=self.w, flux_unit='flam'), ans_flam)
+        assert_quantity_allclose(
+            sp.integrate(wavelengths=self.w, flux_unit='flam',
+                         integration_type='analytical'),
+            ans_flam)
+        assert_quantity_allclose(
+            sp.integrate(wavelengths=nu, flux_unit='flam'), ans_flam)
+        assert_quantity_allclose(
+            sp.integrate(wavelengths=nu, flux_unit='flam',
+                         integration_type='analytical'),
+            ans_flam)
+
+
+def test_bandpass_Const1D():
+    """Test Const1D bandpass and a few optional logic routes."""
+    bp = SpectralElement(Const1D, amplitude=1)
+    w = np.arange(5499, 5600) * u.AA
+    ans = 100 * u.AA
+
+    assert_quantity_allclose(bp.integrate(wavelengths=w), ans)
+
+    # Const1D from astropy has no integrate method, so this silently falls
+    # back to trapezoid integration anyway.
+    assert_quantity_allclose(
+        bp.integrate(wavelengths=w, integration_type='analytical'), ans)
+
+    with pytest.raises(SynphotError, match='waveset is undefined'):
+        bp.integrate()
+
+    with pytest.raises(SynphotError, match='flux_unit cannot be used'):
+        bp.integrate(wavelengths=w, flux_unit='flam')
+
+    with pytest.raises(NotImplementedError,
+                       match='not a supported integration type'):
+        bp.integrate(wavelengths=w, integration_type='notreal')
+
+
+def test_bandpass_Gaussian1D():
+    bp = SpectralElement(
+        models.Gaussian1D, amplitude=0.8, mean=5000, stddev=130)
+    ans = 260.68913672 * u.AA
+    assert_quantity_allclose(bp.integrate(), ans)
+    assert_quantity_allclose(bp.integrate(integration_type='analytical'), ans,
+                             rtol=1e-6)
+
+    # Equivalent width is basically the integral.
+    assert_quantity_allclose(bp.equivwidth(), ans)
+    assert_quantity_allclose(bp.equivwidth(integration_type='analytical'), ans,
+                             rtol=1e-6)
+
+
+def test_bandpass_Lorentz1D():
+    bp = SpectralElement(models.Lorentz1D, amplitude=0.8, x_0=5000, fwhm=130)
+    ans = 161.28101053 * u.AA
+    assert_quantity_allclose(bp.integrate(), ans)
+    assert_quantity_allclose(bp.integrate(integration_type='analytical'), ans)
+
+
+def test_source_Lorentz1D():
+    sp = SourceSpectrum(
+        models.Lorentz1D, amplitude=1 * units.PHOTLAM, x_0=5000, fwhm=130)
+    ans_photlam = 201.60126316 * (u.ph / (u.cm * u.cm * u.s))
+    ans_flam = 8.06013348e-10 * (u.erg / (u.cm * u.cm * u.s))
+
+    assert_quantity_allclose(sp.integrate(), ans_photlam)
+    assert_quantity_allclose(sp.integrate(flux_unit='flam'), ans_flam)
+    assert_quantity_allclose(
+        sp.integrate(integration_type='analytical'), ans_photlam)
+
+    # Flux unit conversion is iffy here because we do not have proper
+    # LorentzFlux1D model.
+    assert_quantity_allclose(
+        sp.integrate(flux_unit='flam', integration_type='analytical'),
+        ans_flam, rtol=0.01)
+
+
+def test_bandpass_RickerWavelet1D():
+    """Does not make much sense when transmission is negative, so this
+    makes more sense as equivalent width calculation, which is also the
+    integral.
+    """
+    bp = SpectralElement(
+        models.RickerWavelet1D, amplitude=1, x_0=5000, sigma=100)
+    ans = 242.20770777 * u.AA
+    assert_quantity_allclose(bp.equivwidth(), ans)
+    assert_quantity_allclose(bp.equivwidth(integration_type='analytical'),
+                             ans, rtol=5e-3)
 
+    with pytest.raises(NotImplementedError, match='Partial analytic'):
+        bp.equivwidth(integration_type='analytical',
+                      wavelengths=np.arange(4900, 5200) * u.AA)
 
-# TODO: Replace this with model with actual analytic integral when it
-#       exists in Astropy.
-class Const1D(Box1D):
-    """Fake source with its own analytic integral."""
-    @property
-    def integral(self):
-        return Polynomial1D(1, c1=1)  # f(x) = x
+    with pytest.raises(NotImplementedError, match='Partial analytic'):
+        bp.equivwidth(integration_type='analytical',
+                      wavelengths=np.arange(4800, 5100) * u.AA)
 
 
-# TODO: Update test when Astropy models have integral.
-def test_analytic_source():
-    """Test integration using built-in analytic integral."""
-    sp = SourceSpectrum(Const1D)
-    assert_allclose(sp.integrate([1000, 1010, 1020]).value, 20)
+def test_source_RickerWavelet1D():
+    pytest.xfail(
+        'Not sure if this test makes sense as negative flux is unphysical')
 
 
-# TODO: Replace model with one that is without analytic integral when
-#       integral is in Astropy.
-def test_trapz_spec():
-    sp = SourceSpectrum(Box1D, amplitude=1, x_0=5000, width=10)
-    assert_allclose(sp.integrate().value, 10)
+def test_bandpass_Trapezoid1D():
+    """Trapezoid integration of a trapezoid. How meta!"""
+    bp = SpectralElement(
+        models.Trapezoid1D, amplitude=0.8, x_0=5000, width=10, slope=0.25)
+    ans = 10.56 * u.AA
+    assert_quantity_allclose(bp.integrate(), ans)
+    assert_quantity_allclose(bp.integrate(integration_type='analytical'), ans)
 
 
 def test_trapz_box1d():
-    m = Box1D(amplitude=1, x_0=5000, width=10)
+    """Test the underlying trapezoid integration."""
+    m = models.Box1D(amplitude=1, x_0=5000, width=10)
 
     # Ascending.
     x = m.sampleset()
diff --git a/synphot/tests/test_spectrum.py b/synphot/tests/test_spectrum.py
index 34f9a54..c6fad75 100644
--- a/synphot/tests/test_spectrum.py
+++ b/synphot/tests/test_spectrum.py
@@ -87,14 +87,12 @@ def test_flux_conversion_vega(in_q, out_u, ans):
 
     """
     result = units.convert_flux(_wave, in_q, out_u, vegaspec=_vspec)
-    np.testing.assert_allclose(result.value, ans.value, rtol=1e-2)
-    assert result.unit == ans.unit
+    assert_quantity_allclose(result, ans, rtol=1e-2)
 
     # Scalar
     i = 0
     result = units.convert_flux(_wave[i], in_q[i], out_u, vegaspec=_vspec)
-    np.testing.assert_allclose(result.value, ans[i].value, rtol=1e-2)
-    assert result.unit == ans[i].unit
+    assert_quantity_allclose(result, ans[i], rtol=1e-2)
 
 
 @pytest.mark.remote_data
@@ -150,17 +148,17 @@ class TestEmpiricalSourceFromFile:
         assert self.sp.meta['header']['SIMPLE']  # From FITS header
         assert self.sp.warnings == {}
         assert self.sp.z == 0
-        np.testing.assert_allclose(
-            self.sp.waverange.value, [3479.99902344, 10500.00097656])
+        assert_quantity_allclose(
+            self.sp.waverange, [3479.99902344, 10500.00097656] * u.AA)
 
     def test_call(self):
         w = self.sp.model.points[0][5000:5004]
-        y = units.convert_flux(w, self.sp(w), units.FLAM)
+        y = self.sp(w, flux_unit=units.FLAM)
+        y_ans = [1.87284130e-15, 1.85656811e-15, 1.84030867e-15,
+                 1.82404183e-15] * units.FLAM
         np.testing.assert_allclose(
             w, [6045.1640625, 6045.83203125, 6046.49951172, 6047.16748047])
-        np.testing.assert_allclose(
-            y.value,
-            [1.87284130e-15, 1.85656811e-15, 1.84030867e-15, 1.82404183e-15])
+        assert_quantity_allclose(y, y_ans)
 
     def test_neg_flux(self):
         w = [1000, 5000, 9000]
@@ -174,20 +172,21 @@ class TestEmpiricalSourceFromFile:
     def test_conversion(self):
         x = 0.60451641 * u.micron
         w, y = self.sp._get_arrays(x, flux_unit=units.FNU)
-        np.testing.assert_allclose(x.value, w.value)
-        np.testing.assert_allclose(y.value, 2.282950185743497e-26, rtol=1e-6)
+        assert_quantity_allclose(x, w)
+        assert_quantity_allclose(y, 2.282950185743497e-26 * units.FNU,
+                                 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
+        assert_quantity_allclose(f, 8.460125829057308e-12 * expected_unit,
+                                 rtol=1e-5)
 
         # 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
+        assert_quantity_allclose(f, 4.810058069909525e-14 * expected_unit,
+                                 rtol=1e-5)
 
         # Unsupported unit
         with pytest.raises(exceptions.SynphotError):
@@ -203,10 +202,11 @@ class TestEmpiricalSourceFromFile:
             Empirical1D, points=_wave, lookup_table=_flux_photlam)
         sp = sp2.taper()
         x, y = sp._get_arrays(None, flux_unit=units.FLAM)
-        np.testing.assert_allclose(
-            x.value, [4954.05152484, 4956.8, 4959.55, 4962.3, 4965.05152484])
-        np.testing.assert_allclose(
-            y.value, [0, 3.9135e-14, 4.0209e-14, 3.9169e-14, 0], rtol=1e-6)
+        assert_quantity_allclose(
+            x, [4954.05152484, 4956.8, 4959.55, 4962.3, 4965.05152484] * u.AA)
+        assert_quantity_allclose(
+            y,
+            [0, 3.9135e-14, 4.0209e-14, 3.9169e-14, 0] * units.FLAM, rtol=1e-6)
 
 
 @pytest.mark.skipif('not HAS_SCIPY')
@@ -233,57 +233,54 @@ class TestEmpiricalBandpassFromFile:
     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
+        assert_quantity_allclose(f, 272.01081629459344 * 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
+        assert_quantity_allclose(f, 1.2062975715374322 * u.AA, rtol=2.5e-6)
 
     def test_avgwave(self):
         """Compare AVGWAVE with old SYNPHOT result."""
         w = self.bp.avgwave()
-        np.testing.assert_allclose(w.value, 5367.9, rtol=1e-5)
+        assert_quantity_allclose(w, 5367.9 * u.AA, rtol=1e-5)
 
     def test_barlam(self):
         """Test BARLAM (no old SYNPHOT result available)."""
         w = self.bp.barlam()
-        np.testing.assert_allclose(w.value, 5331.8945, rtol=1e-5)
+        assert_quantity_allclose(w, 5331.8945 * u.AA, rtol=1e-5)
 
     def test_pivot(self):
         """Compare PIVWV with ASTROLIB PYSYNPHOT result."""
         w = self.bp.pivot()
-        np.testing.assert_allclose(w.value, 5355.863596422962, rtol=1e-6)
+        assert_quantity_allclose(w, 5355.863596422962 * u.AA, rtol=1e-6)
 
     def test_uresp(self):
         """Compare URESP with old SYNPHOT result."""
         uresp = self.bp.unit_response(area=_area)
-        np.testing.assert_allclose(uresp.value, 3.00737e-19, rtol=1e-4)
-        assert uresp.unit == units.FLAM
+        assert_quantity_allclose(uresp, 3.00737e-19 * units.FLAM, rtol=1e-4)
 
     def test_rmswidth(self):
         w = self.bp.rmswidth()
-        np.testing.assert_allclose(w.value, 359.55954282883687, rtol=1e-4)
+        assert_quantity_allclose(w, 359.55954282883687 * u.AA, rtol=1e-4)
 
         w = self.bp.rmswidth(threshold=0.01 * u.dimensionless_unscaled)
-        np.testing.assert_allclose(w.value, 357.43298216917754, rtol=1e-4)
+        assert_quantity_allclose(w, 357.43298216917754 * u.AA, rtol=1e-4)
 
         # Invalid threshold must raise exception
         with pytest.raises(exceptions.SynphotError):
-            w = self.bp.rmswidth(threshold=0.01 * u.AA)
+            self.bp.rmswidth(threshold=0.01 * u.AA)
         with pytest.raises(exceptions.SynphotError):
-            w = self.bp.rmswidth(threshold=[0.01, 0.02])
+            self.bp.rmswidth(threshold=[0.01, 0.02])
         with pytest.raises(exceptions.SynphotError):
-            w = self.bp.rmswidth(threshold='foo')
+            self.bp.rmswidth(threshold='foo')
 
     def test_fwhm(self):
         """This also calls PHOTBW."""
         w = self.bp.fwhm()
-        np.testing.assert_allclose(w.value, 841.09, rtol=2.5e-5)
+        assert_quantity_allclose(w, 841.09 * u.AA, rtol=2.5e-5)
 
         w = self.bp.fwhm(threshold=0.01 * u.dimensionless_unscaled)
-        np.testing.assert_allclose(w.value, 836.2879507505378, rtol=2.5e-5)
+        assert_quantity_allclose(w, 836.2879507505378 * u.AA, rtol=2.5e-5)
 
         # Zero value
         w = self.bp.fwhm(wavelengths=[2e6, 2.1e6])
@@ -291,46 +288,43 @@ class TestEmpiricalBandpassFromFile:
 
         # Invalid threshold must raise exception
         with pytest.raises(exceptions.SynphotError):
-            w = self.bp.fwhm(threshold=0.01 * u.AA)
+            self.bp.fwhm(threshold=0.01 * u.AA)
         with pytest.raises(exceptions.SynphotError):
-            w = self.bp.fwhm(threshold=[0.01, 0.02])
+            self.bp.fwhm(threshold=[0.01, 0.02])
         with pytest.raises(exceptions.SynphotError):
-            w = self.bp.fwhm(threshold='foo')
+            self.bp.fwhm(threshold='foo')
 
     def test_tlambda(self):
         f = self.bp.tlambda()
-        np.testing.assert_allclose(f.value, 0.22808, rtol=1e-4)
+        assert_quantity_allclose(f, 0.22808, rtol=1e-4)
 
     def test_tpeak(self):
         """Compare TPEAK with old SYNPHOT result."""
         f = self.bp.tpeak()
-        np.testing.assert_allclose(f.value, 0.241445)
+        assert_quantity_allclose(f, 0.241445)
 
     def test_wpeak(self):
         w = self.bp.wpeak()
-        np.testing.assert_allclose(w.value, 5059.8, rtol=1e-5)
+        assert_quantity_allclose(w, 5059.8 * u.AA, rtol=1e-5)
 
     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
+        assert_quantity_allclose(w, 272.01081629459344 * u.AA, rtol=1e-6)
 
     def test_rectw(self):
         """Compare RECTW with old SYNPHOT result."""
         w = self.bp.rectwidth()
-        np.testing.assert_allclose(w.value, 1126.588, rtol=1e-5)
+        assert_quantity_allclose(w, 1126.588 * u.AA, rtol=1e-5)
 
     def test_qtlam(self):
         qtlam = self.bp.efficiency()
-        np.testing.assert_allclose(qtlam.value, 0.050901, rtol=1e-4)
-        assert qtlam.unit == u.dimensionless_unscaled
+        assert_quantity_allclose(qtlam, 0.050901, rtol=1e-4)
 
     def test_emflx(self):
         """Compare EMFLX with old SYNPHOT result."""
         f = self.bp.emflx(area=_area)
-        np.testing.assert_allclose(f.value, 3.586622e-16, rtol=2.5e-5)
-        assert f.unit == units.FLAM
+        assert_quantity_allclose(f, 3.586622e-16 * units.FLAM, rtol=2.5e-5)
 
 
 class TestBoxBandpass:
@@ -352,13 +346,23 @@ class TestBoxBandpass:
     def test_fwhm(self):
         # You would think FWHM of a box is the width but
         # not according to IRAF SYNPHOT.
-        np.testing.assert_allclose(self.bp.fwhm().value, 67.977,
-                                   rtol=1e-3)  # 0.1%
+        assert_quantity_allclose(self.bp.fwhm(), 67.977 * u.AA,
+                                 rtol=1e-3)  # 0.1%
 
     def test_taper(self):
         bp2 = self.bp.taper(np.arange(499, 501.01, 0.01) * u.nm)
         y = bp2([498.9, 499, 500, 501, 501.1] * u.nm)
-        np.testing.assert_allclose(y.value, [0, 1, 1, 1, 0])
+        assert_quantity_allclose(y, [0, 1, 1, 1, 0])
+
+    def test_integrate(self):
+        ans = 100 * u.AA
+        assert_quantity_allclose(self.bp.integrate(), ans)
+        assert_quantity_allclose(
+            self.bp.integrate(integration_type='analytical'), ans)
+
+        with pytest.raises(exceptions.SynphotError,
+                           match='flux_unit cannot be used'):
+            self.bp.integrate(integration_type='analytical', flux_unit='flam')
 
     def test_multi_n_models(self):
         """This is not allowed."""
@@ -376,48 +380,80 @@ class TestBlackBodySource:
     def test_eval(self):
         w = np.arange(3000, 3100, 10)
         y = self.sp(w)
-        np.testing.assert_allclose(
-            y.value,
+        assert_quantity_allclose(
+            y,
             [0.00019318, 0.00019623, 0.0001993, 0.00020238, 0.00020549,
-             0.00020861, 0.00021175, 0.00021491, 0.00021809, 0.00022128],
+             0.00020861, 0.00021175, 0.00021491, 0.00021809,
+             0.00022128] * units.PHOTLAM,
             rtol=2.5e-3)
 
+    def test_integrate(self):
+        ans_photlam = 12.39167258 * (u.ph / (u.cm * u.cm * u.s))
+        ans_flam = 2.62716011e-11 * (u.erg / (u.cm * u.cm * u.s))
+
+        assert_quantity_allclose(self.sp.integrate(), ans_photlam, rtol=1e-5)
+        assert_quantity_allclose(
+            self.sp.integrate(flux_unit='flam'), ans_flam, rtol=1e-5)
+        assert_quantity_allclose(
+            self.sp.integrate(integration_type='analytical', flux_unit='flam'),
+            ans_flam, rtol=5e-3)
+
+    @pytest.mark.xfail(reason='Cannot convert unit in analytical mode')
+    def test_integrate_fixme(self):
+        """Merge this into ``test_integrate()`` above when fixed."""
+        ans_photlam = 12.39167258 * (u.ph / (u.cm * u.cm * u.s))
+        assert_quantity_allclose(
+            self.sp.integrate(integration_type='analytical'), ans_photlam)
+
 
 class TestGaussianSource:
     """Test source spectrum with GaussianFlux1D model."""
     def setup_class(self):
-        tf = 4.96611456e-12 * u.erg / (u.cm * u.cm * u.s)
+        tf = 4.96611456e-12 * (u.erg / (u.cm * u.cm * u.s))
         self.sp = SourceSpectrum(
             GaussianFlux1D, total_flux=tf, mean=4000, fwhm=100)
 
     def test_eval(self):
         y = self.sp([3900, 4000, 4060])
-        np.testing.assert_allclose(
-            y.value, [0.00058715, 0.00939437, 0.00346246], rtol=1e-5)
+        assert_quantity_allclose(
+            y, [0.00058715, 0.00939437, 0.00346246] * units.PHOTLAM, rtol=1e-5)
 
     def test_totalflux(self):
+        """Test Gaussian source integration.
+
+        .. note::
+
+            Analytic integral is more accurate because it does not rely on
+            waveset definition.
+
+        """
         # PHOTLAM
-        f = self.sp.integrate()
-        np.testing.assert_allclose(f.value, 1, rtol=1e-5)
-        assert f.unit == u.photon / (u.cm**2 * u.s)
+        f_ans = 1 * (u.ph / (u.cm**2 * u.s))
+        assert_quantity_allclose(self.sp.integrate(), f_ans, rtol=1e-5)
+        assert_quantity_allclose(
+            self.sp.integrate(integration_type='analytical'), f_ans)
 
         # FLAM
-        x0 = (400 * u.nm).to(u.AA).value
-        fwhm = (10 * u.nm).to(u.AA)
+        x0 = 400 * u.nm
+        fwhm = 10 * u.nm
         sp2 = SourceSpectrum(
             GaussianFlux1D, total_flux=1, mean=x0, fwhm=fwhm)
-        val = sp2.integrate(flux_unit=units.FLAM).value
-        np.testing.assert_allclose(val, 1, rtol=1e-3)
+        val_ans = 1 * (u.erg / (u.cm * u.cm * u.s))
+        assert_quantity_allclose(
+            sp2.integrate(flux_unit=units.FLAM), val_ans, rtol=1e-3)
+        assert_quantity_allclose(
+            sp2.integrate(flux_unit=units.FLAM, integration_type='analytical'),
+            val_ans)
 
     def test_symmetry(self):
-        np.testing.assert_allclose(self.sp(3950), self.sp(4050))
+        assert_quantity_allclose(self.sp(3950), self.sp(4050))
 
     def test_fwhm(self):
         """Should round-trip back to the same bandpass FWHM."""
         m = self.sp.model
         bp = SpectralElement(
             Gaussian1D, mean=m.mean, amplitude=m.amplitude, stddev=m.stddev)
-        np.testing.assert_allclose(bp.fwhm().value, 100, rtol=1e-3)  # 0.1%
+        assert_quantity_allclose(bp.fwhm(), 100 * u.AA, rtol=1e-3)  # 0.1%
 
     def test_alt_source(self):
         """Same source, different way to init."""
@@ -425,43 +461,68 @@ class TestGaussianSource:
             GaussianFlux1D, amplitude=self.sp.model.amplitude.value,
             mean=self.sp.model.mean.value, stddev=self.sp.model.stddev.value)
         w = [3900, 4000, 4060] * u.AA
-        np.testing.assert_allclose(sp2(w), self.sp(w))
+        assert_quantity_allclose(sp2(w), self.sp(w))
 
 
 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
+    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)
+    assert_quantity_allclose(tf, flux, rtol=1e-4)
 
 
 class TestPowerLawSource:
     """Test source spectrum with PowerLawFlux1D model."""
     def setup_class(self):
-        self.sp = SourceSpectrum(PowerLawFlux1D, amplitude=1, x_0=6000,
-                                 alpha=4)
+        self.sp = SourceSpectrum(PowerLawFlux1D, amplitude=1 * units.PHOTLAM,
+                                 x_0=6000 * u.AA, alpha=4)
+        self.w = np.arange(3000, 3100, 10) * u.AA
 
     def test_no_default_wave(self):
         assert self.sp.waverange == [None, None]
 
-        with pytest.raises(exceptions.SynphotError):
+        with pytest.raises(exceptions.SynphotError,
+                           match='waveset is undefined'):
             self.sp(None)
 
     def test_eval(self):
-        w = np.arange(3000, 3100, 10)
-        y = self.sp(w)
-        np.testing.assert_allclose(
-            y.value,
+        y = self.sp(self.w)
+        assert_quantity_allclose(
+            y,
             [16, 15.78843266, 15.58035072, 15.37568551, 15.17436992,
-             14.97633838, 14.78152682, 14.5898726, 14.40131453, 14.21579277],
+             14.97633838, 14.78152682, 14.5898726, 14.40131453,
+             14.21579277] * units.PHOTLAM,
             rtol=1e-6)
 
     def test_normalization(self):
-        np.testing.assert_allclose(self.sp(600 * u.nm).value, 1)
+        assert_quantity_allclose(self.sp(600 * u.nm), 1 * units.PHOTLAM)
+
+    def test_integrate(self):
+        ans_photlam = 1357.75787527 * (u.ph / (u.cm * u.cm * u.s))
+        ans_flam = 8.8608168e-09 * (u.erg / (u.cm * u.cm * u.s))
+        assert_quantity_allclose(
+            self.sp.integrate(wavelengths=self.w), ans_photlam)
+        assert_quantity_allclose(
+            self.sp.integrate(wavelengths=self.w, flux_unit='flam'), ans_flam)
+        assert_quantity_allclose(
+            self.sp.integrate(wavelengths=self.w,
+                              integration_type='analytical'),
+            ans_photlam, rtol=1e-4)
+
+    @pytest.mark.xfail(reason='Cannot convert unit of analytic integral')
+    def test_integrate_wontfix(self):
+        """Powerlaw in one flux unit might not be powerlaw anymore in
+        another, so we cannot convert flux unit of analytical integration
+        easily.
+        """
+        ans_flam = 8.8608168e-09 * (u.erg / (u.cm * u.cm * u.s))
+        assert_quantity_allclose(
+            self.sp.integrate(wavelengths=self.w, flux_unit='flam',
+                              integration_type='analytical'), ans_flam)
 
 
 class TestBuildModels:
@@ -471,56 +532,56 @@ class TestBuildModels:
             BrokenPowerLaw1D, amplitude=1, x_break=6000, alpha_1=1,
             alpha_2=4)
         y = sp([5000, 6000, 7000])
-        np.testing.assert_allclose(y.value, [1.2, 1, 0.53977509])
+        assert_quantity_allclose(y, [1.2, 1, 0.53977509] * units.PHOTLAM)
 
     def test_Const1D(self):
         sp = SourceSpectrum(Const1D, amplitude=1)
         y = sp([1, 1000, 1e6])
-        np.testing.assert_array_equal(y.value, 1)
+        assert_quantity_allclose(y, 1 * units.PHOTLAM, rtol=0)
 
     def test_ConstFlux1D(self):
         sp = SourceSpectrum(ConstFlux1D, amplitude=1 * u.Jy)
-        w = [1, 1000, 1e6]
-        y = units.convert_flux(w, sp(w), u.Jy)
-        np.testing.assert_allclose(y.value, 1)
+        w = [1, 1000, 1e6] * u.AA
+        with u.add_enabled_equivalencies(u.spectral_density(w)):
+            assert_quantity_allclose(sp(w), 1 * u.Jy)
 
     def test_ExponentialCutoffPowerLaw1D(self):
         sp = SourceSpectrum(
             ExponentialCutoffPowerLaw1D, amplitude=1, x_0=6000,
             x_cutoff=10000, alpha=4)
         y = sp([5000, 6000, 10000])
-        np.testing.assert_allclose(
-            y.value, [1.25770198, 0.54881164, 0.04767718])
+        assert_quantity_allclose(
+            y, [1.25770198, 0.54881164, 0.04767718] * units.PHOTLAM)
 
     def test_GaussianAbsorption1D(self):
         """This should be unitless, not a source spectrum."""
         bp = SpectralElement(
             GaussianAbsorption1D, amplitude=0.8, mean=5500, stddev=50)
         y = bp([5300, 5500, 5700])
-        np.testing.assert_allclose(y.value, [0.99973163, 0.2, 0.99973163])
+        assert_quantity_allclose(y, [0.99973163, 0.2, 0.99973163])
 
     def test_LogParabola1D(self):
         sp = SourceSpectrum(
             LogParabola1D, amplitude=1, x_0=6000, alpha=1, beta=4)
         y = sp([5000, 6000, 7000])
-        np.testing.assert_allclose(y.value, [1.0505953, 1, 0.77942375])
+        assert_quantity_allclose(y, [1.0505953, 1, 0.77942375] * units.PHOTLAM)
 
     def test_Lorentz1D(self):
         sp = SourceSpectrum(Lorentz1D, amplitude=1, x_0=6000, fwhm=100)
         y = sp([5000, 6000, 7000])
-        np.testing.assert_allclose(
-            y.value, [0.00249377, 1, 0.00249377], rtol=1e-5)
+        assert_quantity_allclose(
+            y, [0.00249377, 1, 0.00249377] * units.PHOTLAM, rtol=1e-5)
 
     def test_RickerWavelet1D(self):
         sp = SourceSpectrum(RickerWavelet1D, amplitude=1, x_0=6000, sigma=100)
         y = sp([5000, 6000, 7000])
-        np.testing.assert_allclose(
-            y.value, [-1.90946235e-20, 1, -1.90946235e-20])
+        assert_quantity_allclose(
+            y, [-1.90946235e-20, 1, -1.90946235e-20] * units.PHOTLAM)
 
     def test_PowerLaw1D(self):
         sp = SourceSpectrum(PowerLaw1D, amplitude=1, x_0=6000, alpha=4)
         y = sp([5000, 6000, 7000])
-        np.testing.assert_allclose(y.value, [2.0736, 1, 0.53977509])
+        assert_quantity_allclose(y, [2.0736, 1, 0.53977509] * units.PHOTLAM)
 
 
 @pytest.mark.skipif('not HAS_SCIPY')
@@ -595,21 +656,21 @@ class TestForceExtrap:
                             lookup_table=[0.5, 0.6, 10.6, 1.5], fill_value=0)
         sp.z = z
         w = [900, 4300]
-        np.testing.assert_allclose(sp(w).value, 0)  # No extrapolation
+        assert_quantity_allclose(sp(w), 0 * units.PHOTLAM)  # No extrapolation
 
         is_forced = sp.force_extrapolation()  # Force extrapolation
         assert is_forced
-        np.testing.assert_allclose(sp(w).value, [0.5, 1.5])
+        assert_quantity_allclose(sp(w), [0.5, 1.5] * units.PHOTLAM)
 
     def test_analytical(self):
         """Forcing is not possible."""
         sp = SourceSpectrum(GaussianFlux1D, mean=5500, total_flux=1, fwhm=10)
         w = [100, 10000]
-        np.testing.assert_allclose(sp(w).value, 0)
+        assert_quantity_allclose(sp(w), 0 * units.PHOTLAM)
 
         is_forced = sp.force_extrapolation()
         assert not is_forced
-        np.testing.assert_allclose(sp(w).value, 0)
+        assert_quantity_allclose(sp(w), 0 * units.PHOTLAM)
 
 
 @pytest.mark.skipif('not HAS_SCIPY')
@@ -653,7 +714,8 @@ class TestNormalize:
         ct_rate = obs.countrate(_area)
 
         # 0.2% agreement with IRAF SYNPHOT COUNTRATE
-        np.testing.assert_allclose(ct_rate.value, ans_countrate, rtol=2e-3)
+        assert_quantity_allclose(
+            ct_rate, ans_countrate * (u.ct / u.s), rtol=2e-3)
 
     @pytest.mark.parametrize(
         ('sp_type', 'rn_val', 'ans_countrate'),
@@ -791,9 +853,10 @@ class TestWaveset:
         g3 = SourceSpectrum(
             GaussianFlux1D, total_flux=totflux, mean=7500, fwhm=5)
         sp = SpectralElement(Box1D, x_0=1000, width=1) * (g1 + g2 + g3)
-        np.testing.assert_allclose(
-            sp.waveset.value[::100],
-            [999.49, 1000.49, 5019.95906231, 6699.59062307, 7509.7672007])
+        assert_quantity_allclose(
+            sp.waveset[::100],
+            [999.49, 1000.49, 5019.95906231, 6699.59062307,
+             7509.7672007] * u.AA)
 
     def test_redshift(self):
         tf_unit = u.erg / (u.cm * u.cm * u.s)
@@ -801,8 +864,8 @@ class TestWaveset:
             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]
-        np.testing.assert_allclose(sp.waveset.value[::25], m(w_step25_z0))
+        w_step25_z0 = [4978.76695499, 4989.3834775, 5000, 5010.6165225] * u.AA
+        assert_quantity_allclose(sp.waveset[::25], m(w_step25_z0))
 
     def test_redshift_none(self):
         sp = SourceSpectrum(Const1D, amplitude=1, z=1.3)
@@ -873,21 +936,21 @@ class TestRedShift:
     def test_composite_redshift(self):
         sp2 = self.sp_z0 + self.sp  # centers: 5000, 11500
         sp2.z = 0.5  # centers: 7500, 17250
-        np.testing.assert_allclose(sp2([7500, 17250]), self.sp_z0(5000))
+        assert_quantity_allclose(sp2([7500, 17250]), self.sp_z0(5000))
 
     def test_const_flux_redshift(self):
         """Constant flux in Jy is not constant in PHOTLAM."""
         sp_z0 = SourceSpectrum(ConstFlux1D, amplitude=1 * u.Jy)
         sp = SourceSpectrum(ConstFlux1D, amplitude=1 * u.Jy, z=1.3)
-        np.testing.assert_allclose(sp_z0(3000), sp(6900))
+        assert_quantity_allclose(sp_z0(3000), sp(6900))
 
     def test_conserve_flux_redshift(self):
         """Test redshift behavior that conserves flux."""
         sp = SourceSpectrum(self.sp_z0.model, z=1.3, z_type='conserve_flux')
         fac = 1 / (1 + sp.z)
         wave = [5000, 11500]
-        np.testing.assert_allclose(sp(wave), self.sp(wave) * fac)
-        np.testing.assert_allclose(sp.integrate(), self.sp_z0.integrate())
+        assert_quantity_allclose(sp(wave), self.sp(wave) * fac)
+        assert_quantity_allclose(sp.integrate(), self.sp_z0.integrate())
 
 
 @pytest.mark.skipif('not HAS_SCIPY')
@@ -907,24 +970,25 @@ class TestMathOperators:
     def test_source_add(self):
         """Compare with ASTROLIB PYSYNPHOT."""
         ans = self.sp_1 + self.sp_2
-        np.testing.assert_allclose(
-            ans(ans.waveset).value,
+        assert_quantity_allclose(
+            ans(ans.waveset),
             [0.00976521, 0.01681283, 0.01970276, 0.01998463, 0.0197387,
-             0.01985257, 0.02337638, 0.00978454], rtol=1e-4)
+             0.01985257, 0.02337638, 0.00978454] * units.PHOTLAM, rtol=1e-4)
 
     def test_source_sub(self):
         """Compare with ASTROLIB PYSYNPHOT."""
         ans = self.sp_1 - self.sp_2
-        np.testing.assert_allclose(
-            ans(ans.waveset).value,
+        assert_quantity_allclose(
+            ans(ans.waveset),
             [-9.76520783e-03, -2.71758275e-03, 1.72346256e-04, -9.29051118e-05,
-             1.69629843e-04, 2.83499328e-04, 3.80731187e-03, -9.78453651e-03],
+             1.69629843e-04, 2.83499328e-04, 3.80731187e-03,
+             -9.78453651e-03] * units.PHOTLAM,
             rtol=1e-4)
 
     def test_source_addsub_circular(self):
         """sp = sp + sp - sp"""
         ans = self.sp_1 + self.sp_1 - self.sp_1
-        np.testing.assert_array_equal(ans(ans.waveset), self.sp_1(ans.waveset))
+        assert_quantity_allclose(ans(ans.waveset), self.sp_1(ans.waveset))
 
     def test_source_addsub_exception(self):
         with pytest.raises(exceptions.IncompatibleSources):
@@ -932,41 +996,44 @@ class TestMathOperators:
 
     @pytest.mark.parametrize('x', [2, 2 * u.dimensionless_unscaled])
     def test_source_mul_div_scalar(self, x):
-        w = self.sp_1.waveset.value
+        w = self.sp_1.waveset
 
         ans1 = self.sp_1 * x
-        np.testing.assert_allclose(
-            ans1(w).value, [0, 0.01409552, 0.02013646, 0.02718424, 0],
+        assert_quantity_allclose(
+            ans1(w),
+            [0, 0.01409552, 0.02013646, 0.02718424, 0] * units.PHOTLAM,
             rtol=1e-6)
 
         # rmul does not work with Quantity
         if not isinstance(x, u.Quantity):
             ans2 = x * self.sp_1
-            np.testing.assert_array_equal(ans1(w), ans2(w))
+            assert_quantity_allclose(ans1(w), ans2(w), rtol=0)
 
         ans3 = self.sp_1 / x
-        np.testing.assert_allclose(
-            ans3(w).value, [0, 0.00352388, 0.00503411, 0.00679606, 0],
-            atol=1e-7)
+        assert_quantity_allclose(
+            ans3(w),
+            [0, 0.00352388, 0.00503411, 0.00679606, 0] * units.PHOTLAM,
+            atol=1e-7 * units.PHOTLAM)
 
     def test_source_mul_div_spec(self):
         """Compare mul with ASTROLIB PYSYNPHOT. Also test bp * sp."""
         ans1 = self.sp_1 * self.bp_1
         ans2 = self.bp_1 * self.sp_1
         w = ans1.waveset[:-1]
-        np.testing.assert_allclose(
-            ans1(w).value,
+        assert_quantity_allclose(
+            ans1(w),
             [0, 3.52381254e-04, 7.04792712e-04, 2.01360717e-03, 3.97184014e-03,
-             4.03718269e-05, 0], rtol=1e-4)
-        np.testing.assert_array_equal(ans1(w), ans2(w))
+             4.03718269e-05, 0] * units.PHOTLAM, rtol=1e-4)
+        assert_quantity_allclose(ans1(w), ans2(w), rtol=0)
 
         ans3 = self.sp_1 / self.bp_1
-        np.testing.assert_allclose(
-            ans3(w).value,
-            [0, 0.14095528, 0.07048066, 0.05034117, 0.04413243, 4.57601236, 0])
+        assert_quantity_allclose(
+            ans3(w),
+            [0, 0.14095528, 0.07048066, 0.05034117, 0.04413243, 4.57601236,
+             0] * units.PHOTLAM)
 
         ans4 = self.sp_1 / self.sp_1
-        np.testing.assert_allclose(
+        assert_quantity_allclose(
             ans4([4000, 5000, 6000]), 1 * u.dimensionless_unscaled)
 
         # Dividing throughput by flux does not make sense.
@@ -1003,27 +1070,27 @@ class TestMathOperators:
 
     @pytest.mark.parametrize('x', [2.0, 2.0 * u.dimensionless_unscaled])
     def test_bandpass_mul_div_scalar(self, x):
-        w = self.bp_1.waveset.value
+        w = self.bp_1.waveset
 
         ans1 = self.bp_1 * x
-        np.testing.assert_allclose(ans1(w).value, [0, 0.2, 0.4, 0.6, 0])
+        assert_quantity_allclose(ans1(w), [0, 0.2, 0.4, 0.6, 0])
 
         # rmul does not work with Quantity
         if not isinstance(x, u.Quantity):
             ans2 = x * self.bp_1
-            np.testing.assert_array_equal(ans1(w), ans2(w))
+            assert_quantity_allclose(ans1(w), ans2(w), rtol=0)
 
         ans3 = self.bp_1 / x
-        np.testing.assert_allclose(ans3(w).value, [0, 0.05, 0.1, 0.15, 0])
+        assert_quantity_allclose(ans3(w), [0, 0.05, 0.1, 0.15, 0])
 
     def test_bandpass_mul_div_bandpass(self):
         ans1 = self.bp_1 * self.bp_1
-        np.testing.assert_allclose(
-            ans1(ans1.waveset).value, [0, 0.01, 0.04, 0.09, 0])
+        assert_quantity_allclose(
+            ans1(ans1.waveset), [0, 0.01, 0.04, 0.09, 0])
 
         w = [4000.1, 5000, 5900]  # Avoid div by zero
         ans2 = self.bp_1 / self.bp_1
-        np.testing.assert_allclose(ans2(w), 1)
+        assert_quantity_allclose(ans2(w), 1)
 
     def test_bandpass_mul_div_exceptions(self):
         """Only mul is tested but truediv uses the same validation."""
diff --git a/synphot/tests/test_units.py b/synphot/tests/test_units.py
index f4844cb..8132847 100644
--- a/synphot/tests/test_units.py
+++ b/synphot/tests/test_units.py
@@ -3,6 +3,8 @@
 
 .. note:: VEGAMAG conversion is tested in test_spectrum.py.
 
+.. note:: spectral_density_integrated is tested in astropy>=4.1.
+
 """
 
 # THIRD-PARTY
diff --git a/synphot/tests/test_utils.py b/synphot/tests/test_utils.py
index b83b133..57311e8 100644
--- a/synphot/tests/test_utils.py
+++ b/synphot/tests/test_utils.py
@@ -12,7 +12,7 @@ import pytest
 from astropy import units as u
 
 # LOCAL
-from .. import exceptions, utils
+from .. import exceptions, utils, units
 
 
 @pytest.mark.parametrize(
@@ -30,18 +30,19 @@ def test_overlap_status(a, b, ans):
 
 def test_validate_totalflux():
     """Test integrated flux validation."""
-    # Valid integrated flux
     utils.validate_totalflux(0.01)
+    utils.validate_totalflux(0.01 * units.FLAM)
 
-    # Invalid integrated flux
-    with pytest.raises(exceptions.SynphotError):
-        utils.validate_totalflux(-0.01)
-    with pytest.raises(exceptions.SynphotError):
-        utils.validate_totalflux(0)
-    with pytest.raises(exceptions.SynphotError):
-        utils.validate_totalflux(np.inf)
+
+@pytest.mark.parametrize(
+    'val', (-0.01, -0.01 * units.FLAM,
+            0, 0 * units.PHOTLAM,
+            np.inf, np.inf * u.Jy,
+            np.nan, np.nan * units.FLAM))
+def test_validate_totalflux_invalid(val):
+    """Invalid integrated flux."""
     with pytest.raises(exceptions.SynphotError):
-        utils.validate_totalflux(np.nan)
+        utils.validate_totalflux(val)
 
 
 def test_validate_wavelengths():