Commit 08c5c6b4 authored by Yannick Roehlly's avatar Yannick Roehlly

Correct the redshifting

The redshifting also lower the fluxes (that's different from the
cosmological dimming).
Also, simplify some methods.
parent 38747f50
......@@ -118,7 +118,7 @@ class SED(object):
else:
return self._luminosities.sum(0)
def lambda_fnu(self, redshift=0, redshift_spectrum=False):
def lambda_fnu(self, redshift=0, apply_redshift=False):
"""
Return the (redshifted if asked) total Fν flux density vs wavelength
spectrum of the SED.
......@@ -126,8 +126,8 @@ class SED(object):
Parameters
----------
redshift : float, default = 0
If 0 (the default), the flux at 10 pc is computed.
redshift_spectrum : boolean, default = None
If 0 (the default), the flux at 10 pc is computed.
apply_redshift : boolean, default = None
If true, the spectrum will be redshifted before computing the
flux. The default is False because we generally use a specific
module to apply the redshift.
......@@ -138,19 +138,18 @@ class SED(object):
The wavelength is in nm and the Fν is in mJy.
"""
# Fλ flux density in W/m²/nm
f_lambda = utils.luminosity_to_flux(self.luminosity, redshift)
wavelength = self.wavelength_grid
luminosity = self.luminosity
# The 1.e29 factor is to convert from W/m²/Hz to mJy
f_nu = (f_lambda * self.wavelength_grid
* 1.e-9 * self.wavelength_grid / c
* 1.e29)
if apply_redshift:
wavelength, luminosity = utils.redshift_spectrum(
wavelength, luminosity, redshift)
if redshift_spectrum:
wavelength = utils.redshift_wavelength(self.wavelength_grid,
redshift)
else:
wavelength = np.copy(self.wavelength_grid)
# Fλ flux density in W/m²/nm
f_lambda = utils.luminosity_to_flux(luminosity, redshift)
# Fν flux density in Jy
f_nu = utils.lambda_flambda_to_fnu(wavelength, f_lambda)
return wavelength, f_nu
......@@ -283,7 +282,7 @@ class SED(object):
return self.luminosities[idx]
def compute_fnu(self, transmission, lambda_eff,
redshift=0, redshift_spectrum=False):
redshift=0, apply_redshift=False):
"""
Compute the Fν flux density corresponding the filter which
transmission is given.
......@@ -296,7 +295,7 @@ class SED(object):
Fλ = luminosity_to_flux( integ( LλT(λ)dλ ) / integ( T(λ)dλ ) )
Fλ is in W/m²/nm. At redshift 0, the flux is computed at 10 pc. Then,
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λ
......@@ -316,9 +315,9 @@ class SED(object):
Effective wavelength of the filter in nm.
redshift : float
The redshift of the galaxy. If 0, the flux is computed at 10 pc.
The redshift of the galaxy. If 0, the flux is computed at 10 pc.
redshift_spectrum : boolean
apply_redshift : boolean
If true, the spectrum will be redshifted before computing the
flux. The default is False because we generally use a specific
module to apply the redshift.
......@@ -338,13 +337,11 @@ class SED(object):
f_nu = -99.
else:
if redshift_spectrum:
wavelength = utils.redshift_wavelength(self.wavelength_grid,
redshift)
else:
wavelength = np.copy(self.wavelength_grid)
wavelength = self.wavelength_grid
l_lambda = self.luminosity
if apply_redshift:
wavelength, l_lambda = utils.redshift_lambda_l_lambda(
(wavelength, l_lambda), redshift)
# We regrid both spectrum and filter to the best wavelength grid
# to avoid interpolating a high wavelength density curve to a low
......
# -*- coding: utf-8 -*-
"""
Copyright (C) 2012 Centre de données Astrophysiques de Marseille
Copyright (C) 2012, 2013 Centre de données Astrophysiques de Marseille
Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
@author: Yannick Roehlly <yannick.roehlly@oamp.fr>
......@@ -148,124 +148,104 @@ def luminosity_to_flux(luminosity, redshift=0):
return luminosity / (4 * pi * np.square(dist))
def redshift_wavelength(wavelength, redshift):
"""Redshift a wavelength grid
Parameters
----------
wavelength : array of floats
Wavelength vector.
redshift : float
Redshift.
Returns
-------
redshifted_wavelength : array of floats
Redshifted wavelength grid.
"""
if redshift < 0:
return wavelength / (1.0 - redshift)
else:
return wavelength * (1.0 + redshift)
def lambda_flambda_to_lambda_fnu(spectrum):
def lambda_flambda_to_fnu(wavelength, flambda):
"""
Convert a Fλ vs λ spectrum to Fν vs λ
Parameters
----------
spectrum : array of floats
spectrum[0] must contain the wavelength in nm and spectrum[1] must
contain the Fλ flux in erg/cm^2/s/nm.
wavelength : list-like of floats
The wavelengths in nm.
flambda : list-like of floats
Fλ flux density in W/m²/nm (or Lλ luminosity density in W/nm).
Returns
-------
lambda_fnu : array of floats
lambda_fnu[0] contains the wavelength in nm and lambda_fnu[1] contains
the Fν flux in Jansky
fnu : array of floats
The Fν flux density in mJy (or the Lν luminosity density in
1.e-29 W/Hz).
"""
wavelength, flambda = spectrum
# Factor 1e+23 is to switch from erg/s/cm^2/Hz to Jy
wavelength = np.array(wavelength, dtype=float)
flambda = np.array(flambda, dtype=float)
# Factor 1e+29 is to switch from W/m²/Hz to mJy
# Factor 1e-9 is to switch from nm to m (only one because the other nm
# wavelength goes with the Fλ in ergs/s/cm^2/nm).
fnu = 1e+23 * 1e-9 * flambda * wavelength * wavelength / c
# wavelength goes with the Fλ in W/m²/nm).
fnu = 1e+29 * 1e-9 * flambda * wavelength * wavelength / c
return np.vstack((wavelength, fnu))
return fnu
def lambda_fnu_to_lambda_flambda(spectrum):
def lambda_fnu_to_flambda(wavelength, fnu):
"""
Convert a Fν vs λ spectrum to Fλ vs λ
Parameters
----------
spectrum : array of floats
spectrum[0] must contain the wavelength in nm and spectrum[1] must
contain the Fν flux in Jansky
wavelength : list-like of floats
The wavelengths in nm.
fnu : list-like of floats
The Fν flux density in mJy (of the Lν luminosity density in
1.e-29 W/Hz).
Returns
-------
lambda_flambda : array of floats
lambda_flambda[0] contains the wavelength in nm and lambda_flambda[1]
contains the Fλ flux in erg/cm^2/s/nm.
flambda : array of floats
Fλ flux density in W/m²/nm (or Lλ luminosity density in W/nm).
"""
wavelength, fnu = spectrum
# Factor 1e-23 is to switch from Jy to erg/s/cm^2/Hz
# Factor 1e+9 is to switch from m to nm
flambda = 1e-23 * 1e+9 * fnu / (wavelength * wavelength) * c
wavelength = np.array(wavelength, dtype=float)
fnu = np.array(fnu, dtype=float)
return np.vstack((wavelength, flambda))
# Factor 1e-29 is to switch from Jy to W/m²/Hz
# Factor 1e+9 is to switch from m to nm
flambda = 1e-29 * 1e+9 * fnu / (wavelength * wavelength) * c
return flambda
def redshift_spectrum(spectrum, redshift, dimming=False, is_fnu=False):
"""
Redshit a spectrum, optionally adding cosmological dimming
FIXME: Is this usefull?
def redshift_spectrum(wavelength, flux, redshift, is_fnu=False):
"""Redshit a spectrum
Parameters
----------
spectrum : array of floats
spectrum[0] must contain the wavelength in nm and spectrum[1] must
contain the flux. The default is to have Fλ in erg/cm^2/s/nm if is_fnu
is set to true, the Fν in Jansky is expected (it's only important when
dimming).
dimming : boolean
If set to true, the cosmological dimming is applied to the fluxes.
wavelength : array like of floats
The wavelength in nm.
flux : array like of floats
The flux or luminosity density.
redshift : float
The redshift.
is_fnu : boolean
If set to true, the flux are Fν fluxes, else they are assumed to be Fλ.
If false (default) the flux is a Fλ density in W/m²/nm (or a Lλ
luminosity density in W/nm). If true, the flux is a Fν density in mJy
(or a Lν luminosity density in 1.e-29 W/Hz).
Results
-------
spectrum : array of floats
The redshifted spectrum with the same kind of fluxes as the input.
wavelength, flux : tuple of numpy arrays of floats
The redshifted spectrum with the same kind of flux (or luminosity)
density as the input.
"""
wavelength = np.array(wavelength, dtype=float)
flux = np.array(flux, dtype=float)
redshift = float(redshift)
wavelength = redshift_wavelength(spectrum[0])
flux = np.copy(spectrum[1])
if redshift < 0:
redshift_factor = 1. / (1. - redshift)
else:
redshift_factor = 1. + redshift
if dimming:
# If the flux is Fnu, we must switch to Flambda to compute the
# dimming.
if is_fnu:
flux = lambda_fnu_to_lambda_flambda(spectrum)[1]
if is_fnu:
# Switch to Fλ
flux = lambda_fnu_to_flambda(wavelength, flux)
# Now flux is Flambda, we can apply cosmological dim.
if redshift < 0:
flux = flux * (1.0 - redshift)
else:
flux = flux / (1.0 + redshift)
wavelength *= redshift_factor
flux /= redshift_factor
# If the initial flux was Fnu, convert it back from Flambda
if is_fnu:
flux = lambda_flambda_to_lambda_fnu(
np.vstack((wavelength, flux)))[:, 1]
if is_fnu:
# Switch back to Fλ
flux = lambda_flambda_to_fnu(wavelength, flux)
return np.vstack((wavelength, flux))
return wavelength, flux
......@@ -305,10 +305,8 @@ class Module(common.AnalysisModule):
ax = figure.add_subplot(111)
plot_x, plot_y = best_sed_lambda_fnu
plot_mask = (
(plot_x >= utils.redshift_wavelength(PLOT_L_MIN,
obs_redshift)) &
(plot_x <= utils.redshift_wavelength(PLOT_L_MAX,
obs_redshift))
(plot_x >= PLOT_L_MIN * (1 + obs_redshift)) &
(plot_x <= PLOT_L_MAX * (1 + obs_redshift))
)
ax.loglog(plot_x[plot_mask],
best_norm_factor * plot_y[plot_mask],
......
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