Commit 63455e66 authored by Médéric Boquien's avatar Médéric Boquien
Browse files

Change the way redshifts are handled. Now all redshifts are specified in...

Change the way redshifts are handled. Now all redshifts are specified in pcigale.ini. To compute the χ², we simply take the models corresponding to the closest χ².
parent f5ebff9c
......@@ -13,93 +13,6 @@ from ...sed.cosmology import cosmology
from ...creation_modules import get_module as get_creation_module
def gen_compute_fluxes_at_redshift(sed, filters, redshifting_module,
igm_module):
""""Generate function to compute the fluxes of a SED at a given redshift
Given a SED, a list of filters and a redshift module name, this generator
returns a function computing the fluxes of the SED in all the filters at
a given redshift. If the SED is older than the Universe at the given
redshift, the fluxes returned by this function will be all -99.
If the redshift module name is None, the SED is not redshifted by the
returned function. This means that it will always return the same fluxes,
whatever its redshift parameter. This can be used for computing
photometric redshifts.
Parameters
----------
sed : pcigale.sed
The pcigale SED object.
filters : list of pcigale.data.filters
List of pcigale filters objects.
redshifting_module : picgale.creation_modules.Module
A pcigale SED creation module to redshift the SED. It must accept a
redshift parameter (or be None).
igm_module : picgale.creation_modules.Module
A pcigale SED creation module to add the IGM attenuation to the SED.
Is not used if the redshifting module is None.
Return
------
gen_fluxes : function
Function computing the fluxes of the SED in all filters at a given
redshift.
"""
# The returned function is memoized.
cache = {}
def gen_fluxes(redshift):
"""Compute the flux of the SED in various filters.
Parameters
----------
redshift : float
Returns
-------
array fo floats
"""
# If the function is generated without a redshift module, it always
# computes the fluxes at the SED redshift (the SED may be already
# redshifted).
if not redshifting_module:
redshift = sed.redshift
if redshift not in cache:
# Age of the Universe at redshift.
# astropy 0.3 cosmology functions return quantities
try:
age_at_redshift = cosmology.age(redshift).value * 1000
except AttributeError:
age_at_redshift = cosmology.age(redshift) * 1000
if "age" in sed.info.keys() and sed.info["age"] > age_at_redshift:
cache[redshift] = -99 * np.ones(len(filters))
else:
if redshifting_module:
red_sed = deepcopy(sed)
redshifting_module.parameters["redshift"] = redshift
redshifting_module.process(red_sed)
if igm_module:
igm_module.process(red_sed)
else:
red_sed = sed
cache[redshift] = np.array(
[red_sed.compute_fnu(filt_.trans_table,
filt_.effective_wavelength)
for filt_ in filters])
return cache[redshift]
return gen_fluxes
def gen_pdf(values, probabilities, grid):
"""Generate a probability density function
......
......@@ -42,16 +42,17 @@ class IgmAttenuation(CreationModule):
sed : pcigale.sed.SED object
"""
if sed.redshift > 0:
redshift = sed.info['redshift']
if redshift > 0:
# To memoize the IGM attenuation computation (with depends on the
# wavelength grid and the redshift), we have to make the
# wavelength grid hashable. We do this by converting to a string.
wave_redshift = (sed.wavelength_grid.tostring(), sed.redshift)
wave_redshift = (sed.wavelength_grid.tostring(), redshift)
if wave_redshift not in self.igm_att:
self.igm_att[wave_redshift] = (
igm_transmission_meiksin(sed.wavelength_grid,
sed.redshift) - 1)
redshift) - 1)
igm_effect = self.igm_att[wave_redshift] * sed .luminosity
......
......@@ -49,9 +49,9 @@ class Redshifting(CreationModule):
redshift = float(self.parameters["redshift"])
# If the SED is already redshifted, raise an error.
if sed.redshift > 0:
if 'redshift' in sed.info.keys() > 0:
raise Exception("The SED is already redshifted <z={}>."
.format(sed.redshift))
.format(sed.info['redshift']))
# Raise an error when applying a negative redshift. This module is
# not for blue-shifting.
......@@ -65,7 +65,7 @@ class Redshifting(CreationModule):
# We modify each luminosity contribution to keep energy constant
sed.luminosities /= 1. + redshift
sed.redshift = redshift
sed.add_info("redshift", redshift)
sed.add_module(self.name, self.parameters)
# CreationModule to be returned by get_module
......
......@@ -65,7 +65,6 @@ class SED(object):
"""
self.sfh = sfh
self.redshift = 0.
self.modules = []
self.wavelength_grid = None
self.contribution_names = []
......@@ -143,7 +142,7 @@ class SED(object):
"""
# Fλ flux density in W/m²/nm
f_lambda = utils.luminosity_to_flux(self.luminosity, self.redshift)
f_lambda = utils.luminosity_to_flux(self.luminosity, self.info['redshift'])
# Fν flux density in mJy
f_nu = utils.lambda_flambda_to_fnu(self.wavelength_grid, f_lambda)
......@@ -388,7 +387,7 @@ class SED(object):
f_lambda = utils.luminosity_to_flux(
(np.trapz(transmission_r * l_lambda_r, wavelength_r) /
np.trapz(transmission_r, wavelength_r)),
self.redshift
self.info['redshift']
)
# Add the Fλ fluxes from the spectral lines.
......
......@@ -77,7 +77,7 @@ def save_sed_to_vo(sed, filename, norm=1.):
# If there is a stellar population then the norm factor is the stellar
# mass.
votable.infos.append(Info(name="Galaxy mass in Msun", value=norm))
votable.infos.append(Info(name="Redshift", value=sed.redshift))
votable.infos.append(Info(name="Redshift", value=sed.info['redshift']))
for name, value in sed.info.items():
if name in sed.mass_proportional_info:
votable.infos.append(Info(name=name, value=norm * value))
......
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