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

Fν was computed by calculating Fλ and then converting to Fν, which led to...

Fν was computed by calculating Fλ and then converting to Fν, which led to typical differences in fluxes of typically 1-2% and a bit more for a handful of pathological filters. Now Fν is computed directly and a bit faster. At the same time rather than computing an effective wavelength (which was more a mean wavelength), we compute the pivot wavelength, which is more appropriate when converting between Fν and Fλ.
parent 4d78de48
...@@ -17,6 +17,7 @@ ...@@ -17,6 +17,7 @@
- To accommodate the analysis of the observations by blocks, the `save_pdf` option has been eliminated. To plot PDF one needs to set `save_chi2` to True and then run `pcigale-plots pdf`. (Médéric Boquien) - To accommodate the analysis of the observations by blocks, the `save_pdf` option has been eliminated. To plot PDF one needs to set `save_chi2` to True and then run `pcigale-plots pdf`. (Médéric Boquien)
- In order to capture rapid evolutionary phases, we assume that in a given period of 1 Myr, 10 small episodes of star formation occurred every 0.1 Myr, rather than one episode every 1 Myr. - In order to capture rapid evolutionary phases, we assume that in a given period of 1 Myr, 10 small episodes of star formation occurred every 0.1 Myr, rather than one episode every 1 Myr.
- When computing the attenuation curve, the bump is now added so that its relative strength does not depend on δ. (Médéric Boquien, issue reported by Samir Salim) - When computing the attenuation curve, the bump is now added so that its relative strength does not depend on δ. (Médéric Boquien, issue reported by Samir Salim)
- Fν was computed by calculating Fλ and then converting to Fν, which led to typical differences in fluxes of typically 1-2% and a bit more for a handful of pathological filters. Now Fν is computed directly and a bit faster. (Médéric Boquien, issue reported by Yannick Roehlly and Wouter Dobbels)
### Fixed ### Fixed
- Corrected a typo that prevented `restframe\_parameters` from being listed among the available modules. (Médéric Boquien) - Corrected a typo that prevented `restframe\_parameters` from being listed among the available modules. (Médéric Boquien)
......
...@@ -165,13 +165,13 @@ def build_filters(base): ...@@ -165,13 +165,13 @@ def build_filters(base):
new_filter = Filter(filter_name, filter_description, filter_table) new_filter = Filter(filter_name, filter_description, filter_table)
# We normalise the filter and compute the effective wavelength. # We normalise the filter and compute the pivot wavelength. If the
# If the filter is a pseudo-filter used to compute line fluxes, it # filter is a pseudo-filter used to compute line fluxes, it should not
# should not be normalised. # be normalised.
if not filter_name.startswith('PSEUDO'): if not filter_name.startswith('PSEUDO'):
new_filter.normalise() new_filter.normalise()
else: else:
new_filter.effective_wavelength = np.mean( new_filter.pivot_wavelength = np.mean(
filter_table[0][filter_table[1] > 0] filter_table[0][filter_table[1] > 0]
) )
filters.append(new_filter) filters.append(new_filter)
...@@ -211,13 +211,13 @@ def build_filters_gazpar(base): ...@@ -211,13 +211,13 @@ def build_filters_gazpar(base):
new_filter = Filter(filter_name, filter_desc, filter_table) new_filter = Filter(filter_name, filter_desc, filter_table)
# We normalise the filter and compute the effective wavelength. # We normalise the filter and compute the pivot wavelength. If the
# If the filter is a pseudo-filter used to compute line fluxes, it # filter is a pseudo-filter used to compute line fluxes, it should not
# should not be normalised. # be normalised.
if not filter_name.startswith('PSEUDO'): if not filter_name.startswith('PSEUDO'):
new_filter.normalise() new_filter.normalise()
else: else:
new_filter.effective_wavelength = np.mean( new_filter.pivot_wavelength = np.mean(
filter_table[0][filter_table[1] > 0] filter_table[0][filter_table[1] > 0]
) )
filters.append(new_filter) filters.append(new_filter)
......
...@@ -63,13 +63,13 @@ class _Filter(BASE): ...@@ -63,13 +63,13 @@ class _Filter(BASE):
name = Column(String, primary_key=True) name = Column(String, primary_key=True)
description = Column(String) description = Column(String)
trans_table = Column(PickleType) trans_table = Column(PickleType)
effective_wavelength = Column(Float) pivot_wavelength = Column(Float)
def __init__(self, f): def __init__(self, f):
self.name = f.name self.name = f.name
self.description = f.description self.description = f.description
self.trans_table = f.trans_table self.trans_table = f.trans_table
self.effective_wavelength = f.effective_wavelength self.pivot_wavelength = f.pivot_wavelength
class _M2005(BASE): class _M2005(BASE):
...@@ -1062,7 +1062,7 @@ class Database(object): ...@@ -1062,7 +1062,7 @@ class Database(object):
first()) first())
if result: if result:
return Filter(result.name, result.description, result.trans_table, return Filter(result.name, result.description, result.trans_table,
result.effective_wavelength) result.pivot_wavelength)
else: else:
raise DatabaseLookupError( raise DatabaseLookupError(
"The filter <{0}> is not in the database".format(name)) "The filter <{0}> is not in the database".format(name))
...@@ -1081,7 +1081,7 @@ class Database(object): ...@@ -1081,7 +1081,7 @@ class Database(object):
"""Generator to parse the filter database.""" """Generator to parse the filter database."""
for filt in self.session.query(_Filter): for filt in self.session.query(_Filter):
yield Filter(filt.name, filt.description, filt.trans_table, yield Filter(filt.name, filt.description, filt.trans_table,
filt.effective_wavelength) filt.pivot_wavelength)
def parse_m2005(self): def parse_m2005(self):
"""Generator to parse the Maraston 2005 SSP database.""" """Generator to parse the Maraston 2005 SSP database."""
......
...@@ -4,6 +4,7 @@ ...@@ -4,6 +4,7 @@
# Author: Yannick Roehlly # Author: Yannick Roehlly
import numpy as np import numpy as np
from scipy.constants import c
class Filter(object): class Filter(object):
...@@ -11,9 +12,9 @@ class Filter(object): ...@@ -11,9 +12,9 @@ class Filter(object):
""" """
def __init__(self, name, description=None, trans_table=None, def __init__(self, name, description=None, trans_table=None,
effective_wavelength=None): pivot_wavelength=None):
"""Create a new filter. If the transmission type, the description """Create a new filter. If the transmission type, the description
the transmission table or the effective wavelength are not specified, the transmission table or the pivot wavelength are not specified,
their value is set to None. their value is set to None.
Parameters Parameters
...@@ -25,14 +26,14 @@ class Filter(object): ...@@ -25,14 +26,14 @@ class Filter(object):
trans_table: array trans_table: array
trans_table[0] is the wavelength in nm, trans_table[0] is the wavelength in nm,
trans_table[1] is the transmission) trans_table[1] is the transmission)
effective_wavelength: float pivot_wavelength: float
Effective wavelength of the filter Pivot wavelength of the filter
""" """
self.name = name self.name = name
self.description = description self.description = description
self.trans_table = trans_table self.trans_table = trans_table
self.effective_wavelength = effective_wavelength self.pivot_wavelength = pivot_wavelength
def __str__(self): def __str__(self):
"""Pretty print the filter information """Pretty print the filter information
...@@ -40,18 +41,24 @@ class Filter(object): ...@@ -40,18 +41,24 @@ class Filter(object):
result = "" result = ""
result += ("Filter name: %s\n" % self.name) result += ("Filter name: %s\n" % self.name)
result += ("Description: %s\n" % self.description) result += ("Description: %s\n" % self.description)
result += ("Effective wavelength: %s nm\n" % result += ("Pivot wavelength: %s nm\n" %
self.effective_wavelength) self.pivot_wavelength)
return result return result
def normalise(self): def normalise(self):
""" """
Normalise the transmission table to 1 and compute the effective Compute the pivot wavelength of the filter and normalise the filter
wavelength of the filter. to compute the flux in Fν (mJy) in cigale.
""" """
self.trans_table[1] = self.trans_table[1] / ( self.pivot_wavelength = np.sqrt(np.trapz(self.trans_table[1],
np.trapz(self.trans_table[1], self.trans_table[0])) self.trans_table[0]) /
np.trapz(self.trans_table[1] /
self.trans_table[0] ** 2,
self.trans_table[0]))
self.effective_wavelength = np.trapz(self.trans_table[1] * # The factor 10²⁰ is so that we get the fluxes directly in mJy when we
self.trans_table[0], # integrate with the wavelength in units of nm and the spectrum in
self.trans_table[0]) # units of W/m²/nm.
self.trans_table[1] = 1e20 * self.trans_table[1] / (
c * np.trapz(self.trans_table[1] / self.trans_table[0]**2,
self.trans_table[0]))
...@@ -243,20 +243,10 @@ class SED(object): ...@@ -243,20 +243,10 @@ class SED(object):
""" """
Compute the Fν flux density in a given filter Compute the Fν flux density in a given filter
As the SED stores the Lλ luminosity density, we first compute the Fλ The filters are stored in the database in such a way that after
flux density. Fλ is the integration of the Lλ luminosity multiplied by integration and conversion from luminosity to flux we directly get the
the filter transmission, normalised to this transmission and corrected latter in units of mJy. If the SED spectrum does not cover all the
by the luminosity distance of the source. filter response table, NaN is returned.
Fλ is in W/m²/nm. At redshift 0, the flux is computed at 10 pc. Then,
to compute Fν, we make the approximation:
Fν = λeff / c . λeff . Fλ
Fν is computed in W/m²/Hz and then converted to mJy.
If the SED spectrum does not cover all the filter response table,
NaN is returned.
Parameters Parameters
---------- ----------
...@@ -273,7 +263,7 @@ class SED(object): ...@@ -273,7 +263,7 @@ class SED(object):
wavelength = self.wavelength_grid wavelength = self.wavelength_grid
# First we try to fetch the filter's wavelength, transmission and # First we try to fetch the filter's wavelength, transmission and
# effective wavelength from the cache. The two keys are the size of the # pivot wavelength from the cache. The two keys are the size of the
# spectrum wavelength grid and the name of the filter. The first key is # spectrum wavelength grid and the name of the filter. The first key is
# necessary because different spectra may have different sampling. To # necessary because different spectra may have different sampling. To
# avoid having the resample the filter every time on the optimal grid # avoid having the resample the filter every time on the optimal grid
...@@ -288,12 +278,12 @@ class SED(object): ...@@ -288,12 +278,12 @@ class SED(object):
dist = 10. * parsec dist = 10. * parsec
if key in self.cache_filters: if key in self.cache_filters:
wavelength_r, transmission_r, lambda_eff = self.cache_filters[key] wavelength_r, transmission_r, lambda_piv = self.cache_filters[key]
else: else:
with Database() as db: with Database() as db:
filter_ = db.get_filter(filter_name) filter_ = db.get_filter(filter_name)
trans_table = filter_.trans_table trans_table = filter_.trans_table
lambda_eff = filter_.effective_wavelength lambda_piv = filter_.pivot_wavelength
lambda_min = filter_.trans_table[0][0] lambda_min = filter_.trans_table[0][0]
lambda_max = filter_.trans_table[0][-1] lambda_max = filter_.trans_table[0][-1]
...@@ -313,17 +303,20 @@ class SED(object): ...@@ -313,17 +303,20 @@ class SED(object):
trans_table[1]) trans_table[1])
self.cache_filters[key] = (wavelength_r, transmission_r, self.cache_filters[key] = (wavelength_r, transmission_r,
lambda_eff) lambda_piv)
l_lambda_r = interp(wavelength_r, wavelength, self.luminosity) l_lambda_r = interp(wavelength_r, wavelength, self.luminosity)
f_lambda = utils.luminosity_to_flux( # We compute directly Fν from ∫T×Fλ×dλ/∫T×c/λ²×dλ. The filter bandpass
# in the database is already normalised so that we do not need to
# compute the denominator (it is a constant that does not depend on the
# spectrum) and the normalisation is such that the results we obtain
# are directly in mJy.
f_nu = utils.luminosity_to_flux(
utils.flux_trapz(transmission_r * l_lambda_r, wavelength_r, key), utils.flux_trapz(transmission_r * l_lambda_r, wavelength_r, key),
dist) dist)
# Return Fν in mJy. The 1e-9 factor is because λ is in nm and 1e29 for return f_nu
# convert from W/m²/Hz to mJy.
return 1e-9 / c * 1e29 * lambda_eff * lambda_eff * f_lambda
def to_votable(self, filename, mass=1.): def to_votable(self, filename, mass=1.):
""" """
......
...@@ -27,14 +27,14 @@ def list_filters(): ...@@ -27,14 +27,14 @@ def list_filters():
name = Column(data=[filters[f].name for f in filters], name='Name') name = Column(data=[filters[f].name for f in filters], name='Name')
description = Column(data=[filters[f].description for f in filters], description = Column(data=[filters[f].description for f in filters],
name='Description') name='Description')
wl = Column(data=[filters[f].effective_wavelength for f in filters], wl = Column(data=[filters[f].pivot_wavelength for f in filters],
name='Effective Wavelength', unit=u.nm, format='%d') name='Pivot Wavelength', unit=u.nm, format='%d')
samples = Column(data=[filters[f].trans_table[0].size for f in filters], samples = Column(data=[filters[f].trans_table[0].size for f in filters],
name="Points") name="Points")
t = Table() t = Table()
t.add_columns([name, description, wl, samples]) t.add_columns([name, description, wl, samples])
t.sort(['Effective Wavelength']) t.sort(['Pivot Wavelength'])
t.pprint(max_lines=-1, max_width=-1) t.pprint(max_lines=-1, max_width=-1)
...@@ -67,13 +67,13 @@ def add_filters(fnames): ...@@ -67,13 +67,13 @@ def add_filters(fnames):
new_filter = Filter(filter_name, filter_description, filter_table) new_filter = Filter(filter_name, filter_description, filter_table)
# We normalise the filter and compute the effective wavelength. # We normalise the filter and compute the pivot wavelength. If the
# If the filter is a pseudo-filter used to compute line fluxes, it # filter is a pseudo-filter used to compute line fluxes, it should
# should not be normalised. # not be normalised.
if not filter_name.startswith('PSEUDO'): if not filter_name.startswith('PSEUDO'):
new_filter.normalise() new_filter.normalise()
else: else:
new_filter.effective_wavelength = np.mean( new_filter.pivot_wavelength = np.mean(
filter_table[0][filter_table[1] > 0] filter_table[0][filter_table[1] > 0]
) )
......
...@@ -144,7 +144,7 @@ def _sed_worker(obs, mod, filters, sed_type, nologo): ...@@ -144,7 +144,7 @@ def _sed_worker(obs, mod, filters, sed_type, nologo):
sed = Table.read("out/{}_best_model.fits".format(obs['id'])) sed = Table.read("out/{}_best_model.fits".format(obs['id']))
filters_wl = np.array([filt.effective_wavelength filters_wl = np.array([filt.pivot_wavelength
for filt in filters.values()]) for filt in filters.values()])
wavelength_spec = sed['wavelength'] wavelength_spec = sed['wavelength']
obs_fluxes = np.array([obs[filt] for filt in filters.keys()]) obs_fluxes = np.array([obs[filt] for filt in filters.keys()])
......
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