Commit 832d27ae authored by Yannick Roehlly's avatar Yannick Roehlly
Browse files

Correct the Maraston 2005 module to use SFH

The Marasont 2005 module is modified to use the output of a star
formation history module. The SFH is now a complete SFR vs time tuple
(taken from the SED information dictionary) and the age of the galaxy is
defined by this SFH (the maximum time).
The Maraston 2005 data class is also changed to make the convolution in
such SFH.
parent 2c0d855c
# -*- coding: utf-8 -*-
"""
Copyright (C) 2012 Centre de données Astrophysiques de Marseille
Copyright (C) 2012, 2013 Centre de données Astrophysiques de Marseille
Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
@author: Yannick Roehlly <yannick.roehlly@oamp.fr>
......@@ -8,7 +8,6 @@ Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
"""
import numpy as np
from scipy import interpolate
class SspM2005(object):
......@@ -18,7 +17,7 @@ class SspM2005(object):
(SSP) as defined in Maraston (2005). Compare to the pristine Maraston's
SSP:
- The age grid used ranges from 1 My to 13.7 Gyr with 1 My step. This
- The time grid used ranges from 1 My to 13.7 Gyr with 1 My step. This
excludes the metallicities -2.25 and 0.67 for which the pristine age
grid starts at 1 Gyr.
......@@ -49,7 +48,7 @@ class SspM2005(object):
* -0.33 (corresponding to 0.5 Zsun)
* -1.35 (corresponding to 1/50 Zsun)
time_grid : array of floats
The age [Gyr] grid used in the mass_table and the spec_table.
The time [Gyr] grid used in the mass_table and the spec_table.
wavelength_grid : array of floats
The wavelength [nm] grid used in the spec_table.
mass_table : (6, n) array of floats
......@@ -63,7 +62,7 @@ class SspM2005(object):
* mass_table[4]: black hole star mass
* mass_table[5]: mass in the turn off
spec_table : (2, n) array of floats
The 2D table giving the luminosity density [W/nm] at various age.
The 2D table giving the luminosity density [W/nm] at various time.
The first axis is the age, base on the time_grid, the second is the
wavelength, base on the wavelength_grid.
......@@ -80,28 +79,27 @@ class SspM2005(object):
self.mass_table = mass_table
self.spec_table = spec_table
def convolve(self, sfh, age, norm=False):
"""
Given a Star Formation History (SFH) and an age, this method convolves
the mass distribution and the spectrum at a given age.
def convolve(self, sfh_time, sfh_sfr, norm=False):
"""Convolve the SSP with a Star Formation History
If 'age' is an array, the full convolution is computed and the method
returns the linear interpolation of the convolution result at each
age. This can take a few minutes.
Given a SFH (an time grid and the corresponding star formation rate
SFR), this method convolves the mass distribution and the SSP spectrum
along the whole SFR.
If 'age' is a float, the computation is a lot faster, based on the
fact that the SFR is given at each age of the SSP age grid.
The time grid of the SFH is expected to be ordered and must not run
beyong 13.7 Gyr (the maximum time for Maraston 2005 SSP).
Parameters
----------
sfh : array of floats
Star Formation Rates in Msun/y for each age of the SSP age grid.
age : float or array of floats
Age(s) in Gyr we want the mass distribution at. If an array is
given, the various outputs will be arrays of the same length.
sfh_time : array of floats
Time grid of the star formation history. It must be increasing and
not run beyond 13.7 Gyr. As the SFH will be regrided to the SSP
time grid, it is better to have a SFH age grid compatible, i.e.
with a precision limited to 1e-3 Gyr.
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 at
'age' (or max(age) if it's an array).
If true, the sfh will be normalised to 1 solar mass produced.
Returns
-------
......@@ -120,83 +118,36 @@ class SspM2005(object):
"""
# We work on a copy of SFH (as we change it)
sfh = np.copy(sfh)
# The SFH must be on the same age grid as the SSP.
if not len(sfh) == len(self.time_grid):
raise ValueError("The star formation rate must be base"
"on the same grid than the SSP.")
# Check if the age parameter is a unique float or an array.
try:
age = float(age)
isAgeUnique = True
except:
age = np.array(age, dtype='float')
isAgeUnique = False
# Step between two item in the age grid in Gyr
step = self.time_grid[1] - self.time_grid[0]
sfh_time, sfh_sfr = np.copy((sfh_time, sfh_sfr))
if isAgeUnique:
# This is the fast convolution technique
# 1. We limit sfh, mass_table and spec_table for time up to
# the age.
# Index of the time nearest to age in the grid.
idx = np.abs(self.time_grid - age).argmin()
sfh = sfh[:idx + 1]
mass_table = self.mass_table[:, :idx + 1]
spec_table = self.spec_table[:, :idx + 1]
# 2. If needed, we normalise the SFH to 1 solar mass formed at
# 'age'.
if norm:
sfh = sfh / (step * sfh.sum())
# 3. We must convolve the mass evolution array with the SFH.
# As both tables share the same time grid, it's just a matter
# of reverting one and computing the sum of the one to one
# product. This is done using the dot product.
# The 1.e9 * step is because the SFH is in solar mass per year.
masses = 1.e9 * step * np.dot(mass_table, sfh[::-1])
#4. We do the same thing for the spectre.
spectra = 1.e9 * step * np.dot(spec_table, sfh[::-1])
# 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()
else:
# If the age parameter is an array, we do the full convolution.
# 1. We limit sfh, mass_table, spec_table and time_grid for time
# up to the maximum age.
# Index of the time nearest to the maximum age in the grid.
idx = np.abs(self.time_grid - np.max(age)).argmin()
sfh = sfh[:idx + 1]
mass_table = self.mass_table[:, :idx + 1]
spec_table = self.spec_table[:, :idx + 1]
time_grid = self.time_grid[:, :idx + 1]
# If needed, we normalise the SFH to 1 solar mass formed at the
# end of the 'age' lapse.
if norm:
sfh = sfh / (1.e6 * sfh.sum())
# We convolve the mass table with the SFH for each kind of mass.
# Because of the way numpy full convolution work, we must take in
# the result the first slice with the length of the age grid. The
# SFH is multiplicated by step*1e9 to convert it in Msun/Gyr
conv_masses = [
np.convolve(table, sfh * step * 1.e9)[0:len(time_grid)]
for table in mass_table]
conv_masses = np.array(conv_masses)
# We then interpolate the convolved mass table to the given age
# and return it.
masses = interpolate.interp1d(time_grid, conv_masses)(age)
# We convolve the spectrum table with the SFH.
conv_spectra = [
np.convolve(table, sfh * step * 1e9)[0:len(time_grid)]
for table in spec_table]
conv_spectra = np.array(conv_spectra)
spectra = interpolate.interp1d(time_grid, conv_spectra)(age)
# 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.)
# Step between two item in the time grid in Gyr
step = self.time_grid[1] - self.time_grid[0]
# If needed, we normalise the SFH to 1 solar mass produced.
if norm:
sfh_sfr = sfh_sfr / np.trapz(sfh_sfr * 1e9,
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.
mass_table = self.mass_table[:, :idx + 1]
spec_table = self.spec_table[:, :idx + 1]
# The 1.e9 * step is because the SFH is in solar mass per year.
masses = 1.e9 * step * np.dot(mass_table, sfh_sfr[::-1])
spectra = 1.e9 * step * np.dot(spec_table, sfh_sfr[::-1])
return masses, spectra
......@@ -12,7 +12,6 @@ import numpy as np
from . import common
from pcigale.data import Database
# Time lapse used to compute the average star formation rate. We use a
# constant to keep it easily changeable for advanced user while limiting the
# number of parameters. The value is in Gyr.
......@@ -58,28 +57,22 @@ class Module(common.SEDCreationModule):
None
),
'sfh': (
'numpy array of floats',
"Star formation history (in solar mass per year) as an array "
"based on the age grid used for the Maraston 2005 SSP, i.e. "
"np.arange(1e-3,13.701,1e-3).",
None
),
'oldest_age': (
'float',
"Age of the oldest stars in the galaxy in Gyr.",
'string',
"Name of the key in the SED info dictionary, associated with the "
"SFH tuple (time, sfr). The SFH must be normalised to 1 solar "
"mass produced.",
None
),
'separation_age': (
'float',
"Age [Gyr] of the separation between the young and the old star "
"populations. The default value in 10^7 years (0.01 Gyr). Set "
"to 0 not to differentiate ages.",
"to 0 not to differentiate ages (only an old population).",
0.01
)
}
out_parameter_list = {
'galaxy_age': 'Age (in Gyr) of the oldest stars in the galaxy.',
'sfr': 'Instantaneous Star Formation Rate in solar mass per year, '
'at the age of the galaxy.',
'average_sfr': 'Average SFR in the last 0.1 Gyr (default) of the '
......@@ -90,7 +83,6 @@ class Module(common.SEDCreationModule):
'mass_neutrino': 'Mass of neutrino stars in solar mass.',
'mass_black_hole': 'Mass of black holes in solar mass.',
'mass_turn_off': 'Mass in the turn-off in solar mass.',
'age': 'Age of the oldest stars in the galaxy.',
'old_young_separation_age': 'Age (in Gyr) separating the old and '
'the young star populations (0 if there '
'is only one population).',
......@@ -132,58 +124,36 @@ class Module(common.SEDCreationModule):
"""
imf = self.parameters['imf']
metallicity = self.parameters['metallicity']
sfh = np.copy(self.parameters['sfh'])
oldest_age = self.parameters['oldest_age']
separation_age = self.parameters['separation_age']
imf = self.parameters["imf"]
metallicity = self.parameters["metallicity"]
separation_age = self.parameters["separation_age"]
sfh_time, sfh_sfr = np.copy(sed.info[parameters["sfh"]])
# Age of the galaxy at each time of the SFH
sfh_age = np.max(sfh_time) - sfh_time
# First, we take the SSP out of the database.
database = Database()
ssp = database.get_ssp_m2005(imf, metallicity)
database.session.close_all()
# We check that the sfh array length corresponds to the one of the
# SSP age grid.
if len(ssp.time_grid) != len(sfh):
raise ValueError("The SFH array must be based on the same age "
"grid as the Maraston 2005 SSPs.")
# We limit the star formation history to the age of the oldest stars
# and set SFR=0 for every age after (we need to keep the same age
# grid for the computations.)
sfh[ssp.time_grid > oldest_age] = 0
# The age of the galaxy is taken from the age grid.
galaxy_age = np.max(ssp.time_grid[ssp.time_grid <= oldest_age])
# We normalise the SFH to have 1 solar mass formed at the age of the
# galaxy.
sfh = sfh / np.trapz(sfh * 1e9, ssp.time_grid)
# First, we process the old population (age greater than the
# separation age).
old_sfh = np.copy(sfh)
old_sfh[(galaxy_age - ssp.time_grid) < separation_age] = 0
old_masses, old_spectrum = ssp.convolve(old_sfh, galaxy_age)
# If there is a separation between young and old galaxies, we process
# the young population else we set everything to 0 for the young
# population (always having a young population can be help full e.g.
# to probe various separation ages including 0).
if separation_age:
young_sfh = np.copy(sfh)
young_sfh[(galaxy_age - ssp.time_grid) >= separation_age] = 0
young_masses, young_spectrum = ssp.convolve(young_sfh, galaxy_age)
else:
young_masses = np.zeros(6)
young_spectrum = np.zeros(len(ssp.wavelength_grid))
# First, we process the young population (age lower than the
# separation age.)
young_sfh = np.copy(sfh_sfr)
young_sfh[sfh_age <= separation_age] = 0
young_masses, young_spectrum = ssp.convolve(sfh_time, young_sfh)
# Then, we process the old population. If the SFH is shorter than the
# separation age then all the arrays will consist only of 0.
old_sfh = np.copy(sfh_sfr)
old_sfh[sfh_age > separation_age] = 0
old_masses, old_spectrum = ssp.convolve(sfh_time, old_sfh)
# SFR of the galaxy
sfr = sfh[ssp.time_grid == galaxy_age][0]
sfr = sfh_sfr[len(sfh_sfr) - 1]
# Average SFR on the last AV_LAPSE Gyr of its history
average_sfr = np.mean(sfh[(ssp.time_grid <= oldest_age) &
(ssp.time_grid >= (oldest_age - AV_LAPSE))])
average_sfr = np.mean(sfh_sfr[sfh_age <= AV_LAPSE])
# Base name for adding information to the SED.
name = self.name or 'm2005_sfh'
......@@ -192,7 +162,6 @@ class Module(common.SEDCreationModule):
sed.add_info('imf', imf)
sed.add_info('metallicity', metallicity)
sed.add_info('age', oldest_age)
sed.add_info('old_young_separation_age', separation_age)
sed.add_info('sfr', sfr, True)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment