Commit 77681cbc authored by Yannick Roehlly's avatar Yannick Roehlly

Add Bruzual and Charlot (2003) SSP

parent 768dc809
# -*- 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>
......@@ -15,16 +15,130 @@ import sys
import os
sys.path.append(os.path.join(os.path.dirname(__file__), '../'))
import glob
import itertools
import numpy as np
from scipy import interpolate
from pcigale.data import Database, Filter, SspM2005
from pcigale.data import Database, Filter, SspM2005, SspBC03
filters_dir = os.path.join(os.path.dirname(__file__), 'filters/')
m2005_dir = os.path.join(os.path.dirname(__file__), 'maraston2005/')
bc03_dir = os.path.join(os.path.dirname(__file__), 'bc03//')
dh2002_dir = os.path.join(os.path.dirname(__file__), 'dh2002/')
def read_bc03_ssp(filename):
"""Read a Bruzual and Charlot 2003 ASCII SSP file
The ASCII SSP files of Bruzual and Charlot 2003 have se special structure.
A vector is stored with the number of values followed by the values
sepateded by a space (or a carriage return). There are the time vector, 5
(for Chabrier IMF) or 6 lines (for Salpeter IMF) that we don't care of,
then the wavelength vector, then the luminosity vectors, each followed by
a 52 value table, then a bunch of other table of information that are also
in the *colors files.
Parameters
----------
filename : string
Returns
-------
time_grid: numpy 1D array of floats
Vector of the time grid of the SSP in Gyr.
wavelength: numpy 1D array of floats
Vector of the wavelength grid of the SSP in nm.
spectra: numpy 2D array of floats
Array containing the SSP spectra, first axis is the wavelength,
second one is the time.
"""
def file_structure_generator():
"""Generator used to identify table lines in the SSP file
In the SSP file, the vectors are store one next to the other, but
there are 5 informational lines after the time vector. We use this
generator to the if we are on lines to read or not.
"""
if "chab" in filename:
bad_line_number = 5
else:
bad_line_number = 6
yield("data")
for i in range(bad_line_number):
yield("bad")
while True:
yield("data")
file_structure = file_structure_generator()
# Are we in a data line or a bad one.
what_line = file_structure.next()
# Variable conting, in reverse order, the number of value still to
# read for the read vector.
counter = 0
time_grid = []
full_table = []
tmp_table = []
with open(filename) as file_:
# We read the file line by line.
for line in file_:
if what_line == "data":
# If we are in a "data" line, we analyse each number.
for item in line.split():
if counter == 0:
# If counter is 0, then we are not reading a vector
# and the first number is the length of the next
# vector.
counter = int(item)
else:
# If counter > 0, we are currently reading a vector.
tmp_table.append(float(item))
counter -= 1
if counter == 0:
# We reached the end of the vector. If we have not
# yet store the time grid (the first table) we are
# currently reading it.
if time_grid == []:
time_grid = tmp_table[:]
# Else, we store the vector in the full table,
# only if its length is superior to 250 to get rid
# of the 52 item unknown vector and the 221 (time
# grid length) item vectors at the end of the
# file.
elif len(tmp_table) > 250:
full_table.append(tmp_table[:])
tmp_table = []
# If at the end of a line, we have finished reading a vector, it's
# time to change to the next structure context.
if counter == 0:
what_line = file_structure.next()
# The time grid is in year, we want Gyr.
time_grid = np.array(time_grid, dtype=float)
time_grid = time_grid * 1.e-9
# The first "long" vector encountered is the wavelength grid. The value
# are in Ångström, we convert it to nano-meter.
wavelength = np.array(full_table.pop(0), dtype=float)
wavelength = wavelength * 0.1
# The luminosities are in Solar luminosity (3.826.10^33 ergs.s-1) per
# Ångström, we convert it to W/nm.
luminosity = np.array(full_table, dtype=float)
luminosity = luminosity * 3.826e27
# Transposition to have the time in the second axis.
luminosity = luminosity.transpose()
# In the SSP, the time grid begins at 0, but not in the *colors file, so
# we remove t=0 from the SSP.
return time_grid[1:], wavelength, luminosity[:, 1:]
def build_base():
base = Database(writable=True)
base.upgrade_base()
......@@ -126,10 +240,66 @@ def build_base():
print("\nDONE\n")
print('#' * 78)
########################################################################
# Bruzual and Charlot SSP insertion #
########################################################################
print("3- Importing Bruzual and Charlot 2003 SSP\n")
# Time grid (1My to 20Gy with 1My step)
time_grid = np.arange(1e-3, 20., 1e-3)
# Metallicities associated to each key
metallicity = {
"m22": 0.0001,
"m32": 0.0004,
"m42": 0.004,
"m52": 0.008,
"m62": 0.02,
"m72": 0.05
}
for key, imf in itertools.product(metallicity, ["salp", "chab"]):
base_filename = bc03_dir + "bc2003_lr_" + key + "_" + imf + "_ssp"
ssp_filename = base_filename + ".ised_ASCII"
color3_filename = base_filename + ".3color"
color4_filename = base_filename + ".4color"
print("Importing %s..." % base_filename)
# Read the desired information from the color files
color_table = []
color3_table = np.genfromtxt(color3_filename).transpose()
color4_table = np.genfromtxt(color4_filename).transpose()
color_table.append(color4_table[6]) # Mstar
color_table.append(color4_table[7]) # Mgas
color_table.append(color3_table[5]) # NLy
color_table.append(color3_table[1]) # B4000
color_table.append(color3_table[2]) # B4_VN
color_table.append(color3_table[3]) # B4_SDSS
color_table.append(color3_table[4]) # B(912)
color_table = np.array(color_table)
ssp_time, ssp_wave, ssp_lumin = read_bc03_ssp(ssp_filename)
# Regrid the SSP data to the evenly spaced time grid.
color_table = interpolate.interp1d(ssp_time, color_table)(time_grid)
ssp_lumin = interpolate.interp1d(ssp_time,
ssp_lumin)(time_grid)
base.add_ssp_bc03(SspBC03(
imf,
metallicity[key],
time_grid,
ssp_wave,
color_table,
ssp_lumin
))
########################################################################
# Dale and Helou 2002 templates insertion #
########################################################################
print("3- Importing Dale and Helou 2002 templates\n")
print("4- Importing Dale and Helou 2002 templates\n")
# Getting the alpha grid for the templates
dhcal = np.genfromtxt(dh2002_dir + 'dhcal.dat')
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
# -*- 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>
......@@ -23,6 +23,7 @@ from sqlalchemy.ext.declarative import declarative_base
from sqlalchemy.orm import sessionmaker
from .filters import Filter
from .ssp_m2005 import SspM2005
from .ssp_bc03 import SspBC03
from .ir_templates_dh2002 import IrTemplatesDH2002
......@@ -75,6 +76,28 @@ class _SspM2005(BASE):
self.spec_table = ssp.spec_table
class _SspBC03(BASE):
"""Storage for Bruzual and Charlot 2003 SSP
"""
__tablename__ = "bc03"
imf = Column(String, primary_key=True)
metallicity = Column(Float, primary_key=True)
time_grid = Column(PickleType)
wavelength_grid = Column(PickleType)
color_table = Column(PickleType)
lumin_table = Column(PickleType)
def __init__(self, ssp):
self.imf = ssp.imf
self.metallicity = ssp.metallicity
self.time_grid = ssp.time_grid
self.wavelength_grid = ssp.wavelength_grid
self.color_table = ssp.color_table
self.lumin_table = ssp.lumin_table
class _DH2002InfraredTemplates(BASE):
"""Storage for Dale and Helou (2002) infra-red templates
......@@ -166,6 +189,26 @@ class Database(object):
else:
raise StandardError('The database is not writable.')
def add_ssp_bc03(self, ssp_bc03):
"""
Add a Bruzual and Charlot 2003 SSP to pcigale database
Parameters
----------
ssp : pcigale.data.SspBC03
"""
if self.is_writable:
ssp = _SspBC03(ssp_bc03)
self.session.add(ssp)
try:
self.session.commit()
except exc.IntegrityError:
self.session.rollback()
raise StandardError('The SSP is yet in the base.')
else:
raise StandardError('The database is not writable.')
def add_dh2002_infrared_templates(self, data):
"""
Add Dale and Helou (2002) templates the collection.
......@@ -222,6 +265,35 @@ class Database(object):
else:
return None
def get_ssp_bc03(self, imf, metallicity):
"""
Query the database for the Bruzual and Charlot 2003 SSP corresponding
to the given initial mass function and metallicity.
Parameters
----------
imf : string
Initial mass function (salp for Salpeter, chab for Chabrier)
metallicity : float
0.02 for Solar metallicity
Returns
-------
ssp : pcigale.data.SspBC03
The SspBC03 object. If no SSP corresponds to the given imf and
metallicity, returns None.
"""
result = self.session.query(_SspBC03)\
.filter(_SspBC03.imf == imf)\
.filter(_SspBC03.metallicity == metallicity)\
.first()
if result:
return SspBC03(result.imf, result.metallicity, result.time_grid,
result.wavelength_grid, result.color_table,
result.lumin_table)
else:
return None
def get_ssp_m2005(self, imf, metallicity):
"""
Query the database for a Maraston 2005 SSP corresponding to the given
......
# -*- coding: utf-8 -*-
"""
Copyright (C) 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>
"""
import numpy as np
class SspBC03(object):
"""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
----------
imf : string
Initial mass function (IMF): either 'salp' for Salpeter (1955) or
'chab' for Chabrier (2003).
metallicity : float
The metallicity. Possible values are 0.0001, 0.0004, 0.004, 0.008,
0.02 (solar metallicity) and 0.05.
time_grid : array of floats
The time [Gyr] grid used in the colors_table and the lumin_table.
wavelength_grid : array of floats
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
* color_table[2]: log rate of H-ionizing photons (s-1)
* 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
lumin_table : 2 axis array of floats
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
----------
sfh_time : array of floats
Time grid of the star formation history. It must be increasing and
not run beyond 20 Gyr. The SFH will be regrided to the SSP time
grid.
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
-------
wavelength : array of floats
Wavelength grid [nm] for the spectrum
luminosity : array of floats
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
- "n_ly": log rate of H-ionizing photons (s-1)
- "b_4000": Amplitude of 4000 Å break (Bruzual 2003)
- "b4_vn": Amplitude of 4000 Å narrow break (Balogh et al. 1999)
- "b4_sdss" : Amplitude of 4000 Å break (Stoughton et al. 2002)
- "b_912": Amplitude of Lyman discontinuity
"""
# 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.)
# 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.
color_table = self.color_table[:, :idx + 1]
lumin_table = self.lumin_table[:, :idx + 1]
# The 1.e9 * step is because the SFH is in solar mass per year.
color_info = 1.e9 * step * np.dot(color_table, sfh_sfr[::-1])
luminosity = 1.e9 * step * np.dot(lumin_table, sfh_sfr[::-1])
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
This diff is collapsed.
Markdown is supported
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