How to use fsps - 10 common examples

# -*- coding: utf-8 -*-

from __future__ import (division, print_function, absolute_import,
                        unicode_literals)
import pyfits
import numpy as np
import matplotlib.pyplot as pl
import fsps

# Measurements of cluster parameters.
tage = 10. ** (8.04 - 9)
logmass = 4.09
dist_mod = 24.5

# Set up the stellar population model.
sp = fsps.StellarPopulation(imf_type=2, dust_type=1, mwr=3.1, dust2=0.3)

# The measured magnitudes from the literature.
data = {"wfc3_f160w": 16.386,
        "wfc3_f275w": 17.398,
        "wfc_acs_f814w": 17.155}

# There are also a few filters that we have data for but aren't included in
# the standard FSPS install:
#   "F110W": 16.877,
#   "F336W": 17.349,
#   "F475W": 17.762,

# Load the observed spectrum.
f = pyfits.open("DAO69.fits")
obs_spec = np.array(f[0].data, dtype=float)
f.close()

:ages: 
        The times at which the sfr is defined for any number of sfhs, H. Shape (H, SFR). 
        
        :time_steps:
        An array of time values for which a spectrum with the given sfr is returned. Shape (T,)

        :zmets:
        An array of model metallicity values (range 1-22, see FSPS details) for which a spectrum with the given sfr is returned. Shape (Z,)
        
    RETURNS:
        :fsps_flux:
        Fluxes at each manga wavelength, M for the number of sfhs input at each time_step and zmet combination. Shape (H*T*Z, M). 

    """

    sp = fsps.StellarPopulation(compute_vega_mags=False, zcontinuous=0, zmet=18, sfh=3, dust_type=2, dust2=0.2, 
        add_neb_emission=True, add_neb_continuum=False, imf_type=1, add_dust_emission=True, min_wave_smooth=3600, max_wave_smooth=7500, smooth_velocity=True, sigma_smooth=77.)

    sp.set_tabular_sfh(age = ages, sfr=sfr)

    def time_spec(params):
        sps = sp.get_spectrum(tage=params[0], zmet=params[1])[1]
        return sps

    fsps_wave = sp.get_spectrum()[0]
    fsps_spec = np.array(list(map(time_spec, list(product(time_steps, zmets)))))

    manga_wave = np.load("manga_wavelengths_AA.npy")

    f = interpolate.interp1d(fsps_wave, fsps_spec)
    fluxes = f(manga_wave)

def newstars_gen(stars_list):
    global sp
    if sp is None:
        sp = fsps.StellarPopulation()

    #the newstars (particle type 4; so, for cosmological runs, this is all
    #stars) are calculated in a separate function with just one argument so that it is can be fed 
    #into pool.map for multithreading.
    #sp = fsps.StellarPopulation()
    sp.params["tage"] = stars_list[0].age
    sp.params["imf_type"] = cfg.par.imf_type
    sp.params["pagb"] = cfg.par.pagb
    sp.params["sfh"] = 0
    sp.params["zmet"] = stars_list[0].fsps_zmet
    sp.params["add_neb_emission"] = False
    sp.params["add_agb_dust_model"] = cfg.par.add_agb_dust_model
    sp.params['gas_logu'] = cfg.par.gas_logu

    if cfg.par.FORCE_gas_logz == False:
        sp.params['gas_logz'] = np.log10(stars_list[0].metals/cfg.par.solar)

if label is not None:
        ax.text(0.63, 0.85, label, transform=ax.transAxes, fontsize=16)
    
    return fig, ax

if __name__ == "__main__":

    pl.rc('text', usetex=True)
    pl.rc('font', family='serif')
    pl.rc('axes', grid=False)
    pl.rc('xtick', direction='in')
    pl.rc('ytick', direction='in')
    pl.rc('xtick', top=True)
    pl.rc('ytick', right=True)

    sps = fsps.StellarPopulation(zcontinuous=1)
    pdf = PdfPages('features.pdf')

    # Basic spectrum
    sps.params['sfh'] = 4
    sps.params['tau'] = 5.0
    sps.params['logzsol'] = 0.0
    sps.params['dust_type'] = 4 # kriek and Conroy
    sps.params['imf_type'] = 2  #kroupa
    sps.params['imf3'] = 2.3
    fig, ax, spec = makefig(sps)
    fig, ax = prettify(fig, ax, label="$\\tau=5$, Age$=13.7$,\n$\log Z/Z_\odot=0.0$")
    pdf.savefig(fig)
    pl.close(fig)
    
    # change IMF
    sps.params['imf3'] = 2.5

def calc_emline(stars_list):
    print ('[SED_gen/calc_emline]: Calculating Emission Line Fluxes')
    #this function is awful and redundant.  itloops over just the new
    #stars in the box and calculates the SEDs for the ones smaller
    #than cfg.par.HII_max_age.  its a bit wasteful since we already
    #calculate these SEDs, but its difficult to get the file to save
    #in a nice format without incl
    global sp
    if sp is None:
        sp = fsps.StellarPopulation()

    #set up arrays 
    #first how many young stars are there?
    
    newstars_idx = []
    for counter,star in enumerate(stars_list):
        if star.age <= cfg.par.HII_max_age:
            newstars_idx.append(counter)
    num_newstars = len(newstars_idx)
 
    #set up a dummy sps model just to get number of wavelengths
    sp.params["tage"] = stars_list[0].age
    sp.params["zmet"] = stars_list[0].fsps_zmet
    sp.params["add_neb_emission"] = True
    wav,spec = sp.get_spectrum()
    n_emlines = len(sp.emline_wavelengths)

# This demo shows how to make a plot of the fractional contribution of differnt
# stellar mass ranges to the total spectrum, for different properties of the
# stellar population.  There is also some use of the filter objects.

from itertools import product
import numpy as np
import matplotlib.pyplot as pl
import fsps

# make an SPS object with interpolation in metallicity enabled, and
# set a few parameters
sps = fsps.StellarPopulation(zcontinuous=1)
sps.params['imf_type'] = 0   # Salpeter IMF
sps.params['sfh'] = 1  # Tau model SFH

# Get a few filter transmission curves
filterlist = {'galex_fuv': 'FUV', 'galex_nuv': 'NUV', 'sdss_u': 'u',
              'sdss_g': 'g', 'sdss_i': 'i', '2mass_j': 'J'}
filters = [fsps.filters.get_filter(f) for f in filterlist]


# This is the main function
def specplots(tage=10.0, const=1.0, tau=1.0, neb=False, z=0.0, savefig=True,
              masslims=[1.0, 3.0, 15.0, 30.0, 120.0], **kwargs):
    """Make a number of plots showing the fraction of the total light due to
    stars within certain stellar mass ranges.

    :param tage:

An array of time values for which a spectrum with the given sfr is returned. Shape (T,)

        :zmets:
        An array of model metallicity values (range 1-22, see FSPS details) for which a spectrum with the given sfr is returned. Shape (Z,)
        
    RETURNS:
        :fsps_flux:
        Fluxes at each manga wavelength, M for the number of sfhs input at each time_step and zmet combination. Shape (H*T*Z, M). 

    """
    time_step = [params[0]]
    zmets = params[1]
    ages = params[2]
    sfr = params[3]

    sp = fsps.StellarPopulation(compute_vega_mags=False, zcontinuous=0, zmet=18, sfh=3, dust_type=2, dust2=0.2, 
        add_neb_emission=True, add_neb_continuum=False, imf_type=1, add_dust_emission=True, min_wave_smooth=3600, max_wave_smooth=7500, smooth_velocity=True, sigma_smooth=77.)

    sp.set_tabular_sfh(age = ages, sfr=sfr)

    def time_spec(params):
        sps = sp.get_spectrum(tage=params[0], zmet=params[1])[1]
        return sps

    fsps_wave = sp.get_spectrum()[0]
    fsps_spec = np.array(list(map(time_spec, list(product(time_step, zmets)))))

    manga_wave = np.load("manga_wavelengths_AA.npy")
    manga_hz = c*(un.km/un.s).to(un.AA/un.s)/manga_wave

    f = interpolate.interp1d(fsps_wave, fsps_spec)
    fluxes = f(manga_wave)*(manga_hz/manga_wave) # Fluxes returned by FSPS are in units L_sol/Hz - MaNGA DAP needs fluxes in per AA flux units

def _get_grid_stats(self):
        if self.params.dirty:
            self._compute_csp()

        if self._stats is None:
            self._stats = driver.get_stats(driver.get_ntfull(),driver.get_nemline())

        return self._stats

fagn=0.0,
            agn_tau=10.0
        )

        # Parse any input options.
        for k, v in self.params.iteritems():
            self.params[k] = kwargs.pop(k, v)

        # Make sure that we didn't get any unknown options.
        if len(kwargs):
            raise TypeError("__init__() got an unexpected keyword argument "
                            "'{0}'".format(list(kwargs)[0]))

        # Before the first time we interact with the FSPS driver, we need to
        # run the ``setup`` method.
        if not driver.is_setup:
            driver.setup(compute_vega_mags, vactoair_flag)
        else:
            cvms, vtaflag = driver.get_setup_vars()
            assert compute_vega_mags == bool(cvms)
            assert vactoair_flag == bool(vtaflag)
        self._zcontinuous = zcontinuous
        # Caching.
        self._wavelengths = None
        self._emwavelengths = None
        self._zlegend = None
        self._ssp_ages = None
        self._stats = None
        self._libraries = None

def check_params(self):
        NZ = driver.get_nz()
        assert self._params["zmet"] in range(1, NZ + 1), \
            "zmet={0} out of range [1, {1}]".format(self._params["zmet"], NZ)
        assert self._params["dust_type"] in range(5), \
            "dust_type={0} out of range [0, 4]".format(
                self._params["dust_type"])
        assert self._params["imf_type"] in range(6), \
            "imf_type={0} out of range [0, 5]".format(self._params["imf_type"])
        assert (self._params["tage"] <= 0) | (self._params["tage"] > self._params["sf_start"]), \
            "sf_start={0} is greater than tage={1}".format(
                self._params["sf_start"], self._params["tage"])
        assert (self._params["const"]+self._params["fburst"]) <= 1, \
            "const + fburst > 1"