Metallicity Histories in Tabular SFH

[bd-j]

This PR adds support for metallicity histories in tabular SFHs through the zcontinuous=3 flag and the set_tabular_sfh() method.

Similar to the existing tabular SFH options (see #63) usage is

import numpy as np
import fsps
sps = fsps.StellarPopulation(zcontinuous=3)
sps.params['sfh'] = 3
age = np.linspace(1, 10, 50)
sfr = np.exp(-age/3.0)
z = age / age.max() * 0.019
sps.set_tabular_sfh(age, sfr, z)
wave, spec = sps.get_spectrum(tage=0)

Note that the zcontinuous=3 option must only be used for sfh=3, Otherwise an unrecoverable fortran error will occur. Furthermore nebular emission cannot be added when multi-metallicity tabular SFH is computed.

This PR also adds a script to the demo/ directory that shows how to successively turn on and use some of the StellarPopulation options. It also addresses a documentation typo (#82)

[bd-j]

For fun, here is the spectrum from the SFH in the default sfh.dat that ships with FSPS, computed with a single mass-weighted metallicity and including the full metallicity history. Produced using the code snippet below.

import os
import numpy as np
import matplotlib.pyplot as pl
import fsps

fn = os.path.join(os.environ['SPS_HOME'], 'data/sfh.dat')
age, sfr, z = np.genfromtxt(fn, unpack=True, skip_header=0)

pop = fsps.StellarPopulation(zcontinuous=1)
pop.params['sfh'] = 3

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=-99)

pop._zcontinuous = 3
pop.set_tabular_sfh(age, sfr, z)
w, spec_multiz = pop.get_spectrum(tage=-99)

fig, axes = pl.subplots(2, 1, sharex=True)
ax = axes[0]
ax.plot(w, spec_onez, label="$Z_{{mass-weighted}}$")
ax.plot(w, spec_multiz, label="$Z_{{history}}$")
ax.legend()
ax.set_ylabel('$f_\\nu$')
ax = axes[1]
ax.plot(w, (spec_onez/spec_multiz - 1) * 100)
ax.set_ylabel("% difference")
ax.set_ylim(-10, 40)
ax.set_xlabel('wavelength')
[a.set_xlim(1.1e3, 0.95e4) for a in axes.flat]

[jrleja]

I’m surprised by the large delta in the UV. I believe it would be smaller if instead the luminosity-weighted metallicities were set equal. I think most of the delta is because the youngest stars have a considerably lower metallicity in the “onez” model.

For comparison, here's a Z(t) versus constant-Z comparison I did, with a constant SFH. Z(t) was set by a closed-box evolution model with the following values:

0-100 Myr, Z = 1 Zsun
100-300 Myr, Z = 0.96 Zsun
300 Myr-1 Gyr, Z = 0.87 Zsun
1 Gyr - 3 Gyr, Z = 0.66 Zsun
3-6 Gyr, Z = 0.41 Zsun
6-13.6 Gyr, Z = 0.13 Zsun

while the constant-Z value was taken to be Zsun. Emission lines are included because I did this the slow way, by generating 6 different populations and summing them.

[bd-j]

@jrleja I agree! Using the luminosity weighted metallicity would be a much fairer comparison and I expect it would show a much smaller UV difference.