Commit 6b56ce3a authored by Médéric Boquien's avatar Médéric Boquien

Initial implementation to handle the Yggdrasil SSP.

parent 3a9d8b8a
......@@ -22,7 +22,8 @@ import scipy.constants as cst
from astropy.table import Table
from pcigale.data import (Database, Filter, M2005, BC03, BC03_SSP, Fritz2006,
Dale2014, DL2007, DL2014, NebularLines,
NebularContinuum, Schreiber2016, THEMIS)
NebularContinuum, Schreiber2016, THEMIS,
Yggdrasil_SSP)
def read_bc03_ssp(filename):
......@@ -453,6 +454,42 @@ def build_bc2003_ssp(base, res):
ssp_lumin
))
def build_yggdrasil_ssp(base):
yggdrasil_dir = os.path.join(os.path.dirname(__file__), 'yggdrasil/')
# Metallicities associated to each key
metallicities = ["0.004", "0.008", "0.02"]
for Z in metallicities:
filename = f"{yggdrasil_dir}Z={Z}_kroupa_IMF_fcov_0.5_SFR_inst_Spectra"
print(f"Importing {filename}...")
with open(filename) as f:
specallages = "".join(f.readlines()[2::]).split('\n\n')[:-1]
for i in range(len(specallages)):
specallages[i] = specallages[i].split('\n')
ssp_time = np.array([float(item[0].split()[-1]) for item in specallages])
ssp_info = np.array([float(item[2].split()[-1]) for item in specallages])
# Normalisation from 10⁶ to 1 Msun
ssp_info *= 1e-6
minwlsize = int(np.min([float(item[3].split()[-1]) for item in specallages]))
ssp_wave = np.array([float(line.split(' ')[0]) for line in specallages[0][7:]])[-minwlsize:]
# Conversion from Å to nm
ssp_wave *= .1
ssp_lumin = np.empty((minwlsize, ssp_time.size))
for i, spec in enumerate(specallages):
ssp_lumin[:, i] = np.array([float(line.split(' ')[-1]) for line in spec[7:]])[-minwlsize:]
# Conversion from erg/s/Å/10⁶ Msun to W/nm/Msun
ssp_lumin *= 1e-7 * 10 * 1e-6
base.add_yggdrasil_ssp(Yggdrasil_SSP(float(Z), ssp_time, ssp_wave,
ssp_info, ssp_lumin))
def build_dale2014(base):
models = []
dale2014_dir = os.path.join(os.path.dirname(__file__), 'dale2014/')
......@@ -942,6 +979,7 @@ def build_base(bc03res='lr'):
print("\nDONE\n")
print('#' * 78)
build_yggdrasil_ssp(base)
base.session.close_all()
......
......@@ -33,6 +33,7 @@ from .nebular_continuum import NebularContinuum
from .nebular_lines import NebularLines
from .schreiber2016 import Schreiber2016
from .themis import THEMIS
from .yggdrasil_ssp import Yggdrasil_SSP
DATABASE_FILE = pkg_resources.resource_filename(__name__, 'data.db')
......@@ -139,6 +140,26 @@ class _BC03_SSP(BASE):
self.spec_table = ssp.spec_table
class _Yggdrasil_SSP(BASE):
"""Storage for the Yggdrasil SSP
"""
__tablename__ = "yggdrasil_ssp"
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.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
"""
......@@ -532,6 +553,42 @@ class Database(object):
"""
return self._get_parameters(_BC03_SSP)
def add_yggdrasil_ssp(self, ssp_yggdrasil):
"""
Add a Bruzual and Charlot 2003 SSP to pcigale database
Parameters
----------
ssp: pcigale.data.SspBC03
"""
if self.is_writable:
ssp = _Yggdrasil_SSP(ssp_yggdrasil)
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_yggdrasil_ssp(self, metallicity):
result = self.session.query(_Yggdrasil_SSP)\
.filter(_Yggdrasil_SSP.metallicity == metallicity)\
.first()
if result:
return Yggdrasil_SSP(result.metallicity, result.time_grid,
result.wavelength_grid, result.info_table,
result.spec_table)
else:
raise DatabaseLookupError(
"The Yggdrasil SSP for metallicity <{}> is not in the "
"database.".format(metallicity))
def get_yggdrasil_ssp_parameters(self):
return self._get_parameters(_Yggdrasil_SSP)
def add_dl2007(self, models):
"""
Add a list of Draine and Li (2007) models to the database.
......
import numpy as np
class Yggdrasil_SSP(object):
def __init__(self, metallicity, time_grid, wavelength_grid,
info_table, spec_table):
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 YggdrasilSSP(SedModule):
parameter_list = OrderedDict([
("metallicity", (
"cigale_list(options=0.004 & 0.008 & 0.02)",
"Metalicity. Possible values are: 0.004, 0.008, and 0.02.",
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.metallicity = float(self.parameters["metallicity"])
self.separation_age = int(self.parameters["separation_age"])
with Database() as database:
self.ssp = database.get_yggdrasil_ssp(self.metallicity)
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
wave = self.ssp.wavelength_grid
lum_young, lum_old = np.trapz([spec_young, spec_old], wave)
sed.add_module(self.name, self.parameters)
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, True)
sed.add_info("stellar.lum_young", lum_young, True)
sed.add_info("stellar.m_star_old", info_old, True)
sed.add_info("stellar.lum_old", lum_old, True)
sed.add_info("stellar.m_star", info_all, 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 = YggdrasilSSP
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