bc03.py 6.05 KB
Newer Older
1
# -*- coding: utf-8 -*-
2 3
# Copyright (C) 2013 Centre de données Astrophysiques de Marseille
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
Yannick Roehlly's avatar
Yannick Roehlly committed
4
# Author: Yannick Roehlly
5 6 7 8

import numpy as np


9
class BC03(object):
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
    """Single Stellar Population as defined in Bruzual and Charlot (2003)

    This class holds the data associated with a single stellar population
    (SSP) as defined in Bruzual and Charlot (2003). Compare to the pristine
    SSP, the wavelength are given in nm (rather than Å), the time is given in

    """

    def __init__(self, imf, metallicity, time_grid, wavelength_grid,
                 color_table, lumin_table):
        """Create a new single stellar population as defined in Bruzual and
        Charlot (2003).

        Parameters
        ----------
25
        imf: string
26 27
            Initial mass function (IMF): either 'salp' for Salpeter (1955) or
            'chab' for Chabrier (2003).
28
        metallicity: float
29 30
            The metallicity. Possible values are 0.0001, 0.0004, 0.004, 0.008,
            0.02 (solar metallicity) and 0.05.
31
        time_grid: array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
32
            The time [Myr] grid used in the colors_table and the lumin_table.
33
        wavelength_grid: array of floats
34 35 36 37 38 39 40
            The wavelength [nm] grid used in the lumin_table.
        color_table: 2 axis array of floats
            Array containing information from some of the *.?color tables from
            Bruzual and Charlot (2003) at each time of the time_grid.
                * color_table[0]: Total mass in stars in solar mass
                * color_table[1]: Mass returned to the ISM by evolved stars
                                   in solar mass
41
                * color_table[2]: rate of H-ionizing photons (s-1)
42 43 44 45 46 47
                * color_table[3]: Amplitude of 4000 Å break (Bruzual 2003)
                * color_table[4]: Amplitude of 4000 Å narrow break (Balogh
                                   et al. 1999)
                * color_table[5]: Amplitude of 4000 Å break (Stoughton
                                   et al. 2002)
                * color_table[6]: Amplitude of Lyman discontinuity
48
        lumin_table: 2 axis array of floats
49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
            Luminosity density in W/nm. The first axis is the wavelength and
            the second the time (index bases on the wavelength and time grids).
        """

        if imf in ['salp', 'chab']:
            self.imf = imf
        else:
            raise ValueError('IMF must be either sal for Salpeter or '
                             'cha for Chabrier.')
        self.metallicity = metallicity
        self.time_grid = time_grid
        self.wavelength_grid = wavelength_grid
        self.color_table = color_table
        self.lumin_table = lumin_table

    def convolve(self, sfh_time, sfh_sfr, norm=False):
        """Convolve the SSP with a Star Formation History

        Given a SFH (an time grid and the corresponding star formation rate
        SFR), this method convolves the color table and the SSP luminosity
        spectrum along the whole SFR.

        The time grid of the SFH is expected to be ordered and must not run
        beyong 20 Gyr (the maximum time for Bruzual and Charlot 2003 SSP).

        Parameters
        ----------
76
        sfh_time: array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
77 78 79
            Time grid [Myr] of the star formation history. It must be
            increasing and not run beyond 20 Gyr. The SFH will be regrided to
            the SSP time.
80 81 82 83 84 85 86
        sfh_sfr: array of floats
            Star Formation Rates in Msun/yr at each time of the SFH time grid.
        norm: boolean
            If true, the sfh will be normalised to 1 solar mass produced.

        Returns
        -------
87
        wavelength: array of floats
88
            Wavelength grid [nm] for the spectrum
89
        luminosity: array of floats
90 91 92 93 94 95
            Luminosity density [W/nm] at each wavelength.
        bc03_info: dictionary
            Dictionary containing various information from the *.?color tables:
            - "m_star": Total mass in stars in solar mass
            - "m_gas": Mass returned to the ISM by evolved stars
                       in solar mass
96
            - "n_ly": rate of H-ionizing photons (s-1)
97 98
            - "b_4000": Amplitude of 4000 Å break (Bruzual 2003)
            - "b4_vn": Amplitude of 4000 Å narrow break (Balogh et al. 1999)
99
            - "b4_sdss": Amplitude of 4000 Å break (Stoughton et al. 2002)
100
            - "b_912": Amplitude of Lyman discontinuity
Yannick Roehlly's avatar
Yannick Roehlly committed
101

102 103 104 105 106 107 108 109 110 111 112 113 114 115 116
        """
        # We work on a copy of SFH (as we change it)
        sfh_time, sfh_sfr = np.copy((sfh_time, sfh_sfr))

        # Index, in the SSP time grid, of the time nearest to the age of
        # the SFH.
        idx = np.abs(self.time_grid - np.max(sfh_time)).argmin()

        # We regrid the SFH to the time grid of the SSP using a linear
        # interpolation. If the SFH does no start at 0, the first SFR values
        # will be set to 0.
        sfh_sfr = np.interp(self.time_grid[:idx + 1],
                            sfh_time, sfh_sfr,
                            left=0., right=0.)

Yannick Roehlly's avatar
Yannick Roehlly committed
117
        # Step between two item in the time grid in Myr
118 119 120 121
        step = self.time_grid[1] - self.time_grid[0]

        # If needed, we normalise the SFH to 1 solar mass produced.
        if norm:
Yannick Roehlly's avatar
Yannick Roehlly committed
122
            sfh_sfr = sfh_sfr / np.trapz(sfh_sfr * 1.e6,
123 124 125 126 127 128
                                         self.time_grid[:idx + 1])

        # As both the SFH and the SSP (limited to the age of the SFH) data now
        # share the same time grid, the convolution is just a matter of
        # reverting one and computing the sum of the one to one product; this
        # is done using the dot product.
129 130
        color_table = self.color_table[:,:idx + 1]
        lumin_table = self.lumin_table[:,:idx + 1]
131

Yannick Roehlly's avatar
Yannick Roehlly committed
132 133 134
        # The 1.e6 * step is because the SFH is in solar mass per year.
        color_info = 1.e6 * step * np.dot(color_table, sfh_sfr[::-1])
        luminosity = 1.e6 * step * np.dot(lumin_table, sfh_sfr[::-1])
135 136 137 138 139 140 141

        bc03_info = dict(zip(
            ["m_star", "m_gas", "n_ly", "b_4000", "b4_vn", "b4_sdss", "b_912"],
            color_info
        ))

        return self.wavelength_grid, luminosity, bc03_info