diff --git a/src/fsps/fsps.f90 b/src/fsps/fsps.f90 index f43c72e8..7eff93be 100644 --- a/src/fsps/fsps.f90 +++ b/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 diff --git a/src/fsps/fsps.py b/src/fsps/fsps.py index e4856bfc..a5729772 100644 --- a/src/fsps/fsps.py +++ b/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""" diff --git a/tests/tests.py b/tests/tests.py index 014354f3..4df7dbe2 100644 --- a/tests/tests.py +++ b/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,6 +324,9 @@ 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) @@ -277,6 +334,31 @@ def test_smoothspec(pop_and_params): 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():