dfm / python-fsps

Python bindings to Charlie Conroy's Flexible Stellar Population Synthesis (FSPS) Fortran code
https://python-fsps.readthedocs.io
MIT License
67 stars 38 forks source link

inconsistency in SEDs generated by tabular and functional SFHs #147

Closed mtakahiro closed 3 years ago

mtakahiro commented 3 years ago

Hi, I was checking the consistency between SEDs generated by a tabular SFH and functional one (here delayed tau model), somewhat similar to what was addressed previously (https://github.com/dfm/python-fsps/issues/91)

I found that resulting SEDs are different when metallicity was provided to the tabular SFHs, though the value was fixed to the solar, to match with the functional one (attached). When the metallicity was not provided, this does not happen.

The code I used was as following, and the version of fsps is 0.3.0. Could you advise me what makes the difference? When zcontinuous=0 for the functional SFH, the output SED changes but still differs from what I get from the tabular SFH. result

import fsps
import matplotlib.pyplot as pl

print(fsps.__version__)

# 1
sps = fsps.StellarPopulation(zcontinuous=1, logzsol=0)
# get the spectrum via the intrinsic SFH convolutions
sps.params["sfh"] = 4  #delayed_exponential
sdss_bands = fsps.find_filter("sdss")
mags = sps.get_mags(bands=sdss_bands)
age = (10 ** sps.log_age) / 1e9
sfr = sps.sfr
wavelength, spectrum = sps.get_spectrum(tage=13.7)

# 2
# get the spectrum from a tabular version
sps1 = fsps.StellarPopulation(zcontinuous=3)
sps1.params["sfh"] = 3
sps1.set_tabular_sfh(age, sfr, Z=sfr*0+1)
wavelength, spectrum1 = sps1.get_spectrum(tage=13.7)

# Plot the SFR
fig, axes = pl.subplots(2, 1)
axes[0].plot(age, sfr, '-o')
axes[0].set_xlabel('Age (Gyr)')
axes[0].set_ylabel('SFR (M$_\odot$/yr')
axes[1].plot(wavelength, spectrum, label="via sfh=4", lw=1)
axes[1].plot(wavelength, spectrum1, label="via sfh=3", linestyle='--', lw=0.5)
axes[1].set_xlim(1e3, 1e4)
axes[1].set_xlabel('$\lambda (\AA)$')
axes[1].set_ylabel('$f_\nu (L_\odot/Hz)$')
axes[1].legend(loc=0)
pl.show()
bd-j commented 3 years ago

When using set_tabular_sfh the units of Z are 'absolute', not Z/Z_sun or log(Z/Z_sun). So with sps.libraries == ('mist', 'miles') you would set Z=np.ones_like(sfr) * 0.014 (for padova Z_sun is 0.019). Perhaps the docstring for set_tabular_sfh could be worded more clearly. I've verified this produces spectra that are consistent within a few % (max difference around 2000AA) and the remaining difference likely due to the linear approximation used when SFH=3.

bd-j commented 3 years ago

also, when zcontinuous=0, logzsol is ignored and the metallicity is set by the integer zmet to one of the pre-computed SSP metallicity values (i.e., Z=sps.zlegend[sps.params["zmet"]]). I believe the default is zmet=1, which is the second lowest metallicity and therefore quite different from the default for zcontinuous=1. Sorry for the changing behavior with different zcontinuous values, it is inherited from fsps and we haven't made the effort to homogenize in the python interface.

mtakahiro commented 3 years ago

You are right! I was providing Z in the units of Z_sun. Now switched to absolute values, the resulting SEDs look identical. And thanks for the note on the case with zcontinuous=0. I was a bit confused with an unchanged outcome despite different metallicity inputs. Your comment here makes sense.

I will close this issue. Thanks a lot!