Commit 3a9d8b8a authored by Médéric Boquien's avatar Médéric Boquien

Add the necessary bits in order to allow to fit directly pure BC03 SSP.

parent b487455e
......@@ -20,7 +20,7 @@ import numpy as np
from scipy import interpolate
import scipy.constants as cst
from astropy.table import Table
from pcigale.data import (Database, Filter, M2005, BC03, Fritz2006,
from pcigale.data import (Database, Filter, M2005, BC03, BC03_SSP, Fritz2006,
Dale2014, DL2007, DL2014, NebularLines,
NebularContinuum, Schreiber2016, THEMIS)
......@@ -410,6 +410,48 @@ def build_bc2003(base, res):
ssp_lumin
))
def build_bc2003_ssp(base, res):
bc03_dir = os.path.join(os.path.dirname(__file__), 'bc03/')
# 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"]):
ssp_filename = "{}bc2003_{}_{}_{}_ssp.ised_ASCII".format(bc03_dir, res,
key, imf)
color3_filename = "{}bc2003_lr_{}_{}_ssp.3color".format(bc03_dir, key,
imf)
color4_filename = "{}bc2003_lr_{}_{}_ssp.4color".format(bc03_dir, key,
imf)
print("Importing {}...".format(ssp_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(10 ** color3_table[5]) # NLy
color_table = np.array(color_table)
ssp_time, ssp_wave, ssp_lumin = read_bc03_ssp(ssp_filename)
base.add_bc03_ssp(BC03_SSP(
imf,
metallicity[key],
ssp_time,
ssp_wave,
color_table,
ssp_lumin
))
def build_dale2014(base):
models = []
......@@ -860,6 +902,7 @@ def build_base(bc03res='lr'):
print('#' * 78)
print("3- Importing Bruzual and Charlot 2003 SSP\n")
build_bc2003_ssp(base, bc03res)
build_bc2003(base, bc03res)
print("\nDONE\n")
print('#' * 78)
......
......@@ -24,6 +24,7 @@ import numpy as np
from .filters import Filter
from .m2005 import M2005
from .bc03 import BC03
from .bc03_ssp import BC03_SSP
from .dale2014 import Dale2014
from .dl2007 import DL2007
from .dl2014 import DL2014
......@@ -116,6 +117,28 @@ class _BC03(BASE):
self.spec_table = ssp.spec_table
class _BC03_SSP(BASE):
"""Storage for Bruzual and Charlot 2003 SSP
"""
__tablename__ = "bc03_ssp"
imf = Column(String, primary_key=True)
metallicity = Column(Float, primary_key=True)
time_grid = Column(PickleType)
wavelength_grid = Column(PickleType)
info_table = Column(PickleType)
spec_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.info_table = ssp.info_table
self.spec_table = ssp.spec_table
class _Dale2014(BASE):
"""Storage for Dale et al (2014) infra-red templates
"""
......@@ -444,6 +467,71 @@ class Database(object):
"""
return self._get_parameters(_BC03)
def add_bc03_ssp(self, ssp_bc03):
"""
Add a Bruzual and Charlot 2003 SSP to pcigale database
Parameters
----------
ssp: pcigale.data.SspBC03
"""
if self.is_writable:
ssp = _BC03_SSP(ssp_bc03)
self.session.add(ssp)
try:
self.session.commit()
except exc.IntegrityError:
self.session.rollback()
raise DatabaseInsertError('The SSP is already in the base.')
else:
raise Exception('The database is not writable.')
def get_bc03_ssp(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.BC03
The BC03 object.
Raises
------
DatabaseLookupError: if the requested SSP is not in the database.
"""
result = self.session.query(_BC03_SSP)\
.filter(_BC03_SSP.imf == imf)\
.filter(_BC03_SSP.metallicity == metallicity)\
.first()
if result:
return BC03_SSP(result.imf, result.metallicity, result.time_grid,
result.wavelength_grid, result.info_table,
result.spec_table)
else:
raise DatabaseLookupError(
"The BC03 SSP for imf <{0}> and metallicity <{1}> is not in "
"the database.".format(imf, metallicity))
def get_bc03_ssp_parameters(self):
"""Get parameters for the Bruzual & Charlot 2003 stellar models.
Returns
-------
paramaters: dictionary
dictionary of parameters and their values
"""
return self._get_parameters(_BC03_SSP)
def add_dl2007(self, models):
"""
Add a list of Draine and Li (2007) models to the database.
......
# -*- 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
import numpy as np
class BC03_SSP(object):
def __init__(self, imf, metallicity, time_grid, wavelength_grid,
info_table, spec_table):
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.info_table = info_table
self.spec_table = spec_table
"""
Bruzual and Charlot (2003) stellar emission module for an SSP
=============================================================
This module implements the Bruzual and Charlot (2003) Single Stellar
Populations.
"""
from collections import OrderedDict
import numpy as np
from . import SedModule
from ..data import Database
class BC03SSP(SedModule):
parameter_list = OrderedDict([
("imf", (
"cigale_list(dtype=int, options=0. & 1.)",
"Initial mass function: 0 (Salpeter) or 1 (Chabrier).",
0
)),
("metallicity", (
"cigale_list(options=0.0001 & 0.0004 & 0.004 & 0.008 & 0.02 & "
"0.05)",
"Metalicity. Possible values are: 0.0001, 0.0004, 0.004, 0.008, "
"0.02, 0.05.",
0.02
)),
("separation_age", (
"cigale_list(dtype=int, minvalue=0)",
"Age [Myr] of the separation between the young and the old star "
"populations. The default value in 10^7 years (10 Myr). Set "
"to 0 not to differentiate ages (only an old population).",
10
))
])
def _init_code(self):
"""Read the SSP from the database."""
self.imf = int(self.parameters["imf"])
self.metallicity = float(self.parameters["metallicity"])
self.separation_age = int(self.parameters["separation_age"])
with Database() as database:
if self.imf == 0:
self.ssp = database.get_bc03_ssp('salp', self.metallicity)
elif self.imf == 1:
self.ssp = database.get_bc03_ssp('chab', self.metallicity)
else:
raise Exception("IMF #{} unknown".format(self.imf))
def process(self, sed):
"""Add the convolution of a Bruzual and Charlot SSP to the SED
Parameters
----------
sed: pcigale.sed.SED
SED object.
"""
if 'ssp.index' in sed.info:
index = sed.info['ssp.index']
else:
raise Exception('The stellar models do not correspond to pure SSP.')
if self.ssp.time_grid[index] <= self.separation_age:
spec_young = self.ssp.spec_table[:, index]
info_young = self.ssp.info_table[:, index]
spec_old = np.zeros_like(spec_young)
info_old = np.zeros_like(info_young)
else:
spec_old = self.ssp.spec_table[:, index]
info_old = self.ssp.info_table[:, index]
spec_young = np.zeros_like(spec_old)
info_young = np.zeros_like(info_old)
info_all = info_young + info_old
info_young = dict(zip(["m_star", "m_gas", "n_ly"], info_young))
info_old = dict(zip(["m_star", "m_gas", "n_ly"], info_old))
info_all = dict(zip(["m_star", "m_gas", "n_ly"], info_all))
# We compute the Lyman continuum luminosity as it is important to
# compute the energy absorbed by the dust before ionising gas.
wave = self.ssp.wavelength_grid
w = np.where(wave <= 91.1)
lum_lyc_young, lum_lyc_old = np.trapz([spec_young[w], spec_old[w]],
wave[w])
# We do similarly for the total stellar luminosity
lum_young, lum_old = np.trapz([spec_young, spec_old], wave)
sed.add_module(self.name, self.parameters)
sed.add_info("stellar.imf", self.imf)
sed.add_info("stellar.metallicity", self.metallicity)
sed.add_info("stellar.old_young_separation_age", self.separation_age)
sed.add_info("stellar.age", self.ssp.time_grid[index])
sed.add_info("stellar.m_star_young", info_young["m_star"], True)
sed.add_info("stellar.m_gas_young", info_young["m_gas"], True)
sed.add_info("stellar.n_ly_young", info_young["n_ly"], True)
sed.add_info("stellar.lum_ly_young", lum_lyc_young, True)
sed.add_info("stellar.lum_young", lum_young, True)
sed.add_info("stellar.m_star_old", info_old["m_star"], True)
sed.add_info("stellar.m_gas_old", info_old["m_gas"], True)
sed.add_info("stellar.n_ly_old", info_old["n_ly"], True)
sed.add_info("stellar.lum_ly_old", lum_lyc_old, True)
sed.add_info("stellar.lum_old", lum_old, True)
sed.add_info("stellar.m_star", info_all["m_star"], True)
sed.add_info("stellar.m_gas", info_all["m_gas"], True)
sed.add_info("stellar.n_ly", info_all["n_ly"], True)
sed.add_info("stellar.lum_ly", lum_lyc_young + lum_lyc_old, True)
sed.add_info("stellar.lum", lum_young + lum_old, True)
sed.add_contribution("stellar.old", wave, spec_old)
sed.add_contribution("stellar.young", wave, spec_young)
# SedModule to be returned by get_module
Module = BC03SSP
"""
Simple module to reduce the SFH to a single SSP
===========================================================
This module implements a star formation history (SFH) through a single SSP.
"""
from collections import OrderedDict
import numpy as np
from . import SedModule
class SSP(SedModule):
"""Instantaneous burst corresponding to a model-provided SSP
This module sets the SED star formation history (SFH) as a single stellar
population
"""
parameter_list = OrderedDict([
("index", (
"cigale_list(dtype=int, minvalue=0)",
"Index of the SSP to use.",
0
))
])
def _init_code(self):
self.index = int(self.parameters["index"])
def process(self, sed):
"""Add a double decreasing exponential Star Formation History.
Parameters
----------
sed: pcigale.sed.SED object
"""
sed.add_module(self.name, self.parameters)
# Add the sfh and the output parameters to the SED.
sed.sfh = np.array([0.])
sed.add_info("ssp.index", self.index)
# SedModule to be returned by get_module
Module = SSP
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