diff --git a/CHANGES.rst b/CHANGES.rst
index 0e44bf9..eba7947 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -1,6 +1,15 @@
-1.3.1 (unreleased)
+1.4.0 (unreleased)
 ==================
 
+- ``read_fits_spec()`` now uses ``astropy.table.QTable.read`` for parsing to
+  ensure that the correct ``TUNITn`` is read. As a result, ``wave_unit`` and
+  ``flux_unit`` keywords are deprecated and no longer used in that function.
+  Additionally, if any ``TUNITn`` in the table is invalid, regardless whether
+  the column is used or not, an exception will now be raised. [#384]
+
+- ``read_spec()`` now detects whether given filename is FITS more consistently
+  w.r.t. ``astropy``. [#384]
+
 1.3.0 (2023-11-28)
 ==================
 
diff --git a/docs/synphot/overview.rst b/docs/synphot/overview.rst
index 0bf3455..6b43962 100644
--- a/docs/synphot/overview.rst
+++ b/docs/synphot/overview.rst
@@ -396,9 +396,7 @@ keywords (unless you overwrite them with non-default values in
 * ``TTYPE2`` set to "FLUX" (for source spectrum) or "THROUGHPUT"
   (for bandpass).
 
-While these were required in ASTROLIB PYSYNPHOT, they are optional here in that
-default units would be applied, where applicable, if they are missing from
-the header. Regardless, setting them is highly recommended:
+These were required in ASTROLIB PYSYNPHOT, as well as this package:
 
 * ``TUNIT1`` set to :ref:`supported wavelength unit name <synphot_units>`.
 * ``TUNIT2`` set to :ref:`supported flux unit name <synphot_units>`
diff --git a/synphot/reddening.py b/synphot/reddening.py
index 6ee9b9f..027e1ee 100644
--- a/synphot/reddening.py
+++ b/synphot/reddening.py
@@ -7,6 +7,7 @@
 # THIRD-PARTY
 import numpy as np
 from astropy import units as u
+from astropy.io.fits.connect import is_fits
 
 # LOCAL
 from synphot import exceptions, specio, units
@@ -137,8 +138,8 @@ def to_fits(self, filename, wavelengths=None, **kwargs):
     def from_file(cls, filename, **kwargs):
         """Create a reddening law from file.
 
-        If filename has 'fits' or 'fit' suffix, it is read as FITS.
-        Otherwise, it is read as ASCII.
+        If filename is recognized by ``astropy.io.fits`` as FITS,
+        it is read as such. Otherwise, it is read as ASCII.
 
         Parameters
         ----------
@@ -156,13 +157,12 @@ def from_file(cls, filename, **kwargs):
             Empirical reddening law.
 
         """
-        if 'flux_unit' not in kwargs:
+        if is_fits("", filename, None):
+            if 'flux_col' not in kwargs:
+                kwargs['flux_col'] = 'Av/E(B-V)'
+        elif 'flux_unit' not in kwargs:  # pragma: no cover
             kwargs['flux_unit'] = cls._internal_flux_unit
 
-        if ((filename.endswith('fits') or filename.endswith('fit')) and
-                'flux_col' not in kwargs):
-            kwargs['flux_col'] = 'Av/E(B-V)'
-
         header, wavelengths, rvs = specio.read_spec(filename, **kwargs)
 
         return cls(Empirical1D, points=wavelengths, lookup_table=rvs,
@@ -217,13 +217,12 @@ def from_extinction_model(cls, modelname, **kwargs):
 
         filename = cfgitem()
 
-        if 'flux_unit' not in kwargs:
+        if is_fits("", filename, None):
+            if 'flux_col' not in kwargs:
+                kwargs['flux_col'] = 'Av/E(B-V)'
+        elif 'flux_unit' not in kwargs:  # pragma: no cover
             kwargs['flux_unit'] = cls._internal_flux_unit
 
-        if ((filename.endswith('fits') or filename.endswith('fit')) and
-                'flux_col' not in kwargs):
-            kwargs['flux_col'] = 'Av/E(B-V)'
-
         header, wavelengths, rvs = specio.read_remote_spec(filename, **kwargs)
         header['filename'] = filename
         header['descrip'] = cfgitem.description
diff --git a/synphot/specio.py b/synphot/specio.py
index 91e2a45..83097b9 100644
--- a/synphot/specio.py
+++ b/synphot/specio.py
@@ -12,7 +12,10 @@
 from astropy import log
 from astropy import units as u
 from astropy.io import ascii, fits
+from astropy.io.fits.connect import is_fits
+from astropy.table import QTable
 from astropy.utils.data import get_readable_fileobj
+from astropy.utils.decorators import deprecated_renamed_argument
 from astropy.utils.exceptions import AstropyUserWarning
 
 # LOCAL
@@ -88,7 +91,7 @@ def read_spec(filename, fname='', **kwargs):
     elif not fname:  # pragma: no cover
         raise exceptions.SynphotError('Cannot determine filename.')
 
-    if fname.endswith('fits') or fname.endswith('fit'):
+    if is_fits("", fname, None):
         read_func = read_fits_spec
     else:
         read_func = read_ascii_spec
@@ -143,12 +146,15 @@ def read_ascii_spec(filename, wave_unit=u.AA, flux_unit=units.FLAM, **kwargs):
     return header, wavelengths, fluxes
 
 
+@deprecated_renamed_argument(
+    ["wave_unit", "flux_unit"], [None, None], ["1.4", "1.4"],
+    alternative='TUNITn as per FITS standards')
 def read_fits_spec(filename, ext=1, wave_col='WAVELENGTH', flux_col='FLUX',
                    wave_unit=u.AA, flux_unit=units.FLAM):
     """Read FITS spectrum.
 
-    Wavelength and flux units are extracted from ``TUNIT1`` and ``TUNIT2``
-    keywords, respectively, from data table (not primary) header.
+    Wavelength and flux units are extracted from respective ``TUNITn``
+    keywords, from data table (not primary) header.
     If these keywords are not present, units are taken from
     ``wave_unit`` and ``flux_unit`` instead.
 
@@ -164,9 +170,11 @@ def read_fits_spec(filename, ext=1, wave_col='WAVELENGTH', flux_col='FLUX',
         Wavelength and flux column names (case-insensitive).
 
     wave_unit, flux_unit : str or `~astropy.units.Unit`
-        Wavelength and flux units, which default to Angstrom and FLAM,
-        respectively. These are *only* used if ``TUNIT1`` and ``TUNIT2``
-        keywords are not present in table (not primary) header.
+        Wavelength and flux units. These are *no longer used*.
+        Define your units in the respective ``TUNITn``
+        keywords in table (not primary) header.
+
+        .. deprecated:: 1.4
 
     Returns
     -------
@@ -177,37 +185,42 @@ def read_fits_spec(filename, ext=1, wave_col='WAVELENGTH', flux_col='FLUX',
         Wavelength and flux of the spectrum.
 
     """
+    wave_col = wave_col.lower()
+    flux_col = flux_col.lower()
+
     try:
         fs = fits.open(filename)
-        header = dict(fs[str('PRIMARY')].header)
-        wave_dat = fs[ext].data.field(wave_col).copy()
-        flux_dat = fs[ext].data.field(flux_col).copy()
-        fits_wave_unit = fs[ext].header.get('TUNIT1')
-        fits_flux_unit = fs[ext].header.get('TUNIT2')
-
-        if fits_wave_unit is not None:
-            try:
-                wave_unit = units.validate_unit(fits_wave_unit)
-            except (exceptions.SynphotError, ValueError) as e:  # pragma: no cover  # noqa: E501
-                warnings.warn(
-                    '{0} from FITS header is not valid wavelength unit, using '
-                    '{1}: {2}'.format(fits_wave_unit, wave_unit, e),
-                    AstropyUserWarning)
-
-        if fits_flux_unit is not None:
-            try:
-                flux_unit = units.validate_unit(fits_flux_unit)
-            except (exceptions.SynphotError, ValueError) as e:  # pragma: no cover  # noqa: E501
-                warnings.warn(
-                    '{0} from FITS header is not valid flux unit, using '
-                    '{1}: {2}'.format(fits_flux_unit, flux_unit, e),
-                    AstropyUserWarning)
-
-        wave_unit = units.validate_unit(wave_unit)
-        flux_unit = units.validate_unit(flux_unit)
-
-        wavelengths = wave_dat * wave_unit
-        fluxes = flux_dat * flux_unit
+        subhdu = fs[ext]
+
+        # Need to fix table units
+        for key in subhdu.header["TUNIT*"]:
+            val = subhdu.header[key]
+            if not val:
+                continue
+            newval = units.validate_unit(val)
+            subhdu.header[key] = newval.to_string()  # Must be generic to handle mag  # noqa: E501
+
+        with warnings.catch_warnings():
+            warnings.filterwarnings("ignore", category=u.UnitsWarning,
+                                    message=".* did not parse as fits unit")
+            t = QTable.read(subhdu)
+        header = dict(fs["PRIMARY"].header)
+
+        # https://github.com/astropy/astropy/issues/16221
+        lower_colnames = [c.lower() for c in t.colnames]
+        t_col_wave = t.columns[lower_colnames.index(wave_col)]
+        if t_col_wave.unit:
+            t_col_wave_unit = units.validate_unit(t_col_wave.unit.to_string())
+        else:
+            t_col_wave_unit = u.dimensionless_unscaled
+        t_col_flux = t.columns[lower_colnames.index(flux_col)]
+        if t_col_flux.unit:
+            t_col_flux_unit = units.validate_unit(t_col_flux.unit.to_string())
+        else:
+            t_col_flux_unit = u.dimensionless_unscaled
+        wavelengths = t_col_wave.value * t_col_wave_unit
+        fluxes = t_col_flux.value * t_col_flux_unit
+
     finally:
         if isinstance(filename, str):
             fs.close()
diff --git a/synphot/spectrum.py b/synphot/spectrum.py
index f1a6173..196649e 100644
--- a/synphot/spectrum.py
+++ b/synphot/spectrum.py
@@ -14,6 +14,7 @@
 # ASTROPY
 from astropy import log
 from astropy import units as u
+from astropy.io.fits.connect import is_fits
 from astropy.modeling import Model
 from astropy.modeling.core import CompoundModel
 from astropy.modeling.models import RedshiftScaleFactor, Scale
@@ -1921,8 +1922,8 @@ def to_fits(self, filename, wavelengths=None, **kwargs):
     def from_file(cls, filename, **kwargs):
         """Creates a bandpass from file.
 
-        If filename has 'fits' or 'fit' suffix, it is read as FITS.
-        Otherwise, it is read as ASCII.
+        If filename is recognized by ``astropy.io.fits`` as FITS,
+        it is read as such. Otherwise, it is read as ASCII.
 
         Parameters
         ----------
@@ -1940,13 +1941,12 @@ def from_file(cls, filename, **kwargs):
             Empirical bandpass.
 
         """
-        if 'flux_unit' not in kwargs:
+        if is_fits("", filename, None):
+            if 'flux_col' not in kwargs:
+                kwargs['flux_col'] = 'THROUGHPUT'
+        elif 'flux_unit' not in kwargs:  # pragma: no cover
             kwargs['flux_unit'] = cls._internal_flux_unit
 
-        if ((filename.endswith('fits') or filename.endswith('fit')) and
-                'flux_col' not in kwargs):
-            kwargs['flux_col'] = 'THROUGHPUT'
-
         header, wavelengths, throughput = specio.read_spec(filename, **kwargs)
         return cls(Empirical1D, points=wavelengths, lookup_table=throughput,
                    keep_neg=True, meta={'header': header})
@@ -2009,13 +2009,12 @@ def from_filter(cls, filtername, **kwargs):
 
         filename = cfgitem()
 
-        if 'flux_unit' not in kwargs:
+        if is_fits("", filename, None):
+            if 'flux_col' not in kwargs:
+                kwargs['flux_col'] = 'THROUGHPUT'
+        elif 'flux_unit' not in kwargs:  # pragma: no cover
             kwargs['flux_unit'] = cls._internal_flux_unit
 
-        if ((filename.endswith('fits') or filename.endswith('fit')) and
-                'flux_col' not in kwargs):
-            kwargs['flux_col'] = 'THROUGHPUT'
-
         header, wavelengths, throughput = specio.read_remote_spec(
             filename, **kwargs)
         header['filename'] = filename
diff --git a/synphot/tests/test_binning.py b/synphot/tests/test_binning.py
index 79dd482..1815342 100644
--- a/synphot/tests/test_binning.py
+++ b/synphot/tests/test_binning.py
@@ -198,7 +198,7 @@ def setup_class(self):
             get_pkg_data_filename(
                 os.path.join('data', 'hst_acs_hrc_f555w.fits'),
                 package='synphot.tests'),
-            flux_col='THROUGHPUT', flux_unit=u.dimensionless_unscaled)
+            flux_col='THROUGHPUT')
 
         # Binned data.
         bins = generate_wavelengths(
diff --git a/synphot/tests/test_reddening.py b/synphot/tests/test_reddening.py
index 336d07c..1343240 100644
--- a/synphot/tests/test_reddening.py
+++ b/synphot/tests/test_reddening.py
@@ -3,8 +3,6 @@
 
 # STDLIB
 import os
-import shutil
-import tempfile
 
 # THIRD-PARTY
 import numpy as np
@@ -156,32 +154,26 @@ def test_redlaw_from_model_exception():
         ReddeningLaw.from_extinction_model('foo')
 
 
-class TestWriteReddeningLaw:
+@pytest.mark.parametrize('ext_hdr', [None, {'foo': 'foo'}])
+def test_write_reddening_law(tmp_path, ext_hdr):
     """Test ReddeningLaw ``to_fits()`` method."""
-    def setup_class(self):
-        self.outdir = tempfile.mkdtemp()
-        self.x = np.linspace(1000, 5000, 5)
-        self.y = np.linspace(1, 5, 5) * 0.1
-        self.redlaw = ReddeningLaw(
-            Empirical1D, points=self.x, lookup_table=self.y)
-
-    @pytest.mark.parametrize('ext_hdr', [None, {'foo': 'foo'}])
-    def test_write(self, ext_hdr):
-        outfile = os.path.join(self.outdir, 'outredlaw.fits')
-
-        if ext_hdr is None:
-            self.redlaw.to_fits(outfile, overwrite=True)
-        else:
-            self.redlaw.to_fits(outfile, overwrite=True, ext_header=ext_hdr)
-
-        # Read it back in and check
-        redlaw2 = ReddeningLaw.from_file(outfile)
-        np.testing.assert_allclose(redlaw2.waveset.value, self.x)
-        np.testing.assert_allclose(redlaw2(self.x).value, self.y)
-
-        if ext_hdr is not None:
-            hdr = fits.getheader(outfile, 1)
-            assert 'foo' in hdr
-
-    def teardown_class(self):
-        shutil.rmtree(self.outdir)
+    x = np.linspace(1000, 5000, 5)
+    y = np.linspace(1, 5, 5) * 0.1
+    redlaw = ReddeningLaw(
+        Empirical1D, points=x, lookup_table=y, meta={"expr": "ebv(test)"})
+
+    outfile = str(tmp_path / 'outredlaw.fits')
+
+    if ext_hdr is None:
+        redlaw.to_fits(outfile, overwrite=True)
+    else:
+        redlaw.to_fits(outfile, overwrite=True, ext_header=ext_hdr)
+
+    # Read it back in and check
+    redlaw2 = ReddeningLaw.from_file(outfile)
+    np.testing.assert_allclose(redlaw2.waveset.value, x)
+    np.testing.assert_allclose(redlaw2(x).value, y)
+
+    if ext_hdr is not None:
+        hdr = fits.getheader(outfile, 1)
+        assert 'foo' in hdr
diff --git a/synphot/tests/test_specio.py b/synphot/tests/test_specio.py
index 7a4e56c..e807ca3 100644
--- a/synphot/tests/test_specio.py
+++ b/synphot/tests/test_specio.py
@@ -3,8 +3,6 @@
 
 # STDLIB
 import os
-import shutil
-import tempfile
 
 # THIRD-PARTY
 import numpy as np
@@ -15,10 +13,12 @@
 from astropy.io import fits
 from astropy.tests.helper import assert_quantity_allclose
 from astropy.utils.data import get_pkg_data_filename
-from astropy.utils.exceptions import AstropyUserWarning
+from astropy.utils.exceptions import (
+    AstropyUserWarning, AstropyDeprecationWarning)
 
 # LOCAL
 from synphot import exceptions, specio, units
+from synphot.spectrum import SpectralElement
 
 
 @pytest.mark.remote_data
@@ -55,7 +55,6 @@ class TestReadWriteFITS:
     """Test read/write FITS spectrum."""
     def setup_class(self):
         self.epsilon = 0.00031
-        self.outdir = tempfile.mkdtemp()
         self.wave = np.array([1000.0, 2000.0, 2000.0 + self.epsilon, 3000.0,
                               4000.0, 5000.0], dtype=np.float64) * u.AA
         self.flux = np.array([0.1, 100.2, 10.0, 0.0, 6.5, 1.2],
@@ -63,9 +62,9 @@ def setup_class(self):
         self.prihdr = {'PEDIGREE': 'DUMMY'}
         self.scihdr = {'SPEC_SRC': 'RANDOM'}
 
-    def test_array_data(self):
+    def test_array_data(self, tmp_path):
         """Data as Numpy array."""
-        outfile = os.path.join(self.outdir, 'outspec1.fits')
+        outfile = str(tmp_path / 'outspec1.fits')
 
         # Write it out
         with pytest.warns(AstropyUserWarning, match=r'rows are thrown out'):
@@ -76,7 +75,8 @@ def test_array_data(self):
                 wave_unit=self.wave.unit, flux_unit=self.flux.unit)
 
         # Read it back in and check values (flux_unit should be ignored)
-        hdr, wave, flux = specio.read_spec(outfile, flux_unit='foo')
+        with pytest.warns(AstropyDeprecationWarning, match=r"\"flux_unit\" was deprecated"):  # noqa: E501
+            hdr, wave, flux = specio.read_spec(outfile, flux_unit='foo')
 
         # Compare data
         np.testing.assert_allclose(
@@ -95,9 +95,9 @@ def test_array_data(self):
         assert sci_hdr['SPEC_SRC'] == 'RANDOM'
         assert sci_hdr['TFORM2'].lower() == 'e'
 
-    def test_quantity_data(self):
+    def test_quantity_data(self, tmp_path):
         """Data as Quantity."""
-        outfile = os.path.join(self.outdir, 'outspec2.fits')
+        outfile = str(tmp_path / 'outspec2.fits')
 
         # Write it out (flux_unit should be ignored)
         specio.write_fits_spec(
@@ -105,7 +105,8 @@ def test_quantity_data(self):
             ext_header=self.scihdr, precision='double', flux_unit='foo')
 
         # Read it back in and check values (flux_unit should be ignored)
-        hdr, wave, flux = specio.read_spec(outfile, flux_unit='foo')
+        with pytest.warns(AstropyDeprecationWarning, match=r"\"flux_unit\" was deprecated"):  # noqa: E501
+            hdr, wave, flux = specio.read_spec(outfile, flux_unit='foo')
 
         # Compare data (trim_zero=True, pad_zero_ends=True)
         np.testing.assert_allclose(
@@ -125,9 +126,9 @@ def test_quantity_data(self):
         assert sci_hdr['SPEC_SRC'] == 'RANDOM'
         assert sci_hdr['TFORM2'].lower() == 'd'
 
-    def test_exceptions(self):
+    def test_exceptions(self, tmp_path):
         """Test for appropriate exceptions."""
-        outfile = os.path.join(self.outdir, 'outspec3.fits')
+        outfile = str(tmp_path / 'outspec3.fits')
 
         # Shape mismatch
         with pytest.raises(exceptions.SynphotError):
@@ -149,5 +150,61 @@ def test_exceptions(self):
             specio.write_fits_spec(
                 outfile, self.wave, np.arange(6), overwrite=True)
 
-    def teardown_class(self):
-        shutil.rmtree(self.outdir)
+
+def test_read_nonstandard_fits_cols_01(tmp_path):
+    """See https://github.com/spacetelescope/synphot_refactor/issues/372"""
+    pix = np.arange(5, dtype=int) + 1
+    wav = (pix * 0.1) * u.micron
+    trace = np.array([0, 0.5, 1, 0.9, 0])
+    coldefs = fits.ColDefs([
+        fits.Column(name="X", format="I", array=pix),
+        fits.Column(name="Wavelength", format="E",
+                    unit=wav.unit.to_string(format="fits"), array=wav.value),
+        fits.Column(name="Trace", format="E", array=trace)])
+    hdulist = fits.HDUList([
+        fits.PrimaryHDU(),
+        fits.BinTableHDU.from_columns(coldefs)])
+    outfile = str(tmp_path / "jwst_niriss_soss_trace.fits")
+    hdulist.writeto(outfile, overwrite=True)
+
+    # Make sure column names are still case insensitive.
+    for (wave_col, flux_col) in (
+            ("Wavelength", "Trace"),
+            ("WAVELENGTH", "TRACE"),
+            ("wavelength", "trace")):
+        tr = SpectralElement.from_file(
+            outfile, wave_col=wave_col, flux_col=flux_col)
+        assert_quantity_allclose(tr.waveset, wav)
+        assert_quantity_allclose(tr(wav), trace, atol=1e-7)
+
+
+def test_read_nonstandard_fits_cols_02(tmp_path):
+    """See https://github.com/spacetelescope/synphot_refactor/issues/372"""
+
+    wav = (np.arange(5) + 1) * u.nm
+    flux_unit_str = "ph/s/m2/micron/arcsec2"  # Invalid but should not matter.
+    flux = np.ones(5)
+    thru = np.array([0, 0.5, 1, 0.9, 0])
+    coldefs = fits.ColDefs([
+        fits.Column(name="lam", format="E",
+                    unit=wav.unit.to_string(format="fits"), array=wav.value),
+        fits.Column(name="flux", format="E",
+                    unit=flux_unit_str, array=flux),
+        fits.Column(name="dflux1", format="E",
+                    unit=flux_unit_str, array=flux),
+        fits.Column(name="dflux2", format="E",
+                    unit=flux_unit_str, array=flux),
+        fits.Column(name="trans", format="E", unit="1", array=thru)])
+    hdulist = fits.HDUList([
+        fits.PrimaryHDU(),
+        fits.BinTableHDU.from_columns(coldefs)])
+    outfile = str(tmp_path / "skytable.fits")
+    hdulist.writeto(outfile, overwrite=True)
+
+    with pytest.warns(u.UnitsWarning, match="'ph/s/m2/micron/arcsec2'"):  # noqa: E501
+        header, wavelengths, transmission = specio.read_spec(
+            outfile, wave_col="lam", flux_col="trans")
+
+    assert header["SIMPLE"]
+    assert_quantity_allclose(wavelengths, wav)
+    assert_quantity_allclose(transmission, thru)
diff --git a/synphot/tests/test_spectrum_misc.py b/synphot/tests/test_spectrum_misc.py
index fc93742..c4a8503 100644
--- a/synphot/tests/test_spectrum_misc.py
+++ b/synphot/tests/test_spectrum_misc.py
@@ -2,11 +2,6 @@
 """Test spectrum.py module and related functionalities that are not covered
 by ``test_spectrum_source.py`` nor ``test_spectrum_bandpass.py``."""
 
-# STDLIB
-import os
-import shutil
-import tempfile
-
 # THIRD-PARTY
 import numpy as np
 import pytest
@@ -390,7 +385,6 @@ class DummyObject:
 class TestWriteSpec:
     """Test spectrum to_fits() method."""
     def setup_class(self):
-        self.outdir = tempfile.mkdtemp()
         self.sp = SourceSpectrum(
             Empirical1D, points=_wave, lookup_table=_flux_photlam,
             meta={'expr': 'Test source'})
@@ -404,8 +398,8 @@ def setup_class(self):
          (True, {'foo': 'foo'}),
          (False, None),
          (False, {'foo': 'foo'})])
-    def test_write(self, is_sp, ext_hdr):
-        outfile = os.path.join(self.outdir, 'outspec.fits')
+    def test_write(self, tmp_path, is_sp, ext_hdr):
+        outfile = str(tmp_path / 'outspec.fits')
 
         if is_sp:
             sp1 = self.sp
@@ -426,6 +420,3 @@ def test_write(self, is_sp, ext_hdr):
         assert 'expr' in hdr
         if ext_hdr is not None:
             assert 'foo' in hdr
-
-    def teardown_class(self):
-        shutil.rmtree(self.outdir)
diff --git a/synphot/tests/test_utils.py b/synphot/tests/test_utils.py
index 6435647..13f42f0 100644
--- a/synphot/tests/test_utils.py
+++ b/synphot/tests/test_utils.py
@@ -116,25 +116,25 @@ def test_merge_same(self):
         np.testing.assert_array_equal(wave, self.wave)
 
 
-def test_download_bad_root(tmpdir):
+def test_download_bad_root(tmp_path):
     """Test data download helper when input dir is invalid."""
-    ptr = tmpdir.join('bad_cdbs')
-    ptr.write('content')
+    ptr = tmp_path / 'bad_cdbs'
+    ptr.write_text("something")
     cdbs_root = str(ptr)
 
-    with pytest.raises(OSError):
+    with pytest.raises(OSError, match=".* must be a directory"):
         utils.download_data(cdbs_root, verbose=False)
 
     with pytest.raises(FileNotFoundError):
         utils.download_data('', verbose=False)
 
 
-def test_download_data(tmpdir):
+def test_download_data(tmp_path):
     """Test data download helper in dry run mode."""
     from synphot.config import conf
 
     # Use case where user downloads all data into new dir.
-    cdbs_root = os.path.join(tmpdir.strpath, 'cdbs')
+    cdbs_root = str(tmp_path / 'cdbs')
     file_list_1 = utils.download_data(cdbs_root, verbose=False, dry_run=True)
     filename = file_list_1[0]
     assert len(file_list_1) == 21
diff --git a/synphot/thermal.py b/synphot/thermal.py
index 093d412..8ad4ff7 100644
--- a/synphot/thermal.py
+++ b/synphot/thermal.py
@@ -4,6 +4,7 @@
 # ASTROPY
 from astropy import units as u
 from astropy.io import fits
+from astropy.io.fits.connect import is_fits
 
 # LOCAL
 from synphot import exceptions, specio, units
@@ -124,7 +125,7 @@ def from_file(cls, filename, temperature_key='DEFT',
             Invalid inputs.
 
         """
-        if not (filename.endswith('fits') or filename.endswith('fit')):
+        if not is_fits("", filename, None):
             raise exceptions.SynphotError('Only FITS format is supported.')
 
         # Extra info from table header
@@ -138,9 +139,6 @@ def from_file(cls, filename, temperature_key='DEFT',
 
         beam_fill_factor = tab_hdr.get('BEAMFILL', 1)
 
-        if 'flux_unit' not in kwargs:
-            kwargs['flux_unit'] = cls._internal_flux_unit
-
         if 'flux_col' not in kwargs:
             kwargs['flux_col'] = 'EMISSIVITY'
 
diff --git a/synphot/units.py b/synphot/units.py
index e3c4f9f..d9c3ed4 100644
--- a/synphot/units.py
+++ b/synphot/units.py
@@ -291,10 +291,22 @@ def validate_unit(input_unit):
         elif input_unit_lowcase == 'inversemicrons':
             output_unit = u.micron ** -1
         elif input_unit_lowcase in ('transmission', 'extinction',
-                                    'emissivity'):
+                                    'emissivity', 'throughput'):
             output_unit = THROUGHPUT
         elif input_unit_lowcase == 'jy':
             output_unit = u.Jy
+        elif input_unit_lowcase == "flam":
+            output_unit = FLAM
+        elif input_unit_lowcase == "fnu":
+            output_unit = FNU
+        elif input_unit_lowcase == "photlam":
+            output_unit = PHOTLAM
+        elif input_unit_lowcase == "photnu":
+            output_unit = PHOTNU
+        elif input_unit_lowcase == "none":
+            output_unit = u.dimensionless_unscaled
+        elif input_unit_lowcase == "sec":
+            output_unit = u.s
 
         # Work around mag unit limitations
         elif input_unit_lowcase in ('stmag', 'mag(st)'):