Repository: spacetelescope/synphot_refactor

Issue:
TST: Move bbfunc test from astropy to here
https://github.com/astropy/astropy/blob/09b7796bbd4eadc793bb48fc5c8554bb5bd4af0d/astropy/modeling/tests/test_blackbody.py#L97-L111

```python
def test_blackbody_synphot():
    """Test that it is consistent with IRAF SYNPHOT BBFUNC."""
    # Solid angle of solar radius at 1 kpc
    fac = np.pi * (const.R_sun / const.kpc) ** 2 * u.sr

    with np.errstate(all='ignore'):
        flux = blackbody_nu([100, 1, 1000, 1e4, 1e5] * u.AA, 5000) * fac
    assert flux.unit == FNU

    # Special check for overflow value (SYNPHOT gives 0)
    assert np.log10(flux[0].value) < -143

    np.testing.assert_allclose(
        flux.value[1:], [0, 2.01950807e-34, 3.78584515e-26, 1.90431881e-27],
        rtol=0.01)  # 1% accuracy
```

This test will not run for the refactored `BlackBody` code after astropy/astropy#9282 is merged.

Discussion:

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index eee74f6..218ccd8 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -6,6 +6,8 @@
   object, respectively. [#243]
 - Dropped support for Python 3.5 and ``astropy`` 2.x. This version is only
   compatible with Python 3.6 or later and ``astropy`` 3.x or later. [#243]
+- Added support for ``RickerWavelet1D`` model that is the renamed version
+  of ``MexicanHat1D`` model to be consistent with ``astropy`` 4.0. [#250]
 
 0.2.1 (2019-12-20)
 ==================
diff --git a/docs/index.rst b/docs/index.rst
index 0356759..9427f58 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -315,6 +315,9 @@ Also imports this C-extension to local namespace:
 
    synphot/c_ext
 
+.. automodapi:: synphot.blackbody
+   :no-inheritance-diagram:
+
 .. automodapi:: synphot.config
    :no-inheritance-diagram:
 
diff --git a/docs/synphot/overview.rst b/docs/synphot/overview.rst
index a555964..4bca019 100644
--- a/docs/synphot/overview.rst
+++ b/docs/synphot/overview.rst
@@ -109,7 +109,7 @@ and special notes:
 |`~synphot.models.Lorentz1D`                              |**synphot**|Like `astropy.modeling.functional_models.Lorentz1D`           |
 |                                                         |           |but with ``sampleset`` and ``bounding_box``.                  |
 +---------------------------------------------------------+-----------+--------------------------------------------------------------+
-|`~synphot.models.MexicanHat1D`                           |**synphot**|Like `astropy.modeling.functional_models.MexicanHat1D`        |
+|`~synphot.models.RickerWavelet1D`                        |**synphot**|Like `astropy.modeling.functional_models.RickerWavelet1D`     |
 |                                                         |           |but with ``sampleset`` and ``bounding_box``.                  |
 +---------------------------------------------------------+-----------+--------------------------------------------------------------+
 |`~astropy.modeling.powerlaws.PowerLaw1D`                 |Astropy    ||note_flux_conv_incorrect|                                    |
diff --git a/setup.cfg b/setup.cfg
index 742e077..8e36b72 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -10,8 +10,6 @@ xfail_strict = true
 filterwarnings =
     ignore:can't resolve package from __spec__
     ignore:numpy.ufunc size changed:RuntimeWarning
-    ignore:BlackBody provides the same capabilities
-    ignore:The MexicanHat1D class is deprecated
 
 [metadata]
 name = synphot
diff --git a/synphot/blackbody.py b/synphot/blackbody.py
new file mode 100644
index 0000000..13c6b7d
--- /dev/null
+++ b/synphot/blackbody.py
@@ -0,0 +1,260 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+"""Model and functions related to blackbody radiation.
+
+.. note::
+
+    This was ``astropy.modeling.blackbody`` module that was
+    deprecated in ``astropy`` 4.0. The content is copied here
+    so we can still use it without deprecation warning.
+    Eventually, when we can fully pass in unit handling directly
+    to Astropy models, we can remove this module. See
+    https://github.com/astropy/astropy/pull/9282 and
+    https://github.com/spacetelescope/synphot_refactor/pull/224 .
+
+"""
+import warnings
+
+import numpy as np
+
+from astropy import constants as const
+from astropy import units as u
+from astropy.modeling.core import Fittable1DModel
+from astropy.modeling.parameters import Parameter
+from astropy.utils.exceptions import AstropyUserWarning
+
+from .units import FNU, FLAM
+
+__all__ = ['BlackBody1D', 'blackbody_nu', 'blackbody_lambda']
+
+with warnings.catch_warnings():
+    warnings.simplefilter('ignore', RuntimeWarning)
+    _has_buggy_expm1 = np.isnan(np.expm1(1000)) or np.isnan(np.expm1(1e10))
+
+
+class BlackBody1D(Fittable1DModel):
+    """
+    One dimensional blackbody model.
+
+    Parameters
+    ----------
+    temperature : :class:`~astropy.units.Quantity`
+        Blackbody temperature.
+    bolometric_flux : :class:`~astropy.units.Quantity`
+        The bolometric flux of the blackbody (i.e., the integral over the
+        spectral axis).
+
+    Notes
+    -----
+
+    Model formula:
+
+        .. math:: f(x) = \\pi B_{\\nu} f_{\\text{bolometric}} / (\\sigma  T^{4})
+
+    Examples
+    --------
+    >>> from astropy import units as u
+    >>> from synphot.blackbody import BlackBody1D
+    >>> bb = BlackBody1D()
+    >>> bb(6000 * u.AA)  # doctest: +FLOAT_CMP
+    <Quantity 1.3585381201978953e-15 erg / (cm2 Hz s)>
+
+    .. plot::
+        :include-source:
+
+        import numpy as np
+        import matplotlib.pyplot as plt
+        from astropy import units as u
+        from astropy.visualization import quantity_support
+        from synphot.blackbody import BlackBody1D
+        from synphot.units import FLAM
+
+        bb = BlackBody1D(temperature=5778*u.K)
+        wav = np.arange(1000, 110000) * u.AA
+        flux = bb(wav).to(FLAM, u.spectral_density(wav))
+
+        with quantity_support():
+            plt.figure()
+            plt.semilogx(wav, flux)
+            plt.axvline(bb.lambda_max.to(u.AA).value, ls='--')
+            plt.show()
+
+    """  # noqa
+
+    # We parametrize this model with a temperature and a bolometric flux. The
+    # bolometric flux is the integral of the model over the spectral axis. This
+    # is more useful than simply having an amplitude parameter.
+    temperature = Parameter(default=5000, min=0, unit=u.K)
+    bolometric_flux = Parameter(default=1, min=0, unit=u.erg / u.cm ** 2 / u.s)
+
+    # We allow values without units to be passed when evaluating the model, and
+    # in this case the input x values are assumed to be frequencies in Hz.
+    _input_units_allow_dimensionless = True
+
+    # We enable the spectral equivalency by default for the spectral axis
+    input_units_equivalencies = {'x': u.spectral()}
+
+    def evaluate(self, x, temperature, bolometric_flux):
+        """Evaluate the model.
+
+        Parameters
+        ----------
+        x : float, `~numpy.ndarray`, or `~astropy.units.Quantity`
+            Frequency at which to compute the blackbody. If no units are given,
+            this defaults to Hz.
+
+        temperature : float, `~numpy.ndarray`, or `~astropy.units.Quantity`
+            Temperature of the blackbody. If no units are given, this defaults
+            to Kelvin.
+
+        bolometric_flux : float, `~numpy.ndarray`, or `~astropy.units.Quantity`
+            Desired integral for the blackbody.
+
+        Returns
+        -------
+        y : number or ndarray
+            Blackbody spectrum. The units are determined from the units of
+            ``bolometric_flux``.
+        """
+
+        # We need to make sure that we attach units to the temperature if it
+        # doesn't have any units. We do this because even though blackbody_nu
+        # can take temperature values without units, the / temperature ** 4
+        # factor needs units to be defined.
+        if isinstance(temperature, u.Quantity):
+            temperature = temperature.to(u.K, equivalencies=u.temperature())
+        else:
+            temperature = u.Quantity(temperature, u.K)
+
+        # We normalize the returned blackbody so that the integral would be
+        # unity, and we then multiply by the bolometric flux. A normalized
+        # blackbody has f_nu = pi * B_nu / (sigma * T^4), which is what we
+        # calculate here. We convert to 1/Hz to make sure the units are
+        # simplified as much as possible, then we multiply by the bolometric
+        # flux to get the normalization right.
+        fnu = ((np.pi * u.sr * blackbody_nu(x, temperature) /
+                const.sigma_sb / temperature ** 4).to(1 / u.Hz) *
+               bolometric_flux)
+
+        # If the bolometric_flux parameter has no unit, we should drop the /Hz
+        # and return a unitless value. This occurs for instance during fitting,
+        # since we drop the units temporarily.
+        if hasattr(bolometric_flux, 'unit'):
+            return fnu
+        else:
+            return fnu.value
+
+    @property
+    def input_units(self):
+        # The input units are those of the 'x' value, which should always be
+        # Hz. Because we do this, and because input_units_allow_dimensionless
+        # is set to True, dimensionless values are assumed to be in Hz.
+        return {'x': u.Hz}
+
+    def _parameter_units_for_data_units(self, inputs_unit, outputs_unit):
+        return {'temperature': u.K,
+                'bolometric_flux': outputs_unit['y'] * u.Hz}
+
+    @property
+    def lambda_max(self):
+        """Peak wavelength when the curve is expressed as power density."""
+        return const.b_wien / self.temperature
+
+
+def blackbody_nu(in_x, temperature):
+    """Calculate blackbody flux per steradian, :math:`B_{\\nu}(T)`.
+
+    .. note::
+
+        Use `numpy.errstate` to suppress Numpy warnings, if desired.
+
+    .. warning::
+
+        Output values might contain ``nan`` and ``inf``.
+
+    Parameters
+    ----------
+    in_x : number, array_like, or `~astropy.units.Quantity`
+        Frequency, wavelength, or wave number.
+        If not a Quantity, it is assumed to be in Hz.
+
+    temperature : number, array_like, or `~astropy.units.Quantity`
+        Blackbody temperature.
+        If not a Quantity, it is assumed to be in Kelvin.
+
+    Returns
+    -------
+    flux : `~astropy.units.Quantity`
+        Blackbody monochromatic flux in
+        :math:`erg \\; cm^{-2} s^{-1} Hz^{-1} sr^{-1}`.
+
+    Raises
+    ------
+    ValueError
+        Invalid temperature.
+
+    ZeroDivisionError
+        Wavelength is zero (when converting to frequency).
+
+    """
+    # Convert to units for calculations, also force double precision
+    with u.add_enabled_equivalencies(u.spectral() + u.temperature()):
+        freq = u.Quantity(in_x, u.Hz, dtype=np.float64)
+        temp = u.Quantity(temperature, u.K, dtype=np.float64)
+
+    # Check if input values are physically possible
+    if np.any(temp < 0):
+        raise ValueError(f'Temperature should be positive: {temp}')
+    if not np.all(np.isfinite(freq)) or np.any(freq <= 0):
+        warnings.warn('Input contains invalid wavelength/frequency value(s)',
+                      AstropyUserWarning)
+
+    log_boltz = const.h * freq / (const.k_B * temp)
+    boltzm1 = np.expm1(log_boltz)
+
+    if _has_buggy_expm1:
+        # Replace incorrect nan results with infs--any result of 'nan' is
+        # incorrect unless the input (in log_boltz) happened to be nan to begin
+        # with.  (As noted in #4393 ideally this would be replaced by a version
+        # of expm1 that doesn't have this bug, rather than fixing incorrect
+        # results after the fact...)
+        boltzm1_nans = np.isnan(boltzm1)
+        if np.any(boltzm1_nans):
+            if boltzm1.isscalar and not np.isnan(log_boltz):
+                boltzm1 = np.inf
+            else:
+                boltzm1[np.where(~np.isnan(log_boltz) & boltzm1_nans)] = np.inf
+
+    # Calculate blackbody flux
+    bb_nu = (2.0 * const.h * freq ** 3 / (const.c ** 2 * boltzm1))
+    flux = bb_nu.to(FNU, u.spectral_density(freq))
+
+    return flux / u.sr  # Add per steradian to output flux unit
+
+
+def blackbody_lambda(in_x, temperature):
+    """Like :func:`blackbody_nu` but for :math:`B_{\\lambda}(T)`.
+
+    Parameters
+    ----------
+    in_x : number, array_like, or `~astropy.units.Quantity`
+        Frequency, wavelength, or wave number.
+        If not a Quantity, it is assumed to be in Angstrom.
+
+    temperature : number, array_like, or `~astropy.units.Quantity`
+        Blackbody temperature.
+        If not a Quantity, it is assumed to be in Kelvin.
+
+    Returns
+    -------
+    flux : `~astropy.units.Quantity`
+        Blackbody monochromatic flux in
+        :math:`erg \\; cm^{-2} s^{-1} \\mathring{A}^{-1} sr^{-1}`.
+
+    """
+    if getattr(in_x, 'unit', None) is None:
+        in_x = u.Quantity(in_x, u.AA)
+
+    bb_nu = blackbody_nu(in_x, temperature) * u.sr  # Remove sr for conversion
+    flux = bb_nu.to(FLAM, u.spectral_density(in_x))
+
+    return flux / u.sr  # Add per steradian to output flux unit
diff --git a/synphot/compat.py b/synphot/compat.py
index 66c37bb..bcd3c94 100644
--- a/synphot/compat.py
+++ b/synphot/compat.py
@@ -11,7 +11,6 @@ except ImportError:
 else:
     HAS_SPECUTILS = True
 
-__all__ = ['ASTROPY_LT_2_0', 'ASTROPY_LT_4_0', 'HAS_SPECUTILS']
+__all__ = ['ASTROPY_LT_4_0', 'HAS_SPECUTILS']
 
-ASTROPY_LT_2_0 = not minversion(astropy, '2.0')
 ASTROPY_LT_4_0 = not minversion(astropy, '4.0')
diff --git a/synphot/models.py b/synphot/models.py
index 23249c4..0ec273a 100644
--- a/synphot/models.py
+++ b/synphot/models.py
@@ -21,19 +21,21 @@ from astropy.utils.exceptions import AstropyUserWarning
 
 # LOCAL
 from . import units
-from .compat import ASTROPY_LT_2_0, ASTROPY_LT_4_0
+from .compat import ASTROPY_LT_4_0
 from .exceptions import SynphotError
 from .utils import merge_wavelengths
 
 if ASTROPY_LT_4_0:
     from astropy.modeling.core import _CompoundModel as CompoundModel
+    from astropy.modeling.models import MexicanHat1D as _RickerWavelet1D
 else:
     from astropy.modeling.core import CompoundModel
+    from astropy.modeling.models import RickerWavelet1D as _RickerWavelet1D
 
 __all__ = ['BlackBody1D', 'BlackBodyNorm1D', 'Box1D', 'ConstFlux1D',
            'Empirical1D', 'Gaussian1D', 'GaussianAbsorption1D',
-           'GaussianFlux1D', 'Lorentz1D', 'MexicanHat1D', 'PowerLawFlux1D',
-           'Trapezoid1D', 'get_waveset', 'get_metadata']
+           'GaussianFlux1D', 'Lorentz1D', 'MexicanHat1D', 'RickerWavelet1D',
+           'PowerLawFlux1D', 'Trapezoid1D', 'get_waveset', 'get_metadata']
 
 
 class BlackBody1D(Fittable1DModel):
@@ -117,10 +119,7 @@ class BlackBody1D(Fittable1DModel):
             Blackbody radiation in PHOTLAM per steradian.
 
         """
-        if ASTROPY_LT_2_0:
-            from astropy.analytic_functions.blackbody import blackbody_nu
-        else:
-            from astropy.modeling.blackbody import blackbody_nu
+        from synphot.blackbody import blackbody_nu
 
         # Silence Numpy
         old_np_err_cfg = np.seterr(all='ignore')
@@ -554,9 +553,9 @@ class Lorentz1D(_models.Lorentz1D):
         return np.asarray(w)
 
 
-class MexicanHat1D(_models.MexicanHat1D):
-    """Same as `astropy.modeling.models.MexicanHat1D`, except with
-    ``sampleset`` defined.
+class RickerWavelet1D(_RickerWavelet1D):
+    """Same as `astropy.modeling.functional_models.RickerWavelet1D`, except
+    with ``sampleset`` defined.
 
     """
     def sampleset(self, factor_step=0.1, **kwargs):
@@ -583,6 +582,12 @@ class MexicanHat1D(_models.MexicanHat1D):
         return np.asarray(w)
 
 
+# TODO: Emit proper deprecation warning.
+# https://github.com/spacetelescope/synphot_refactor/issues/249
+class MexicanHat1D(RickerWavelet1D):
+    """This is the deprecated name for `RickerWavelet1D`."""
+
+
 class PowerLawFlux1D(_models.PowerLaw1D):
     """One dimensional power law model with proper flux handling.
 
diff --git a/synphot/spectrum.py b/synphot/spectrum.py
index 3c70218..fcb464c 100644
--- a/synphot/spectrum.py
+++ b/synphot/spectrum.py
@@ -98,7 +98,10 @@ class BaseSpectrum(object):
             'alpha': u.dimensionless_unscaled,
             'beta': u.dimensionless_unscaled},
         'Lorentz1D': {'amplitude': 'flux', 'x_0': 'wave', 'fwhm': 'wave'},
-        'MexicanHat1D': {'amplitude': 'flux', 'x_0': 'wave', 'sigma': 'wave'},
+        'RickerWavelet1D': {
+            'amplitude': 'flux', 'x_0': 'wave', 'sigma': 'wave'},
+        'MexicanHat1D': {
+            'amplitude': 'flux', 'x_0': 'wave', 'sigma': 'wave'},
         'PowerLaw1D': {
             'amplitude': 'flux', 'x_0': 'wave',
             'alpha': u.dimensionless_unscaled},
@@ -120,6 +123,7 @@ class BaseSpectrum(object):
         'GaussianFlux1D': 'mean',
         'LogParabola1D': 'x_0',
         'Lorentz1D': 'x_0',
+        'RickerWavelet1D': 'x_0',
         'MexicanHat1D': 'x_0',
         'PowerLaw1D': 'x_0',
         'Trapezoid1D': 'x_0'}


Test patch:
diff --git a/synphot/tests/test_blackbody.py b/synphot/tests/test_blackbody.py
new file mode 100644
index 0000000..df1aaa2
--- /dev/null
+++ b/synphot/tests/test_blackbody.py
@@ -0,0 +1,28 @@
+"""This should be updated if blackbody.py is removed from synphot."""
+
+import numpy as np
+from astropy import constants as const
+from astropy import units as u
+
+from ..blackbody import blackbody_nu
+from ..units import FNU
+
+
+# This test was removed from astropy in
+# https://github.com/astropy/astropy/pull/9282
+# but the rest of blackbody.py tests are still over there on astropy.
+def test_blackbody_synphot():
+    """Test that it is consistent with IRAF SYNPHOT BBFUNC."""
+    # Solid angle of solar radius at 1 kpc
+    fac = np.pi * (const.R_sun / const.kpc) ** 2 * u.sr
+
+    with np.errstate(all='ignore'):
+        flux = blackbody_nu([100, 1, 1000, 1e4, 1e5] * u.AA, 5000) * fac
+    assert flux.unit == FNU
+
+    # Special check for overflow value (SYNPHOT gives 0)
+    assert np.log10(flux[0].value) < -143
+
+    np.testing.assert_allclose(
+        flux.value[1:], [0, 2.01950807e-34, 3.78584515e-26, 1.90431881e-27],
+        rtol=0.01)  # 1% accuracy
diff --git a/synphot/tests/test_models.py b/synphot/tests/test_models.py
index 15286aa..5f88dfb 100644
--- a/synphot/tests/test_models.py
+++ b/synphot/tests/test_models.py
@@ -19,7 +19,6 @@ import numpy as np
 import pytest
 
 # ASTROPY
-import astropy
 from astropy import units as u
 from astropy.modeling.models import Const1D
 from astropy.utils import minversion
@@ -39,28 +38,24 @@ else:
     HAS_SCIPY = True
 
 HAS_SCIPY = HAS_SCIPY and minversion(scipy, '0.14')
-ASTROPY_LT_20 = not minversion(astropy, '2.0')
 
 
 def setup_module(module):
-    # https://github.com/astropy/astropy/issues/6383
-    if not ASTROPY_LT_20:
-        import astropy.constants as const
-        from astropy.constants import si, astropyconst13
+    import astropy.constants as const
+    from astropy.constants import si, astropyconst13
 
-        const.sigma_sb = si.sigma_sb = astropyconst13.sigma_sb
-        const.h = si.h = astropyconst13.h
-        const.k_B = si.k_B = astropyconst13.k_B
+    const.sigma_sb = si.sigma_sb = astropyconst13.sigma_sb
+    const.h = si.h = astropyconst13.h
+    const.k_B = si.k_B = astropyconst13.k_B
 
 
 def teardown_module(module):
-    # https://github.com/astropy/astropy/issues/6383
-    if not ASTROPY_LT_20:
-        import astropy.constants as const
-        from astropy.constants import si, astropyconst20
-        const.sigma_sb = si.sigma_sb = astropyconst20.sigma_sb
-        const.h = si.h = astropyconst20.h
-        const.k_B = si.k_B = astropyconst20.k_B
+    import astropy.constants as const
+    from astropy.constants import si, astropyconst20
+
+    const.sigma_sb = si.sigma_sb = astropyconst20.sigma_sb
+    const.h = si.h = astropyconst20.h
+    const.k_B = si.k_B = astropyconst20.k_B
 
 
 class TestBlackBody1D(object):
diff --git a/synphot/tests/test_spectrum.py b/synphot/tests/test_spectrum.py
index c0cc99f..6664fca 100644
--- a/synphot/tests/test_spectrum.py
+++ b/synphot/tests/test_spectrum.py
@@ -12,7 +12,6 @@ import numpy as np
 import pytest
 
 # ASTROPY
-import astropy
 from astropy import modeling
 from astropy import units as u
 from astropy.io import fits
@@ -30,7 +29,7 @@ from .. import exceptions, units
 from ..compat import ASTROPY_LT_4_0, HAS_SPECUTILS  # noqa
 from ..models import (
     BlackBodyNorm1D, Box1D, ConstFlux1D, Empirical1D, Gaussian1D,
-    GaussianAbsorption1D, GaussianFlux1D, Lorentz1D, MexicanHat1D,
+    GaussianAbsorption1D, GaussianFlux1D, Lorentz1D, RickerWavelet1D,
     PowerLawFlux1D, get_waveset)
 from ..observation import Observation
 from ..spectrum import SourceSpectrum, SpectralElement
@@ -43,31 +42,26 @@ else:
     HAS_SCIPY = True
 
 HAS_SCIPY = HAS_SCIPY and minversion(scipy, '0.14')
-ASTROPY_LT_20 = not minversion(astropy, '2.0')
 
 # GLOBAL VARIABLES
 _vspec = None  # Loaded in test_load_vspec()
 
 
 def setup_module(module):
-    # https://github.com/astropy/astropy/issues/6383
-    if not ASTROPY_LT_20:
-        import astropy.constants as const
-        from astropy.constants import si, astropyconst13
+    import astropy.constants as const
+    from astropy.constants import si, astropyconst13
 
-        const.sigma_sb = si.sigma_sb = astropyconst13.sigma_sb
-        const.h = si.h = astropyconst13.h
-        const.k_B = si.k_B = astropyconst13.k_B
+    const.sigma_sb = si.sigma_sb = astropyconst13.sigma_sb
+    const.h = si.h = astropyconst13.h
+    const.k_B = si.k_B = astropyconst13.k_B
 
 
 def teardown_module(module):
-    # https://github.com/astropy/astropy/issues/6383
-    if not ASTROPY_LT_20:
-        import astropy.constants as const
-        from astropy.constants import si, astropyconst20
-        const.sigma_sb = si.sigma_sb = astropyconst20.sigma_sb
-        const.h = si.h = astropyconst20.h
-        const.k_B = si.k_B = astropyconst20.k_B
+    import astropy.constants as const
+    from astropy.constants import si, astropyconst20
+    const.sigma_sb = si.sigma_sb = astropyconst20.sigma_sb
+    const.h = si.h = astropyconst20.h
+    const.k_B = si.k_B = astropyconst20.k_B
 
 
 @pytest.mark.remote_data
@@ -517,8 +511,8 @@ class TestBuildModels:
         np.testing.assert_allclose(
             y.value, [0.00249377, 1, 0.00249377], rtol=1e-5)
 
-    def test_MexicanHat1D(self):
-        sp = SourceSpectrum(MexicanHat1D, amplitude=1, x_0=6000, sigma=100)
+    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])
diff --git a/synphot/tests/test_thermal.py b/synphot/tests/test_thermal.py
index bd6edd8..a5c0fea 100644
--- a/synphot/tests/test_thermal.py
+++ b/synphot/tests/test_thermal.py
@@ -9,7 +9,6 @@ import numpy as np
 import pytest
 
 # ASTROPY
-import astropy
 from astropy import units as u
 from astropy.utils import minversion
 from astropy.utils.data import get_pkg_data_filename
@@ -26,28 +25,24 @@ else:
     HAS_SCIPY = True
 
 HAS_SCIPY = HAS_SCIPY and minversion(scipy, '0.14')
-ASTROPY_LT_20 = not minversion(astropy, '2.0')
 
 
 def setup_module(module):
-    # https://github.com/astropy/astropy/issues/6383
-    if not ASTROPY_LT_20:
-        import astropy.constants as const
-        from astropy.constants import si, astropyconst13
+    import astropy.constants as const
+    from astropy.constants import si, astropyconst13
 
-        const.sigma_sb = si.sigma_sb = astropyconst13.sigma_sb
-        const.h = si.h = astropyconst13.h
-        const.k_B = si.k_B = astropyconst13.k_B
+    const.sigma_sb = si.sigma_sb = astropyconst13.sigma_sb
+    const.h = si.h = astropyconst13.h
+    const.k_B = si.k_B = astropyconst13.k_B
 
 
 def teardown_module(module):
-    # https://github.com/astropy/astropy/issues/6383
-    if not ASTROPY_LT_20:
-        import astropy.constants as const
-        from astropy.constants import si, astropyconst20
-        const.sigma_sb = si.sigma_sb = astropyconst20.sigma_sb
-        const.h = si.h = astropyconst20.h
-        const.k_B = si.k_B = astropyconst20.k_B
+    import astropy.constants as const
+    from astropy.constants import si, astropyconst20
+
+    const.sigma_sb = si.sigma_sb = astropyconst20.sigma_sb
+    const.h = si.h = astropyconst20.h
+    const.k_B = si.k_B = astropyconst20.k_B
 
 
 @pytest.mark.skipif('not HAS_SCIPY')