How to use the fsps.StellarPopulation function in fsps

# -*- 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