Commit e40d9852 authored by Médéric Boquien's avatar Médéric Boquien

Modules were a bit inconsistent how the input parameters are used. It was not...

Modules were a bit inconsistent how the input parameters are used. It was not always ensured they were the right type, they were sometimes read several times, and finally they were not always obtained from the _init_code() method. This is an attempt to clean that up as much as possible.
parent 4f8c2043
...@@ -51,15 +51,17 @@ class BC03(CreationModule): ...@@ -51,15 +51,17 @@ class BC03(CreationModule):
def _init_code(self): def _init_code(self):
"""Read the SSP from the database.""" """Read the SSP from the database."""
if self.parameters["imf"] == 0: self.imf = int(self.parameters["imf"])
imf = 'salp' self.metallicity = float(self.parameters["metallicity"])
elif self.parameters["imf"] == 1: self.separation_age = int(self.parameters["separation_age"])
imf = 'chab'
else:
raise Exception("IMF #{} unknown".format(self.parameters["imf"]))
metallicity = float(self.parameters["metallicity"])
with Database() as database: with Database() as database:
self.ssp = database.get_bc03(imf, metallicity) if self.imf == 0:
self.ssp = database.get_bc03('salp', self.metallicity)
elif self.imf == 1:
self.ssp = database.get_bc03('chab', self.metallicity)
else:
raise Exception("IMF #{} unknown".format(self.imf))
def process(self, sed): def process(self, sed):
"""Add the convolution of a Bruzual and Charlot SSP to the SED """Add the convolution of a Bruzual and Charlot SSP to the SED
...@@ -70,9 +72,6 @@ class BC03(CreationModule): ...@@ -70,9 +72,6 @@ class BC03(CreationModule):
SED object. SED object.
""" """
imf = self.parameters["imf"]
metallicity = float(self.parameters["metallicity"])
separation_age = int(self.parameters["separation_age"])
sfh_time, sfh_sfr = sed.sfh sfh_time, sfh_sfr = sed.sfh
ssp = self.ssp ssp = self.ssp
...@@ -82,13 +81,13 @@ class BC03(CreationModule): ...@@ -82,13 +81,13 @@ class BC03(CreationModule):
# First, we process the young population (age lower than the # First, we process the young population (age lower than the
# separation age.) # separation age.)
young_sfh = np.copy(sfh_sfr) young_sfh = np.copy(sfh_sfr)
young_sfh[sfh_age > separation_age] = 0 young_sfh[sfh_age > self.separation_age] = 0
young_wave, young_lumin, young_info = ssp.convolve(sfh_time, young_sfh) young_wave, young_lumin, young_info = ssp.convolve(sfh_time, young_sfh)
# Then, we process the old population. If the SFH is shorter than the # Then, we process the old population. If the SFH is shorter than the
# separation age then all the arrays will consist only of 0. # separation age then all the arrays will consist only of 0.
old_sfh = np.copy(sfh_sfr) old_sfh = np.copy(sfh_sfr)
old_sfh[sfh_age <= separation_age] = 0 old_sfh[sfh_age <= self.separation_age] = 0
old_wave, old_lumin, old_info = ssp.convolve(sfh_time, old_sfh) old_wave, old_lumin, old_info = ssp.convolve(sfh_time, old_sfh)
# We compute the Lyman continuum luminosity as it is important to # We compute the Lyman continuum luminosity as it is important to
...@@ -99,9 +98,9 @@ class BC03(CreationModule): ...@@ -99,9 +98,9 @@ class BC03(CreationModule):
sed.add_module(self.name, self.parameters) sed.add_module(self.name, self.parameters)
sed.add_info("stellar.imf", imf) sed.add_info("stellar.imf", self.imf)
sed.add_info("stellar.metallicity", metallicity) sed.add_info("stellar.metallicity", self.metallicity)
sed.add_info("stellar.old_young_separation_age", separation_age) sed.add_info("stellar.old_young_separation_age", self.separation_age)
sed.add_info("stellar.m_star_young", young_info["m_star"], True) sed.add_info("stellar.m_star_young", young_info["m_star"], True)
sed.add_info("stellar.m_gas_young", young_info["m_gas"], True) sed.add_info("stellar.m_gas_young", young_info["m_gas"], True)
......
...@@ -49,7 +49,8 @@ class Casey2012(CreationModule): ...@@ -49,7 +49,8 @@ class Casey2012(CreationModule):
def _init_code(self): def _init_code(self):
"""Build the model for a given set of parameters.""" """Build the model for a given set of parameters."""
# To compactify the following equations we only assign them to self at
# the end of the method
T = float(self.parameters["temperature"]) T = float(self.parameters["temperature"])
beta = float(self.parameters["beta"]) beta = float(self.parameters["beta"])
alpha = float(self.parameters["alpha"]) alpha = float(self.parameters["alpha"])
...@@ -83,6 +84,10 @@ class Casey2012(CreationModule): ...@@ -83,6 +84,10 @@ class Casey2012(CreationModule):
self.lumin_blackbody /= norm self.lumin_blackbody /= norm
self.lumin = self.lumin_powerlaw + self.lumin_blackbody self.lumin = self.lumin_powerlaw + self.lumin_blackbody
self.temperature = T
self.beta = beta
self.alpha = alpha
def process(self, sed): def process(self, sed):
"""Add the IR re-emission contributions. """Add the IR re-emission contributions.
...@@ -96,9 +101,10 @@ class Casey2012(CreationModule): ...@@ -96,9 +101,10 @@ class Casey2012(CreationModule):
luminosity = sed.info['dust.luminosity'] luminosity = sed.info['dust.luminosity']
sed.add_module(self.name, self.parameters) sed.add_module(self.name, self.parameters)
sed.add_info("dust.temperature", self.parameters["temperature"]) sed.add_info("dust.temperature", self.temperature)
sed.add_info("dust.alpha", self.parameters["alpha"]) sed.add_info("dust.beta", self.beta)
sed.add_info("dust.beta", self.parameters["beta"]) sed.add_info("dust.alpha", self.alpha)
sed.add_contribution('dust.powerlaw', self.wave, sed.add_contribution('dust.powerlaw', self.wave,
luminosity * self.lumin_powerlaw) luminosity * self.lumin_powerlaw)
......
...@@ -60,10 +60,11 @@ class Dale2014(CreationModule): ...@@ -60,10 +60,11 @@ class Dale2014(CreationModule):
the quasar the quasar
The energy attenuated is re-injected in model_sb only. The energy attenuated is re-injected in model_sb only.
""" """
alpha = self.parameters["alpha"] self.fracAGN = float(self.parameters["fracAGN"])
self.alpha = float(self.parameters["alpha"])
with Database() as database: with Database() as database:
self.model_sb = database.get_dale2014(0.00, alpha) self.model_sb = database.get_dale2014(0.00, self.alpha)
self.model_quasar = database.get_dale2014(1.00, 0.0) self.model_quasar = database.get_dale2014(1.00, 0.0)
def process(self, sed): def process(self, sed):
...@@ -79,21 +80,18 @@ class Dale2014(CreationModule): ...@@ -79,21 +80,18 @@ class Dale2014(CreationModule):
sed.add_info('dust.luminosity', 1., True) sed.add_info('dust.luminosity', 1., True)
luminosity = sed.info['dust.luminosity'] luminosity = sed.info['dust.luminosity']
frac_agn = self.parameters["fracAGN"] if self.fracAGN < 1.:
L_AGN = luminosity * (1./(1.-self.fracAGN) - 1.)
if frac_agn < 1.:
L_AGN = luminosity * (1./(1.-frac_agn) - 1.)
else: else:
raise Exception("AGN fraction is exactly 1. Behaviour " raise Exception("AGN fraction is exactly 1. Behaviour undefined.")
"undefined.")
sed.add_module(self.name, self.parameters) sed.add_module(self.name, self.parameters)
sed.add_info("agn.fracAGN_dale2014", self.parameters["fracAGN"]) sed.add_info("agn.fracAGN_dale2014", self.fracAGN)
sed.add_info("dust.alpha", self.parameters["alpha"]) sed.add_info("dust.alpha", self.alpha)
sed.add_contribution('dust', self.model_sb.wave, sed.add_contribution('dust', self.model_sb.wave,
luminosity * self.model_sb.lumin) luminosity * self.model_sb.lumin)
if frac_agn != 0.: if self.fracAGN != 0.:
sed.add_contribution('agn', self.model_quasar.wave, sed.add_contribution('agn', self.model_quasar.wave,
L_AGN * self.model_quasar.lumin) L_AGN * self.model_quasar.lumin)
......
...@@ -63,10 +63,10 @@ class DL2007(CreationModule): ...@@ -63,10 +63,10 @@ class DL2007(CreationModule):
def _init_code(self): def _init_code(self):
"""Get the model out of the database""" """Get the model out of the database"""
self.qpah = self.parameters["qpah"] self.qpah = float(self.parameters["qpah"])
self.umin = self.parameters["umin"] self.umin = float(self.parameters["umin"])
self.umax = self.parameters["umax"] self.umax = float(self.parameters["umax"])
self.gamma = self.parameters["gamma"] self.gamma = float(self.parameters["gamma"])
with Database() as database: with Database() as database:
self.model_minmin = database.get_dl2007(self.qpah, self.umin, self.model_minmin = database.get_dl2007(self.qpah, self.umin,
......
...@@ -66,10 +66,10 @@ class DL2014(CreationModule): ...@@ -66,10 +66,10 @@ class DL2014(CreationModule):
def _init_code(self): def _init_code(self):
"""Get the model out of the database""" """Get the model out of the database"""
self.qpah = self.parameters["qpah"] self.qpah = float(self.parameters["qpah"])
self.umin = self.parameters["umin"] self.umin = float(self.parameters["umin"])
self.alpha = self.parameters["alpha"] self.alpha = float(self.parameters["alpha"])
self.gamma = self.parameters["gamma"] self.gamma = float(self.parameters["gamma"])
self.umax = 1e7 self.umax = 1e7
with Database() as database: with Database() as database:
......
...@@ -236,6 +236,15 @@ class CalzLeit(CreationModule): ...@@ -236,6 +236,15 @@ class CalzLeit(CreationModule):
def _init_code(self): def _init_code(self):
"""Get the filters from the database""" """Get the filters from the database"""
self.ebvs = {}
self.ebvs['young'] = float(self.parameters["E_BVs_young"])
self.ebvs_old_factor = float(self.parameters["E_BVs_old_factor"])
self.ebvs['old'] = self.ebvs_old_factor * self.ebvs['young']
self.uv_bump_wavelength = float(self.parameters["uv_bump_wavelength"])
self.uv_bump_width = float(self.parameters["uv_bump_width"])
self.uv_bump_amplitude = float(self.parameters["uv_bump_amplitude"])
self.powerlaw_slope = float(self.parameters["powerlaw_slope"])
self.filter_list = [item.strip() for item in self.filter_list = [item.strip() for item in
self.parameters["filters"].split("&")] self.parameters["filters"].split("&")]
# We cannot compute the attenuation until we know the wavelengths. Yet, # We cannot compute the attenuation until we know the wavelengths. Yet,
...@@ -250,24 +259,18 @@ class CalzLeit(CreationModule): ...@@ -250,24 +259,18 @@ class CalzLeit(CreationModule):
sed: pcigale.sed.SED object sed: pcigale.sed.SED object
""" """
ebvs = {}
wavelength = sed.wavelength_grid wavelength = sed.wavelength_grid
ebvs['young'] = float(self.parameters["E_BVs_young"])
ebvs_old_factor = float(self.parameters["E_BVs_old_factor"])
ebvs['old'] = ebvs_old_factor * ebvs['young']
uv_bump_wavelength = float(self.parameters["uv_bump_wavelength"])
uv_bump_width = float(self.parameters["uv_bump_width"])
uv_bump_amplitude = float(self.parameters["uv_bump_amplitude"])
powerlaw_slope = float(self.parameters["powerlaw_slope"])
# Fλ fluxes (only from continuum) in each filter before attenuation. # Fλ fluxes (only from continuum) in each filter before attenuation.
flux_noatt = {filt: sed.compute_fnu(filt) for filt in self.filter_list} flux_noatt = {filt: sed.compute_fnu(filt) for filt in self.filter_list}
# Compute attenuation curve # Compute attenuation curve
if self.sel_attenuation is None: if self.sel_attenuation is None:
self.sel_attenuation = a_vs_ebv(wavelength, uv_bump_wavelength, self.sel_attenuation = a_vs_ebv(wavelength,
uv_bump_width, uv_bump_amplitude, self.uv_bump_wavelength,
powerlaw_slope) self.uv_bump_width,
self.uv_bump_amplitude,
self.powerlaw_slope)
attenuation_total = 0. attenuation_total = 0.
contribs = [contrib for contrib in sed.contribution_names if contribs = [contrib for contrib in sed.contribution_names if
...@@ -275,16 +278,16 @@ class CalzLeit(CreationModule): ...@@ -275,16 +278,16 @@ class CalzLeit(CreationModule):
for contrib in contribs: for contrib in contribs:
age = contrib.split('.')[-1].split('_')[-1] age = contrib.split('.')[-1].split('_')[-1]
luminosity = sed.get_lumin_contribution(contrib) luminosity = sed.get_lumin_contribution(contrib)
attenuated_luminosity = (luminosity * 10 ** attenuated_luminosity = (luminosity * 10. ** (self.ebvs[age] *
(ebvs[age] * self.sel_attenuation / -2.5)) self.sel_attenuation / -2.5))
attenuation_spectrum = attenuated_luminosity - luminosity attenuation_spectrum = attenuated_luminosity - luminosity
# We integrate the amount of luminosity attenuated (-1 because the # We integrate the amount of luminosity attenuated (-1 because the
# spectrum is negative). # spectrum is negative).
attenuation = -1 * np.trapz(attenuation_spectrum, wavelength) attenuation = -1. * np.trapz(attenuation_spectrum, wavelength)
attenuation_total += attenuation attenuation_total += attenuation
sed.add_module(self.name, self.parameters) sed.add_module(self.name, self.parameters)
sed.add_info("attenuation.E_BVs." + contrib, ebvs[age]) sed.add_info("attenuation.E_BVs." + contrib, self.ebvs[age])
sed.add_info("attenuation." + contrib, attenuation, True) sed.add_info("attenuation." + contrib, attenuation, True)
sed.add_contribution("attenuation." + contrib, wavelength, sed.add_contribution("attenuation." + contrib, wavelength,
attenuation_spectrum) attenuation_spectrum)
...@@ -305,11 +308,11 @@ class CalzLeit(CreationModule): ...@@ -305,11 +308,11 @@ class CalzLeit(CreationModule):
sed.add_info("attenuation." + filt, sed.add_info("attenuation." + filt,
-2.5 * np.log10(flux_att[filt] / flux_noatt[filt])) -2.5 * np.log10(flux_att[filt] / flux_noatt[filt]))
sed.add_info('attenuation.ebvs_old_factor', ebvs_old_factor) sed.add_info('attenuation.ebvs_old_factor', self.ebvs_old_factor)
sed.add_info('attenuation.uv_bump_wavelength', uv_bump_wavelength) sed.add_info('attenuation.uv_bump_wavelength', self.uv_bump_wavelength)
sed.add_info('attenuation.uv_bump_width', uv_bump_width) sed.add_info('attenuation.uv_bump_width', self.uv_bump_width)
sed.add_info('attenuation.uv_bump_amplitude', uv_bump_amplitude) sed.add_info('attenuation.uv_bump_amplitude', self.uv_bump_amplitude)
sed.add_info('attenuation.powerlaw_slope', powerlaw_slope) sed.add_info('attenuation.powerlaw_slope', self.powerlaw_slope)
# CreationModule to be returned by get_module # CreationModule to be returned by get_module
Module = CalzLeit Module = CalzLeit
...@@ -152,6 +152,15 @@ class PowerLawAtt(CreationModule): ...@@ -152,6 +152,15 @@ class PowerLawAtt(CreationModule):
]) ])
def _init_code(self): def _init_code(self):
self.av = {}
self.av['young'] = float(self.parameters["Av_young"])
self.av_old_factor = float(self.parameters["Av_old_factor"])
self.av['old'] = self.av_old_factor * self.av['young']
self.uv_bump_wavelength = float(self.parameters["uv_bump_wavelength"])
self.uv_bump_width = float(self.parameters["uv_bump_width"])
self.uv_bump_amplitude = float(self.parameters["uv_bump_amplitude"])
self.powerlaw_slope = float(self.parameters["powerlaw_slope"])
self.filter_list = [item.strip() for item in self.filter_list = [item.strip() for item in
self.parameters["filters"].split("&")] self.parameters["filters"].split("&")]
# We cannot compute the attenuation until we know the wavelengths. Yet, # We cannot compute the attenuation until we know the wavelengths. Yet,
...@@ -166,24 +175,17 @@ class PowerLawAtt(CreationModule): ...@@ -166,24 +175,17 @@ class PowerLawAtt(CreationModule):
sed: pcigale.sed.SED object sed: pcigale.sed.SED object
""" """
av = {}
wavelength = sed.wavelength_grid wavelength = sed.wavelength_grid
av['young'] = float(self.parameters["Av_young"])
av_old_factor = float(self.parameters["Av_old_factor"])
av['old'] = av_old_factor * av['young']
uv_bump_wavelength = float(self.parameters["uv_bump_wavelength"])
uv_bump_width = float(self.parameters["uv_bump_width"])
uv_bump_amplitude = float(self.parameters["uv_bump_amplitude"])
powerlaw_slope = float(self.parameters["powerlaw_slope"])
# Fλ fluxes (only from continuum)) in each filter before attenuation. # Fλ fluxes (only from continuum)) in each filter before attenuation.
flux_noatt = {filt: sed.compute_fnu(filt) for filt in self.filter_list} flux_noatt = {filt: sed.compute_fnu(filt) for filt in self.filter_list}
# Compute attenuation curve # Compute attenuation curve
if self.sel_attenuation is None: if self.sel_attenuation is None:
self.sel_attenuation = alambda_av(wavelength, powerlaw_slope, self.sel_attenuation = alambda_av(wavelength, self.powerlaw_slope,
uv_bump_wavelength, self.uv_bump_wavelength,
uv_bump_width, uv_bump_amplitude) self.uv_bump_width,
self.uv_bump_amplitude)
attenuation_total = 0. attenuation_total = 0.
contribs = [contrib for contrib in sed.contribution_names if contribs = [contrib for contrib in sed.contribution_names if
...@@ -191,8 +193,8 @@ class PowerLawAtt(CreationModule): ...@@ -191,8 +193,8 @@ class PowerLawAtt(CreationModule):
for contrib in contribs: for contrib in contribs:
age = contrib.split('.')[-1].split('_')[-1] age = contrib.split('.')[-1].split('_')[-1]
luminosity = sed.get_lumin_contribution(contrib) luminosity = sed.get_lumin_contribution(contrib)
attenuated_luminosity = (luminosity * 10 ** attenuated_luminosity = (luminosity * 10. **
(av[age] * self.sel_attenuation / -2.5)) (self.av[age]*self.sel_attenuation/-2.5))
attenuation_spectrum = attenuated_luminosity - luminosity attenuation_spectrum = attenuated_luminosity - luminosity
# We integrate the amount of luminosity attenuated (-1 because the # We integrate the amount of luminosity attenuated (-1 because the
# spectrum is negative). # spectrum is negative).
...@@ -200,16 +202,16 @@ class PowerLawAtt(CreationModule): ...@@ -200,16 +202,16 @@ class PowerLawAtt(CreationModule):
attenuation_total += attenuation attenuation_total += attenuation
sed.add_module(self.name, self.parameters) sed.add_module(self.name, self.parameters)
sed.add_info("attenuation.Av." + contrib, av[age]) sed.add_info("attenuation.Av." + contrib, self.av[age])
sed.add_info("attenuation." + contrib, attenuation, True) sed.add_info("attenuation." + contrib, attenuation, True)
sed.add_contribution("attenuation." + contrib, wavelength, sed.add_contribution("attenuation." + contrib, wavelength,
attenuation_spectrum) attenuation_spectrum)
sed.add_info("attenuation.av_old_factor", av_old_factor) sed.add_info("attenuation.av_old_factor", self.av_old_factor)
sed.add_info('attenuation.uv_bump_wavelength', uv_bump_wavelength) sed.add_info('attenuation.uv_bump_wavelength', self.uv_bump_wavelength)
sed.add_info('attenuation.uv_bump_width', uv_bump_width) sed.add_info('attenuation.uv_bump_width', self.uv_bump_width)
sed.add_info("attenuation.uv_bump_amplitude", uv_bump_amplitude) sed.add_info("attenuation.uv_bump_amplitude", self.uv_bump_amplitude)
sed.add_info("attenuation.powerlaw_slope", powerlaw_slope) sed.add_info("attenuation.powerlaw_slope", self.powerlaw_slope)
# Total attenuation # Total attenuation
if 'dust.luminosity' in sed.info: if 'dust.luminosity' in sed.info:
......
...@@ -86,16 +86,18 @@ class Fritz2006(CreationModule): ...@@ -86,16 +86,18 @@ class Fritz2006(CreationModule):
def _init_code(self): def _init_code(self):
"""Get the template set out of the database""" """Get the template set out of the database"""
r_ratio = self.parameters["r_ratio"] self.r_ratio = float(self.parameters["r_ratio"])
tau = self.parameters["tau"] self.tau = float(self.parameters["tau"])
beta = self.parameters["beta"] self.beta = float(self.parameters["beta"])
gamma = self.parameters["gamma"] self.gamma = float(self.parameters["gamma"])
opening_angle = (180. - self.parameters["opening_angle"]) / 2. self.opening_angle = (180. - self.parameters["opening_angle"]) / 2.
psy = self.parameters["psy"] self.psy = float(self.parameters["psy"])
self.fracAGN = float(self.parameters["fracAGN"])
with Database() as base: with Database() as base:
self.fritz2006 = base.get_fritz2006(r_ratio, tau, beta, gamma, self.fritz2006 = base.get_fritz2006(self.r_ratio, self.tau,
opening_angle, psy) self.beta, self.gamma,
self.opening_angle, self.psy)
def process(self, sed): def process(self, sed):
"""Add the IR re-emission contributions """Add the IR re-emission contributions
...@@ -111,20 +113,20 @@ class Fritz2006(CreationModule): ...@@ -111,20 +113,20 @@ class Fritz2006(CreationModule):
sed.add_info('dust.luminosity', 1., True) sed.add_info('dust.luminosity', 1., True)
luminosity = sed.info['dust.luminosity'] luminosity = sed.info['dust.luminosity']
fracAGN = self.parameters["fracAGN"]
sed.add_module(self.name, self.parameters) sed.add_module(self.name, self.parameters)
sed.add_info('agn.r_ratio', self.parameters["r_ratio"]) sed.add_info('agn.r_ratio', self.r_ratio)
sed.add_info('agn.tau', self.parameters["tau"]) sed.add_info('agn.tau', self.tau)
sed.add_info('agn.beta', self.parameters["beta"]) sed.add_info('agn.beta', self.beta)
sed.add_info('agn.gamma', self.parameters["gamma"]) sed.add_info('agn.gamma', self.gamma)
sed.add_info('agn.opening_angle', self.parameters["opening_angle"]) sed.add_info('agn.opening_angle', self.parameters["opening_angle"])
sed.add_info('agn.psy', self.parameters["psy"]) sed.add_info('agn.psy', self.psy)
sed.add_info('agn.fracAGN', self.parameters["fracAGN"]) sed.add_info('agn.fracAGN', self.fracAGN)
# Compute the AGN luminosity # Compute the AGN luminosity
if fracAGN < 1.: if self.fracAGN < 1.:
agn_power = luminosity * (1./(1.-fracAGN) - 1.) agn_power = luminosity * (1./(1.-self.fracAGN) - 1.)
l_agn_therm = agn_power l_agn_therm = agn_power
l_agn_scatt = np.trapz(agn_power * self.fritz2006.lumin_scatt, l_agn_scatt = np.trapz(agn_power * self.fritz2006.lumin_scatt,
x=self.fritz2006.wave) x=self.fritz2006.wave)
......
...@@ -63,15 +63,18 @@ class M2005(CreationModule): ...@@ -63,15 +63,18 @@ class M2005(CreationModule):
def _init_code(self): def _init_code(self):
"""Read the SSP from the database.""" """Read the SSP from the database."""
if self.parameters["imf"] == 0: self.imf = int(self.parameters["imf"])
imf = 'salp' self.metallicity = float(self.parameters["metallicity"])
elif self.parameters["imf"] == 1: self.separation_age = int(self.parameters["separation_age"])
imf = 'krou'
if self.imf == 0:
with Database() as database:
self.ssp = database.get_m2005('salp', self.metallicity)
elif self.imf == 1:
with Database() as database:
self.ssp = database.get_m2005('krou', self.metallicity)
else: else:
raise Exception("IMF #{} unknown".format(self.parameters["imf"])) raise Exception("IMF #{} unknown".format(self.imf))
metallicity = float(self.parameters["metallicity"])
with Database() as database:
self.ssp = database.get_m2005(imf, metallicity)
def process(self, sed): def process(self, sed):
"""Add the convolution of a Maraston 2005 SSP to the SED """Add the convolution of a Maraston 2005 SSP to the SED
...@@ -82,9 +85,6 @@ class M2005(CreationModule): ...@@ -82,9 +85,6 @@ class M2005(CreationModule):
SED object. SED object.
""" """
imf = self.parameters["imf"]
metallicity = float(self.parameters["metallicity"])
separation_age = int(self.parameters["separation_age"])
sfh_time, sfh_sfr = sed.sfh sfh_time, sfh_sfr = sed.sfh
ssp = self.ssp ssp = self.ssp
...@@ -94,20 +94,20 @@ class M2005(CreationModule): ...@@ -94,20 +94,20 @@ class M2005(CreationModule):
# First, we process the young population (age lower than the # First, we process the young population (age lower than the
# separation age.) # separation age.)
young_sfh = np.copy(sfh_sfr) young_sfh = np.copy(sfh_sfr)
young_sfh[sfh_age > separation_age] = 0 young_sfh[sfh_age > self.separation_age] = 0
young_masses, young_spectrum = ssp.convolve(sfh_time, young_sfh) young_masses, young_spectrum = ssp.convolve(sfh_time, young_sfh)
# Then, we process the old population. If the SFH is shorter than the # Then, we process the old population. If the SFH is shorter than the
# separation age then all the arrays will consist only of 0. # separation age then all the arrays will consist only of 0.
old_sfh = np.copy(sfh_sfr) old_sfh = np.copy(sfh_sfr)
old_sfh[sfh_age <= separation_age] = 0 old_sfh[sfh_age <= self.separation_age] = 0
old_masses, old_spectrum = ssp.convolve(sfh_time, old_sfh) old_masses, old_spectrum = ssp.convolve(sfh_time, old_sfh)
sed.add_module(self.name, self.parameters) sed.add_module(self.name, self.parameters)
sed.add_info('stellar.imf', imf) sed.add_info('stellar.imf', self.imf)
sed.add_info('stellar.metallicity', metallicity) sed.add_info('stellar.metallicity', self.metallicity)
sed.add_info('stellar.old_young_separation_age', separation_age) sed.add_info('stellar.old_young_separation_age', self.separation_age)
sed.add_info('stellar.mass_total_old', old_masses[0], True) sed.add_info('stellar.mass_total_old', old_masses[0], True)
sed.add_info('stellar.mass_alive_old', old_masses[1], True) sed.add_info('stellar.mass_alive_old', old_masses[1], True)
......
...@@ -61,13 +61,13 @@ class MBB(CreationModule): ...@@ -61,13 +61,13 @@ class MBB(CreationModule):
def _init_code(self): def _init_code(self):
"""Build the model for a given set of parameters.""" """Build the model for a given set of parameters."""
epsilon = float(self.parameters["epsilon_mbb"]) self.epsilon = float(self.parameters["epsilon_mbb"])
T = float(self.parameters["t_mbb"]) self.T = float(self.parameters["t_mbb"])
beta = float(self.parameters["beta_mbb"]) self.beta = float(self.parameters["beta_mbb"])
self.energy_balance = bool(self.parameters["energy_balance"])
if epsilon < 0: if self.epsilon < 0.:
epsilon = 0. raise Exception("Error, epsilon_mbb must be ≥ 0.")
print("epsilon_mbb must >= 0, we set epsilon_mbb = 0.0")
# We define various constants necessary to compute the model. For # We define various constants necessary to compute the model. For
# consistency, we define the speed of light in nm s¯¹ rather than in # consistency, we define the speed of light in nm s¯¹ rather than in
...@@ -79,8 +79,8 @@ class MBB(CreationModule): ...@@ -79,8 +79,8 @@ class MBB(CreationModule):
conv = c / (self.wave * self.wave) conv = c / (self.wave * self.wave)
self.lumin_mbb = (conv * (1. - np.exp(-(lambda_0 / self.wave) self.lumin_mbb = (conv * (1. - np.exp(-(lambda_0 / self.wave)
** beta)) * (c / self.wave) ** 3. / (np.exp( ** self.beta)) * (c / self.wave) ** 3. / (np.exp(
cst.h * c / (self.wave * cst.k * T)) - 1.)) cst.h * c / (self.wave * cst.k * self.T)) - 1.))
# TODO, save the right normalisation factor to retrieve the dust mass # TODO, save the right normalisation factor to retrieve the dust mass
norm = np.trapz(self.lumin_mbb, x=self.wave) norm = np.trapz(self.lumin_mbb, x=self.wave)
...@@ -99,13 +99,11 @@ class MBB(CreationModule): ...@@ -99,13 +99,11 @@ class MBB(CreationModule):
luminosity = sed.info['dust.luminosity'] luminosity = sed.info['dust.luminosity']
sed.add_module(self.name, self.parameters) sed.add_module(self.name, self.parameters)
sed.add_info("dust.t_mbb", self.parameters["t_mbb"]) sed.add_info("dust.t_mbb", self.T)
sed.add_info("dust.beta_mbb", self.parameters["beta_mbb"]) sed.add_info("dust.beta_mbb", self.beta)
sed.add_info("dust.epsilon_mbb", self.parameters["epsilon_mbb"]) sed.add_info("dust.epsilon_mbb", self.epsilon)
epsilon = float(self.parameters["epsilon_mbb"])
energy_balance = (self.parameters["energy_balance"].lower() == "true")
if energy_balance: if self.energy_balance:
# Since we can have another contribution to L_dust and the modified # Since we can have another contribution to L_dust and the modified
# black body enters into the energy budget, energy balance, we have # black body enters into the energy budget, energy balance, we have
# to save a new negative component for each one, previously # to save a new negative component for each one, previously
...@@ -121,7 +119,7 @@ class MBB(CreationModule): ...@@ -121,7 +119,7 @@ class MBB(CreationModule):
wavelength = sed.wavelength_grid wavelength = sed.wavelength_grid
sed.add_info(item_balance, 1., True) sed.add_info(item_balance, 1., True)
sed.add_contribution(item_balance, wavelength, -lumin * sed.add_contribution(item_balance, wavelength, -lumin *
epsilon) self.epsilon)