Merge pull request #371 from mkelley/calspec-sun-202304

Add CALSPEC model solar spectrum

Author: mkelley

docs/sbpy/calib.rst

@@ -5,14 +5,15 @@ Spectral Standards and Photometric Calibration (`sbpy.calib`)
 
 sbpy's photometric calibration is based on spectra of the Sun and Vega.  For example, they are used to convert between :ref:`reflectance, cross-section, and magnitude <reflectance-equivalencies>`, between :ref:`Afρ and spectral flux density <afrho-to-from-flux-density>`, and between :ref:`Vega-based and other magnitude systems <vega-magnitudes>`.   sbpy has built-in spectra for each, and users may provide their own.
 
-The spectrum of `Bohlin (2014) <https://dx.doi.org/10.1088/0004-6256/147/6/127>`_ is the default and only built-in spectrum for Vega.  It is distributed with sbpy.  Four solar spectra are built-in:
+The spectrum of `Bohlin (2014) <https://dx.doi.org/10.1088/0004-6256/147/6/127>`_ is the default and only built-in spectrum for Vega.  It is distributed with sbpy.  Five solar spectra are built-in:
 
   * Castelli1996 - Castelli model from Colina et al. (1996).
   * E490_2014 - E490 (2014) standard.
   * E490_2014LR - A low resolution version of the E490 standard.
   * Kurucz1993 - Kurucz (1993) model.
+  * calspec - R=5000, created by R. Bohlin from Kurucz Special Model
 
-The E490 spectra are included with sbpy, and the Kurucz and Castelli spectra are downloaded as needed from `STScI's astronomical catalog <https://www.stsci.edu/hst/instrumentation/reference-data-for-calibration-and-tools/astronomical-catalogs>`_.
+The E490 spectra are included with sbpy, and the others are downloaded as needed from MAST's `Spectral Atlas Files for Synphot Software (REFERENCE-ATLASES) <https://archive.stsci.edu/hlsp/reference-atlases>`_ or STScI's `CALSPEC Database <https://www.stsci.edu/hst/instrumentation/reference-data-for-calibration-and-tools/astronomical-catalogs/calspec>`_.
 
 Each star has a class for use within sbpy.  The classes can be initialized with the default spectrum using :func:`~sbpy.calib.SpectralStandard.from_default`:
 
@@ -32,9 +33,10 @@ The names of the built-in sources are stored as an internal array.  They can be
         name                                description
     ------------ -----------------------------------------------------------------
     Castelli1996      Castelli model, scaled and presented by Colina et al. (1996)
-        E490_2014                E490-00a (2014) reference solar spectrum (Table 3)
-    E490_2014LR E490-00a (2014) low resolution reference solar spectrum (Table 4)
-        Kurucz1993               Kurucz (1993) model, scaled by Colina et al. (1996)
+       E490_2014                E490-00a (2014) reference solar spectrum (Table 3)
+     E490_2014LR E490-00a (2014) low resolution reference solar spectrum (Table 4)
+      Kurucz1993               Kurucz (1993) model, scaled by Colina et al. (1996)
+         calspec            R=5000, created by R. Bohlin from Kurucz Special Model
     >>> sun = Sun.from_builtin('E490_2014LR')
     >>> print(sun)
     <Sun: E490-00a (2014) low resolution reference solar spectrum (Table 4)>

sbpy/calib/solar_sources.py

@@ -52,6 +52,13 @@ class SolarSpectra:
         'bibcode': '1996AJ....112..307C'
     }
 
+    calspec = {
+        "filename": ("https://archive.stsci.edu/hlsps/reference-atlases/cdbs/"
+                     "current_calspec/sun_mod_001.fits"),
+        "description": "R=5000, created by R. Bohlin from Kurucz Special Model",
+        "bibcode": "2014PASP..126..711B"
+    }
+
 
 class SolarPhotometry:
     """Built-in solar photometry.

sbpy/calib/tests/test_sun.py

@@ -10,56 +10,64 @@
 
 try:
     import scipy
+
     HAS_SCIPY = True
 except ImportError:
     HAS_SCIPY = False
 
 
 class TestSun:
     def test___repr__(self):
-        with solar_spectrum.set('E490_2014LR'):
-            assert (repr(Sun.from_default()) ==
-                    ('<Sun: E490-00a (2014) low resolution reference '
-                     'solar spectrum (Table 4)>'))
+        with solar_spectrum.set("E490_2014LR"):
+            assert repr(Sun.from_default()) == (
+                "<Sun: E490-00a (2014) low resolution reference "
+                "solar spectrum (Table 4)>"
+            )
 
         sun = Sun.from_array([1, 2] * u.um, [1, 2] * u.Jy)
-        assert repr(sun) == '<Sun>'
+        assert repr(sun) == "<Sun>"
 
     def test_from_builtin(self):
-        sun = Sun.from_builtin('E490_2014LR')
-        assert sun.description == solar_sources.SolarSpectra.E490_2014LR['description']
+        sun = Sun.from_builtin("E490_2014LR")
+        assert (
+            sun.description
+            == solar_sources.SolarSpectra.E490_2014LR["description"]
+        )
 
     def test_from_builtin_unknown(self):
         with pytest.raises(UndefinedSourceError):
-            Sun.from_builtin('not a solar spectrum')
+            Sun.from_builtin("not a solar spectrum")
 
     def test_from_default(self):
-        with solar_spectrum.set('E490_2014LR'):
+        with solar_spectrum.set("E490_2014LR"):
             sun = Sun.from_default()
-            assert sun.description == solar_sources.SolarSpectra.E490_2014LR['description']
+            assert (
+                sun.description
+                == solar_sources.SolarSpectra.E490_2014LR["description"]
+            )
 
     def test_call_single_wavelength(self):
-        with solar_spectrum.set('E490_2014'):
+        with solar_spectrum.set("E490_2014"):
             sun = solar_spectrum.get()
             f = sun(0.5555 * u.um)
             assert np.isclose(f.value, 1897)
 
     def test_call_single_frequency(self):
-        with solar_spectrum.set('E490_2014'):
+        with solar_spectrum.set("E490_2014"):
             sun = solar_spectrum.get()
             f = sun(3e14 * u.Hz)
-            assert np.isclose(f.value, 2.49484251e+14)
+            assert np.isclose(f.value, 2.49484251e14)
 
-    @pytest.mark.skipif('not HAS_SCIPY')
+    @pytest.mark.skipif("not HAS_SCIPY")
     def test_sun_observe_wavelength_array(self):
         from scipy.integrate import trapz
 
-        unit = 'W/(m2 um)'
+        unit = "W/(m2 um)"
 
         # compare Sun's rebinning with an integration over the spectrum
-        sun = Sun.from_builtin('E490_2014')
+        sun = Sun.from_builtin("E490_2014")
 
-        wave0 = sun.wave.to('um').value
+        wave0 = sun.wave.to("um").value
         fluxd0 = sun.fluxd.to(unit).value
 
         wave = np.linspace(0.35, 0.55, 6)
@@ -71,20 +79,21 @@ def test_sun_observe_wavelength_array(self):
         fluxd1 = np.zeros(len(wave))
         for i in range(len(wave)):
             j = (wave0 >= left_bins[i]) * (wave0 <= right_bins[i])
-            fluxd1[i] = (trapz(fluxd0[j] * wave0[j], wave0[j]) /
-                         trapz(wave0[j], wave0[j]))
+            fluxd1[i] = trapz(fluxd0[j] * wave0[j], wave0[j]) / trapz(
+                wave0[j], wave0[j]
+            )
 
         fluxd2 = sun.observe(wave * u.um, unit=unit).value
 
         assert np.allclose(fluxd1, fluxd2, rtol=0.005)
 
     def test_filt_units(self):
         """Colina et al. V=-26.75 mag, for zero-point flux density
-           36.7e-10 ergs/s/cm2/Å.
+        36.7e-10 ergs/s/cm2/Å.
         """
-        sun = Sun.from_builtin('E490_2014')
-        V = bandpass('johnson v')
-        weff, fluxd = sun.observe_bandpass(V, unit='erg/(s cm2 AA)')
+        sun = Sun.from_builtin("E490_2014")
+        V = bandpass("johnson v")
+        weff, fluxd = sun.observe_bandpass(V, unit="erg/(s cm2 AA)")
         assert np.isclose(weff.value, 5502, rtol=0.001)
         assert np.isclose(fluxd.value, 183.94, rtol=0.0003)
 
@@ -95,8 +104,8 @@ def test_filt_vegamag(self):
         agreement is good.
 
         """
-        sun = Sun.from_builtin('E490_2014')
-        V = bandpass('johnson v')
+        sun = Sun.from_builtin("E490_2014")
+        V = bandpass("johnson v")
         fluxd = sun.observe(V, unit=JMmag)
         assert np.isclose(fluxd.value, -26.75, atol=0.006)
 
@@ -107,8 +116,8 @@ def test_filt_abmag(self):
         optical.
 
         """
-        sun = Sun.from_builtin('E490_2014')
-        V = bandpass('johnson v')
+        sun = Sun.from_builtin("E490_2014")
+        V = bandpass("johnson v")
         fluxd = sun.observe(V, unit=u.ABmag)
         assert np.isclose(fluxd.value, -26.77, atol=0.007)
 
@@ -119,19 +128,19 @@ def test_filt_stmag(self):
         optical.
 
         """
-        sun = Sun.from_builtin('E490_2014')
-        V = bandpass('johnson v')
+        sun = Sun.from_builtin("E490_2014")
+        V = bandpass("johnson v")
         fluxd = sun.observe(V, unit=u.STmag)
         assert np.isclose(fluxd.value, -26.76, atol=0.003)
 
     def test_filt_solar_fluxd(self):
-        with solar_fluxd.set({'V': -26.76 * VEGAmag}):
+        with solar_fluxd.set({"V": -26.76 * VEGAmag}):
             sun = Sun(None)
-            fluxd = sun.observe('V', unit=VEGAmag)
+            fluxd = sun.observe("V", unit=VEGAmag)
         assert np.isclose(fluxd.value, -26.76)
 
     def test_meta(self):
-        sun = Sun.from_builtin('E490_2014')
+        sun = Sun.from_builtin("E490_2014")
         assert sun.meta is None
 
     @pytest.mark.remote_data
@@ -143,8 +152,8 @@ def test_kurucz_nan_error(self):
         NaNs in Kurucz file should not affect this calculation.
 
         """
-        sun = Sun.from_builtin('Kurucz1993')
-        V = bandpass('johnson v')
+        sun = Sun.from_builtin("Kurucz1993")
+        V = bandpass("johnson v")
         fluxd = sun.observe(V, unit=u.ABmag)
         assert np.isclose(fluxd.value, -26.77, atol=0.005)
 
@@ -163,15 +172,25 @@ def test_castelli96(self):
         2022-06-05: sbpy calculates 184.5 ergs/s/cm^2/A; agreement within 0.2%
         """
 
-        sun = Sun.from_builtin('Castelli1996')
-        V = bandpass('johnson v')
-        fluxd = sun.observe(V, unit='erg/(s cm2 AA)')
+        sun = Sun.from_builtin("Castelli1996")
+        V = bandpass("johnson v")
+        fluxd = sun.observe(V, unit="erg/(s cm2 AA)")
         assert np.isclose(fluxd.value, 184.2, rtol=0.002)
 
+    @pytest.mark.remote_data
+    def test_calspec(self):
+        """Verify CALSPEC solar model calibration."""
+
+        sun = Sun.from_builtin("calspec")
+        V = bandpass("johnson v")
+        fluxd = sun.observe(V, unit=VEGAmag)
+        assert np.isclose(fluxd.value, -26.75, rtol=0.002)
+
     def test_show_builtin(self, capsys):
         Sun.show_builtin()
         captured = capsys.readouterr()
         sources = inspect.getmembers(
-            Sun._sources, lambda v: isinstance(v, dict))
+            Sun._sources, lambda v: isinstance(v, dict)
+        )
         for k, v in sources:
             assert k in captured.out