Commit 20f67dd5 authored by BURGARELLA Denis's avatar BURGARELLA Denis Committed by Yannick Roehlly
Browse files

New SED creation modules

- mbb: modified black body
- param: parameterless module to compute some parameters and add them to
  the information dictionary
- sfhcomb_constant and sfhcomb_delayed: two star formation histories
parent e9487032
# -*- coding: utf-8 -*-
# Copyright (C) 2015 Centre de données Astrophysiques de Marseille
# Copyright (C) 2015 Institute of Astronomy, University of Cambridge
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
# Author: Denis Burgarella
"""
Modified Black Body module
======================================
This module implements a modified black body (MBB). This MMB can be
a stand-alone IR emission or it can come as an additional component.
The energy balance can be set to check the presence of a component
or not to account for the emission of a region completely embedded
in dust and not visible in the wavelength range.
"""
import numpy as np
import scipy.constants as cst
from collections import OrderedDict
from . import CreationModule
class MBB(CreationModule):
"""One modified black body IR re-emission
Given an amount of attenuation (e.g. resulting from the action of a dust
attenuation module) this module normalises MBB plus any previous IR
contribution to this amount of energy. The final SED allows to keep the
energy balance or not..
"""
parameter_list = OrderedDict([
("epsilon_mbb", (
"float",
"Fraction [>= Ø] of L_dust(energy balance) in the MBB",
0.5
)),
("t_mbb", (
"float",
"Temperature of black body in K.",
50.
)),
("beta_mbb", (
"float",
"Emissivity index of modified black body.",
1.5
)),
("energy_balance", (
"boolean",
"Energy balance checked?"
"If False, Lum[MBB] not taken into account in energy balance",
False
)),
])
out_parameter_list = OrderedDict([
("t_mbb", "Temperature of the modified black body in K."),
("beta_mbb", "Emissivity index of the modified black body."),
("epsilon_mbb", "Fraction of L_dust in the modified black body.")
])
def _init_code(self):
"""Build the model for a given set of parameters."""
epsilon = float(self.parameters["epsilon_mbb"])
T = float(self.parameters["t_mbb"])
beta = float(self.parameters["beta_mbb"])
if epsilon < 0:
epsilon = 0.
print("epsilon_mbb must >= 0, we set epsilon_mbb = 0.0")
# We define various constants necessary to compute the model. For
# consistency, we define the speed of light in nm s¯¹ rather than in
# m s¯¹.
c = cst.c * 1e9
lambda_0 = 200e3
self.wave = np.logspace(3., 6., 1000.)
conv = c / (self.wave * self.wave)
self.lumin_mbb = (conv * (1. - np.exp(-(lambda_0 / self.wave)
** beta)) * (c / self.wave) ** 3. / (np.exp(
cst.h * c / (self.wave * cst.k * T)) - 1.))
# TODO, save the right normalisation factor to retrieve the dust mass
norm = np.trapz(self.lumin_mbb, x=self.wave)
self.lumin_mbb /= norm
def process(self, sed):
"""Add the IR re-emission contributions.
Parameters
----------
sed: pcigale.sed.SED object
"""
if 'dust.luminosity' not in sed.info.keys():
sed.add_info('dust.luminosity', 1., True)
luminosity = sed.info['dust.luminosity']
sed.add_module(self.name, self.parameters)
sed.add_info("dust.t_mbb", self.parameters["t_mbb"])
sed.add_info("dust.beta_mbb", self.parameters["beta_mbb"])
sed.add_info("dust.epsilon_mbb", self.parameters["epsilon_mbb"])
epsilon = float(self.parameters["epsilon_mbb"])
energy_balance = (self.parameters["energy_balance"].lower() == "true")
if energy_balance:
# Since we can have another contribution to L_dust and the modified black body
# enters into the energy budget, energy balance, we have to save a new negative
# component for each one, previously existing as:
# luminosity_other_balance = -luminosity_other * (1-epsilon);
# epsilon being the contribution of the present MBB to L_dust.
other_dust_contributions = [contrib for contrib in sed.contribution_names
if "dust" in contrib]
for item in other_dust_contributions:
item_balance = item + '_balance'
lumin = sed.get_lumin_contribution(item)
wavelength = sed.wavelength_grid
sed.add_info(item_balance, 1., True)
sed.add_contribution(item_balance, wavelength, -lumin * epsilon)
# We read again the dust luminosities to check everything is OK.
# lumin2 = sed.get_lumin_contribution(item_balance)
# print("residual", lumin, lumin2)
# If the modified black body does not enter into the energy budget,
# we do not change the luminosity of other dust contributions.
# The luminosity of the modified black body L_MBB = epsilon * L_dust.
# The total dust luminosity will be : L_dust(inside energy budget) + L_MBB.
# We add the contribution of the MBB to L_dust.
sed.add_contribution('dust.mbb', self.wave,
luminosity * epsilon * self.lumin_mbb)
#
# CreationModule to be returned by get_module
Module = MBB
# -*- coding: utf-8 -*-
# Copyright (C) 2015 Laboratoire d'Astrophysique de Marseille
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
# Author: Denis Burgarella
"""
Module that estimates other parameters, e.g., UV slope, Lick indices, etc.
==========================================================================
This module estimates additional parameters of interest
close to the observation, e.g., the ultraviolet slope (beta), the rest-frame
far-ultraviolet luminosity, any type of indices (Lick), etc.
This module will be the last one (after redshifting) to account
for all the physical processes at play to build the received total emission.
"""
import numpy as np
from collections import OrderedDict
from . import CreationModule
class Param(CreationModule):
"""Other parameters
This module does not need any input.
It computes some parameters from the models:
beta_UV, Lum_UV, etc.
"""
parameter_list = OrderedDict()
out_parameter_list = OrderedDict([
("beta_calz94", "The UV slope in the range 125-180&240-260nm"),
("FUV_luminosity", "The rest-frame FUV luminosity"),
("D_4000", "The D_4000 index")
])
def process(self, sed):
"""Computes the parameters for each model.
Parameters
----------
sed: pcigale.sed.SED object
"""
# Retrieve the final computed SED using all the previous modules
# including the IGM and the redshifting. In other words,
# this module must be the last one. Note that it does require
# an SFH and an SSP module but nothing else (except redshifting)
redshift = sed.info['redshift']
# Wavelengths are in nanometers.
wavelength = sed.wavelength_grid
# Luminosity is is W/nm.
luminosity = sed.luminosity
# Attenuated (observed) UV slopes beta
# as defined in Calzetti et al. (1994, ApJ 429, 582, Tab. 2)
# that excludes the 217.5 nm bump wavelength range and other spectral features
w_calz94 = np.where((wavelength >= 126.8 * (1. + redshift)) &
(wavelength <= 128.4 * (1. + redshift)) |
(wavelength >= 130.9 * (1. + redshift)) &
(wavelength <= 131.6 * (1. + redshift)) |
(wavelength >= 134.2 * (1. + redshift)) &
(wavelength <= 137.1 * (1. + redshift)) |
(wavelength >= 140.7 * (1. + redshift)) &
(wavelength <= 151.5 * (1. + redshift)) |
(wavelength >= 156.2 * (1. + redshift)) &
(wavelength <= 158.3 * (1. + redshift)) |
(wavelength >= 167.7 * (1. + redshift)) &
(wavelength <= 174.0 * (1. + redshift)) |
(wavelength >= 176.0 * (1. + redshift)) &
(wavelength <= 183.3 * (1. + redshift)) |
(wavelength >= 186.6 * (1. + redshift)) &
(wavelength <= 189.0 * (1. + redshift)) |
(wavelength >= 193.0 * (1. + redshift)) &
(wavelength <= 195.0 * (1. + redshift)) |
(wavelength >= 240.0 * (1. + redshift)) &
(wavelength <= 258.0 * (1. + redshift)))
# Attenuated (observed) FUV luminosity L_FUV, in the GALEX FUV band,
# i.e., lambda_eff = 152.8 nm and effective bandwidth = 11.4 nm
w_FUV = np.where((wavelength >= (152.8 - 11.4) * (1. + redshift)) &
(wavelength <= (152.8 + 11.4) * (1. + redshift)) )
# Strength of the D_4000 break using Balogh et al. (1999, ApJ 527, 54), i.e.,
# ratio of the flux in the red continuum to that in the blue continuum:
# Blue continuum: 385.0-395.0 nm& red continuum: 410.0-410.0 nm.
w_D4000blue = np.where((wavelength >= 385.0 * (1. + redshift)) &
(wavelength <= 395.0 * (1. + redshift)))
w_D4000red = np.where((wavelength >= 400.0 * (1. + redshift)) &
(wavelength <= 410.0 * (1. + redshift)))
regression_calz94 = np.polyfit(np.log10(10.*wavelength[w_calz94]),
np.log10(1e7/10.*luminosity[w_calz94]), 1)
beta_calz94 = regression_calz94[0]
L_FUV = 152.8*(1. + redshift)*np.mean(luminosity[w_FUV])
D_4000 = np.mean(luminosity[w_D4000red]) / np.mean(luminosity[w_D4000blue])
#plt.scatter(np.log10(10.*wavelength[w_calz94]),
# np.log10(1e7/10.*luminosity[w_calz94]),
# alpha=0.5, color="b", marker="v")
#plt.plot(np.log10(10.*wavelength[w_calz94]),
# regression_calz94[0]*np.log10(10.*wavelength[w_calz94])+regression_calz94[1],
# alpha=0.5, color="r", linestyle="-")
#plt.show()
sed.add_info("param.beta_calz94", beta_calz94)
sed.add_info("param.FUV_luminosity", L_FUV, True)
sed.add_info("param.D_4000", D_4000)
# CreationModule to be returned by get_module
Module = Param
# -*- 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: Denis Burgarella
"""
SFH resembling a Dirac Comb formed by several regularly-spaced short constant SF events.
========================================================================================
This module implements a star formation history (SFH) formed by several
regularly-spaced short constant SF events.
"""
import numpy as np
from collections import OrderedDict
from . import CreationModule
# Time lapse used in the age grid in Myr. If should be consistent with the
# time lapse in the SSP modules.
AGE_LAPSE = 1
class SfhComb(CreationModule):
"""Several regularly-spaced short constant SF events
This module sets the SED star formation history (SFH) as a combination of
several regularly-spaced short constant SF events.
"""
parameter_list = OrderedDict([
("N_events", (
"integer",
"Number of individual star formation events. ",
5
)),
("t_duration", (
"integer",
"Length of each individual star formation event [Myr]. "
"Warning: if too long => continuous SFH",
20
)),
("age", (
"integer",
"Age of the main stellar population in the galaxy in Myr."
"The precision is 1 Myr.",
1000
)),
("age_last", (
"integer",
"Time since the end of the last SF event [Myr]. "
"Warning: Depending on the parameters, the last burst might not be finished. "
"The precision is 1 Myr.",
100
)),
("sfr_A", (
"float",
"Strength of each SF event in M_sun/yr (useful if not normalised)",
1.
)),
("normalise", (
"boolean",
"Normalise the SFH to produce one solar mass.",
"True"
)),
])
out_parameter_list = OrderedDict([
("N_events", "Number of individual star formation events"),
("t_duration", "Length of each individual star formation event "
"in Myr."),
("age", "Age of the stellar population in the galaxy in Myr."),
("age_last", "Time since the end of the last SF event in Myr."),
("sfr_A", "Height of each SF event in M_sun/yr.")
])
def process(self, sed):
"""Add a star formation history formed by several regularly-spaced short SF events.
** Parameters **
sed: pcigale.sed.SED object
"""
N_events = int(self.parameters["N_events"])
t_duration = int(self.parameters["t_duration"])
age = int(self.parameters["age"])
age_last = int(self.parameters["age_last"])
sfr_A = int(self.parameters["sfr_A"])
normalise = (self.parameters["normalise"].lower() == "true")
# Time grid and age. If needed, the age is rounded to the inferior Myr
time_grid = np.arange(AGE_LAPSE, age + AGE_LAPSE, AGE_LAPSE)
age = np.max(time_grid)
sfr = np.zeros(len(time_grid))
# Build the Dirac comb of events.
# Ages of the sequence of SF events in Myr,
# N_events over "age" with the last one starting at "age_last"
time_event = np.linspace(0., age-age_last, num = N_events)
#print("time_event", time_event)
for i_time in range(N_events):
sfr[time_event[i_time]:time_event[i_time]+t_duration] = sfr_A
#print("sfr", sfr)
# Compute the galaxy mass and normalise the SFH to 1 solar mass
# produced if asked to.
galaxy_mass = np.trapz(sfr * 1e6, time_grid)
if normalise:
sfr = sfr / galaxy_mass
galaxy_mass = 1.
sed.add_module(self.name, self.parameters)
# Add the sfh and the output parameters to the SED.
sed.sfh = (time_grid, sfr)
sed.add_info("galaxy_mass", galaxy_mass, True)
sed.add_info("sfh.N_events", N_events)
sed.add_info("sfh.t_duration", t_duration)
sed.add_info("sfh.age_last", age_last)
#sed.add_info("age", age)
# CreationModule to be returned by get_module
Module = SfhComb
# -*- 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: Denis Burgarella
"""
SFH resembling a Dirac Comb with regularly-spaced delayed-SFH events
=====================================================================
This module implements a star formation history (SFH) formed by several
regularly-spaced short regularly-spaced delayed-SFH SF events.
"""
import numpy as np
from collections import OrderedDict
from . import CreationModule
# Time lapse used in the age grid in Myr. If should be consistent with the
# time lapse in the SSP modules.
AGE_LAPSE = 1
class SfhComb(CreationModule):
"""Several regularly-spaced short delayed-SFH SF events
This module sets the SED star formation history (SFH) as a combination of
several regularly-spaced short SF events.
"""
parameter_list = OrderedDict([
("N_events", (
"integer",
"Number of individual star formation events. ",
5
)),
("tau_events", (
"float",
"e-folding time of all short events in Myr.",
20.
)),
("age", (
"integer",
"Age of the main stellar population in the galaxy in Myr."
"The precision is 1 Myr.",
1000
)),
("age_last", (
"integer",
"Time since the end of the last SF event [Myr]. "
"Warning: Depending on the parameters, the last burst might not be finished. "
"The precision is 1 Myr.",
100
)),
("sfr_A", (
"float",
"Multiplicative factor controlling the amplitude of SFR (valid for each event).",
1.
)),
("normalise", (
"boolean",
"Normalise the SFH to produce one solar mass.",
"True"
)),
])
out_parameter_list = OrderedDict([
("N_events", "Number of individual star formation events"),
("tau_events", "e-folding time of all short events in Myr.."),
("age", "Age of the stellar population in the galaxy in Myr."),
("age_last", "Time since the end of the last SF event in Myr."),
("sfr_A", "Height of each SF event in M_sun/yr.")
])
def process(self, sed):
"""Add a star formation history formed by several
regularly-spaced delayed-SFH short SF events.
** Parameters **
sed: pcigale.sed.SED object
"""
N_events = int(self.parameters["N_events"])
tau_events = float(self.parameters["tau_events"])
age = int(self.parameters["age"])
age_last = int(self.parameters["age_last"])
sfr_A = int(self.parameters["sfr_A"])
normalise = (self.parameters["normalise"].lower() == "true")
# Time grid and age. If needed, the age is rounded to the inferior Myr
time_grid = np.arange(AGE_LAPSE, age + AGE_LAPSE, AGE_LAPSE)
age = np.max(time_grid)
sfr = np.zeros(len(time_grid))
# Build the Dirac comb of events.
# Ages of the sequence of SF events in Myr,
# N_events over "age" with the last one starting at "age_last"
time_event = np.linspace(0., age-age_last, num = N_events)
#print("time_events", time_event)
# Add each delayed SFH
for i_time in range(N_events):
sfr[time_event[i_time]:] = \
sfr_A * (time_grid[time_event[i_time]:]-time_event[i_time]) / tau_events**2 * \
np.exp(-(time_grid[time_event[i_time]:]-time_event[i_time]) / tau_events)
# Compute the galaxy mass and normalise the SFH to 1 solar mass
# produced if asked to.
galaxy_mass = np.trapz(sfr * 1e6, time_grid)
if normalise:
sfr = sfr / galaxy_mass
galaxy_mass = 1.
sed.add_module(self.name, self.parameters)
# Add the sfh and the output parameters to the SED.
sed.sfh = (time_grid, sfr)
sed.add_info("galaxy_mass", galaxy_mass, True)
sed.add_info("sfh.N_events", N_events)
sed.add_info("sfh.tau_events", tau_events)
sed.add_info("sfh.age_last", age_last)
#sed.add_info("age", age)
# CreationModule to be returned by get_module
Module = SfhComb
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