Access more FSPS quantities as attributes

src/fsps/fsps.f90

@@ -366,6 +366,17 @@ subroutine stellar_spectrum(ns,mact,logt,lbol,logg,phase,ffco,lmdot,wght,spec_ou
 
   end subroutine
 
+  subroutine get_ssp_weights(n_age, n_z, ssp_wghts_out)
+
+    ! Return the weights of each SSP in the CSP
+
+    implicit none
+    integer, intent(in) :: n_age,n_z
+    double precision, dimension(n_age,n_z), intent(inout) :: ssp_wghts_out
+    ssp_wghts_out = weight_ssp
+
+  end subroutine
+
   subroutine get_ssp_spec(ns,n_age,n_z,ssp_spec_out,ssp_mass_out,ssp_lbol_out)
 
     ! Return the contents of the ssp spectral array,
@@ -400,7 +411,7 @@ subroutine set_sfh_tab(ntab, age, sfr, met)
     sfh_tab(2, 1:ntabsfh) = sfr
     sfh_tab(3, 1:ntabsfh) = met
 
-  end subroutine set_sfh_tab
+  end subroutine
 
   subroutine set_ssp_lsf(nsv, sigma, wlo, whi)
 
@@ -414,7 +425,7 @@ subroutine set_ssp_lsf(nsv, sigma, wlo, whi)
     lsfinfo%maxlam = whi
     lsfinfo%lsf = sigma
 
-  end subroutine set_ssp_lsf
+  end subroutine
 
   subroutine get_setup_vars(cvms, vta_flag)
 
@@ -529,6 +540,16 @@ subroutine get_lambda(ns,lambda)
 
   end subroutine
 
+  subroutine get_res(ns,res)
+
+    ! Get the resolution array of the spectral library
+    implicit none
+    integer, intent(in) :: ns
+    double precision, dimension(ns), intent(out) :: res
+    res = spec_res
+
+  end subroutine
+
   subroutine get_libraries(isocname,specname,dustname)
 
     implicit none

src/fsps/fsps.py

@@ -534,6 +534,7 @@ def __init__(
         self._zcontinuous = zcontinuous
         # Caching.
         self._wavelengths = None
+        self._resolutions = None
         self._emwavelengths = None
         self._zlegend = None
         self._solar_metallicity = None
@@ -766,6 +767,21 @@ def _all_ssp_spec(self, update=True, peraa=False):
 
         return spec, mass, lbol
 
+    def _ssp_weights(self):
+         """Get the weights of the SSPs for the CSP
+
+         :returns weights:
+             The weights ``w`` of each SSP s.t. the total spectrum is the sum
+             :math:`L_{\lambda} = \sum_i,j w_i,j \, S_{i,j,\lambda}` where
+             math:`i,j` run over age and metallicity.
+         """
+
+         NTFULL = driver.get_ntfull()
+         NZ = driver.get_nz()
+         weights = np.zeros([NTFULL, NZ], order="F")
+         driver.get_ssp_weights(weights)
+         return weights
+
     def _get_stellar_spectrum(
         self,
         mact,
@@ -1032,6 +1048,16 @@ def wavelengths(self):
             self._wavelengths = driver.get_lambda(NSPEC)
         return self._wavelengths.copy()
 
+    @property
+    def resolutions(self):
+        r"""The resolution array, in km/s dispersion. Negative numbers indicate
+         poorly defined, approximate, resolution (based on coarse opacity
+         binning in theoretical spectra)"""
+        if self._resolutions is None:
+            NSPEC = driver.get_nspec()
+            self._resolutions = driver.get_res(NSPEC)
+        return self._resolutions.copy()
+
     @property
     def emline_wavelengths(self):
         r"""Emission line wavelengths, in Angstroms"""

tests/tests.py

@@ -20,56 +20,22 @@ def _reset_default_params(pop, params):
 
 
 def test_isochrones(pop_and_params):
-    # Just test that `isochrones()` method runs
-    pop, params = pop_and_params
-    _reset_default_params(pop, params)
-    pop.params["imf_type"] = 0
-    iso = pop.isochrones()
+    """Just test that `isochrones()` method runs"""
 
+    # recomputes SSPs
 
-def test_tabular(pop_and_params):
     pop, params = pop_and_params
     _reset_default_params(pop, params)
+    pop.params["imf_type"] = 0
+    iso = pop.isochrones()
 
-    import os
-
-    fn = os.path.join(os.environ["SPS_HOME"], "data/sfh.dat")
-    age, sfr, z = np.genfromtxt(fn, unpack=True, skip_header=0)
-
-    # Mono-metallicity
-    pop.params["sfh"] = 3
-    pop.set_tabular_sfh(age, sfr)
-    w, spec = pop.get_spectrum(tage=0)
-    pop.set_tabular_sfh(age, sfr)
-    assert pop.params.dirty
-    w, spec = pop.get_spectrum(tage=0)
-    assert spec.shape[0] == len(pop.ssp_ages)
-    assert pop.params["sfh"] == 3
-    w, spec_last = pop.get_spectrum(tage=-99)
-    assert spec_last.ndim == 1
-    w, spec = pop.get_spectrum(tage=age.max())
-    assert np.allclose(spec / spec_last - 1.0, 0.0)
-    pop.params["logzsol"] = -1
-    w, spec_lowz = pop.get_spectrum(tage=age.max())
-    assert not np.allclose(spec / spec_lowz - 1.0, 0.0)
-
-    # Multi-metallicity
-    pop._zcontinuous = 3
-    pop.set_tabular_sfh(age, sfr, z)
-    w, spec_multiz = pop.get_spectrum(tage=age.max())
-    assert not np.allclose(spec_lowz / spec_multiz - 1.0, 0.0)
 
-    pop._zcontinuous = 1
-    pop.set_tabular_sfh(age, sfr)
-    # get mass weighted metallicity
-    mbin = np.gradient(age) * sfr
-    mwz = (z * mbin).sum() / mbin.sum()
-    pop.params["logzsol"] = np.log10(mwz / 0.019)
-    w, spec_onez = pop.get_spectrum(tage=age.max())
-    assert not np.allclose(spec_onez / spec_multiz - 1.0, 0.0)
+def test_imf3(pop_and_params):
+    """Make sure that changing the (upper) imf changes the parameter dirtiness
+    and also the SSP spectrum"""
 
+    # recomputes SSPs
 
-def test_imf3(pop_and_params):
     pop, params = pop_and_params
     _reset_default_params(pop, params)
     pop.params["imf_type"] = 2
@@ -92,6 +58,9 @@ def test_imf3(pop_and_params):
 
 
 def test_param_checks(pop_and_params):
+
+    # recomputes SSPs
+
     pop, params = pop_and_params
     _reset_default_params(pop, params)
     pop.params["sfh"] = 1
@@ -113,6 +82,9 @@ def test_param_checks(pop_and_params):
 
 
 def test_smooth_lsf(pop_and_params):
+
+    # recomputes SSPs
+
     pop, params = pop_and_params
     _reset_default_params(pop, params)
     tmax = 1.0
@@ -133,8 +105,62 @@ def test_smooth_lsf(pop_and_params):
     assert pop.params.dirtiness == 2
 
 
+def test_tabular(pop_and_params):
+    """Test that you get the right shape spectral arrays for tabular SFHs, that
+    the parameter dirtiness is appropriately managed for changing tabular SFH,
+    and that multi-metallicity SFH work."""
+
+    # uses default SSPs, but makes them for every metallicity
+
+    pop, params = pop_and_params
+    _reset_default_params(pop, params)
+
+    import os
+
+    fn = os.path.join(os.environ["SPS_HOME"], "data/sfh.dat")
+    age, sfr, z = np.genfromtxt(fn, unpack=True, skip_header=0)
+
+    # Mono-metallicity
+    pop.params["sfh"] = 3
+    pop.set_tabular_sfh(age, sfr)
+    w, spec = pop.get_spectrum(tage=0)
+    pop.set_tabular_sfh(age, sfr)
+    assert pop.params.dirty
+    w, spec = pop.get_spectrum(tage=0)
+    assert spec.shape[0] == len(pop.ssp_ages)
+    assert pop.params["sfh"] == 3
+    w, spec_last = pop.get_spectrum(tage=-99)
+    assert spec_last.ndim == 1
+    w, spec = pop.get_spectrum(tage=age.max())
+    assert np.allclose(spec / spec_last - 1.0, 0.0)
+    pop.params["logzsol"] = -1
+    w, spec_lowz = pop.get_spectrum(tage=age.max())
+    assert not np.allclose(spec / spec_lowz - 1.0, 0.0)
+
+    # test the formed mass for single age
+    assert np.allclose(np.trapz(sfr, age) * 1e9, pop.formed_mass)
+
+    # Multi-metallicity
+    pop._zcontinuous = 3
+    pop.set_tabular_sfh(age, sfr, z)
+    w, spec_multiz = pop.get_spectrum(tage=age.max())
+    assert not np.allclose(spec_lowz / spec_multiz - 1.0, 0.0)
+
+    pop._zcontinuous = 1
+    pop.set_tabular_sfh(age, sfr)
+    # get mass weighted metallicity
+    mbin = np.gradient(age) * sfr
+    mwz = (z * mbin).sum() / mbin.sum()
+    pop.params["logzsol"] = np.log10(mwz / pop.solar_metallicity)
+    w, spec_onez = pop.get_spectrum(tage=age.max())
+    assert not np.allclose(spec_onez / spec_multiz - 1.0, 0.0)
+
+
 def test_get_mags(pop_and_params):
     """Basic test for supplying filter names to get_mags"""
+
+    # uses default SSPs
+
     pop, params = pop_and_params
     _reset_default_params(pop, params)
     fuv1 = pop.get_mags(bands=["galex_fuv"])[:, 0]
@@ -146,6 +172,11 @@ def test_get_mags(pop_and_params):
 
 
 def test_ssp(pop_and_params):
+    """Test that you get a sensible wavelength array, and that you get a
+    sensible V-band magnitude for a 1-Gyr SSP"""
+
+    # uses default SSPs
+
     pop, params = pop_and_params
     _reset_default_params(pop, params)
     pop.params["sfh"] = 0
@@ -162,6 +193,9 @@ def test_ssp(pop_and_params):
 
 def test_libraries(pop_and_params):
     """This does not require or build clean SSPs"""
+
+    # uses default SSPs
+
     pop, params = pop_and_params
     _reset_default_params(pop, params)
     ilib, splib, dlib = pop.libraries
@@ -172,6 +206,9 @@ def test_libraries(pop_and_params):
 
 def test_filters():
     """Test all the filters got transmission data loaded."""
+
+    # uses default SSPs
+
     flist = filters.list_filters()
     # force trans cache to load
     filters.FILTERS[flist[0]]._load_transmission_cache()
@@ -180,6 +217,9 @@ def test_filters():
 
 
 def test_csp_dirtiness(pop_and_params):
+    """Make sure that changing CSP parameters increases dirtiness to 1"""
+    # uses default SSPs
+
     pop, params = pop_and_params
     _reset_default_params(pop, params)
     pop.params["sfh"] = 1
@@ -195,6 +235,9 @@ def test_redshift(pop_and_params):
     1. redshifting does not persist in cached arrays
     2. specifying redshift via get_mags keyword or param key are consistent
     """
+
+    # uses default SSPs
+
     pop, params = pop_and_params
     _reset_default_params(pop, params)
     pop.params["sfh"] = 0
@@ -216,6 +259,9 @@ def test_redshift(pop_and_params):
 
 def test_nebemlineinspec(pop_and_params):
     """Make sure nebular lines are actually added."""
+
+    # uses default SSPs
+
     pop, params = pop_and_params
     _reset_default_params(pop, params)
     pop.params["sfh"] = 4
@@ -233,6 +279,9 @@ def test_nebemlineinspec(pop_and_params):
 
 
 def test_mformed(pop_and_params):
+    """Make sure formed mass integrates to 1 for parameteric SFH"""
+    # uses default SSPs
+
     pop, params = pop_and_params
     _reset_default_params(pop, params)
     pop.params["sfh"] = 1
@@ -246,6 +295,11 @@ def test_mformed(pop_and_params):
 
 
 def test_light_ages(pop_and_params):
+    """Make sure light weighting works, and gives sensible answers for the
+    light-weighted age in the FUV and mass-weighted age stored in
+    `stellar_mass`"""
+    # uses default SSPs
+
     pop, params = pop_and_params
     _reset_default_params(pop, params)
     tmax = 5.0
@@ -270,13 +324,41 @@ def test_light_ages(pop_and_params):
 
 def test_smoothspec(pop_and_params):
     # FIXME: This is not very stringent
+
+    # uses default SSPs
+
     pop, params = pop_and_params
     _reset_default_params(pop, params)
     wave, spec = pop.get_spectrum(tage=1, peraa=True)
     spec2 = pop.smoothspec(wave, spec, 160.0, minw=1e3, maxw=1e4)
     assert (spec - spec2 == 0.0).sum() > 0
 
 
+def test_ssp_weights(pop_and_params):
+    """Check that weights dotted into ssp is the same as the returned spectrum
+    when there's no dust or emission lines and zcontinuous=0"""
+
+    # uses default SSPs
+
+    pop, params = pop_and_params
+    _reset_default_params(pop, params)
+
+    import os
+    fn = os.path.join(os.environ["SPS_HOME"], "data/sfh.dat")
+    age, sfr, z = np.genfromtxt(fn, unpack=True, skip_header=0)
+    pop.params["sfh"] = 3
+    pop.set_tabular_sfh(age, sfr)
+    zind = -3
+    pop.params["logzsol"] = np.log10(pop.zlegend[zind]/pop.solar_metallicity)
+
+    wave, spec = pop.get_spectrum(tage=age.max())
+    mstar = pop.stellar_mass
+    wght = pop._ssp_weights()
+    ssp, smass, slbol = pop._all_ssp_spec()
+
+    assert np.allclose((smass[:, zind] * wght[:, 0]).sum(), mstar)
+
+
 # Requires scipy
 # def test_sfr_avg():