Commit 3bb7e018 authored by Médéric Boquien's avatar Médéric Boquien
Browse files

Revert "Correction of an error in the IGM attenuation in redshifting.py"

This reverts commit c160b8f8.
parent c160b8f8
# -*- coding: utf-8 -*-
# Copyright (C) 2014 Yannick Roehlly, Médéric Boquien, Denis Burgarella
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
# Author: Yannick Roehlly, Médéric Boquien, Denis Burgarella
"""
Redshifting module
==================
This module implements the redshifting of a SED. The SED must be rest-frame
or the module will raise en exception when processing it.
Note that this module, contrary to the other SED creation modules, actually
changes the individual luminosity contributions as it redshifts everyone.
Also note that doing this, this module does not rely on the SED object
interface but on its inner implementations. That means that if the SED object
is changed, this module may need to be adapted.
"""
from collections import OrderedDict
import numpy as np
from scipy.constants import parsec
from scipy.misc import factorial
from ..creation_modules import CreationModule
from pcigale.sed.cosmology import cosmology
def igm_transmission(wavelength, redshift):
"""Intergalactic transmission (Meiksin, 2006)
Compute the intergalactic transmission as described in Meiksin, 2006.
Parameters
----------
wavelength: array like of floats
The wavelength(s) in nm.
redshift: float
The redshift. Must be strictly positive.
Returns
-------
igm_transmission: numpy array of floats
The intergalactic transmission at each input wavelength.
"""
n_transitions_low = 10
n_transitions_max = 31
gamma = 0.2788 # Gamma(0.5,1) i.e., Gamma(2-beta,1) with beta = 1.5
n0 = 0.25
lambda_limit = 91.2 # Lyman limit in nm
lambda_n = np.empty(n_transitions_max)
z_n = np.empty((n_transitions_max, len(wavelength)))
for n in range(2, n_transitions_max):
lambda_n[n] = lambda_limit / (1. - 1. / float(n*n))
z_n[n, :] = (wavelength / lambda_n[n]) - 1.
# From Table 1 in Meiksin (2006), only n >= 3 are relevant.
# fact has a length equal to n_transitions_low.
fact = np.array([1., 1., 1., 0.348, 0.179, 0.109, 0.0722, 0.0508, 0.0373,
0.0283])
# First, tau_alpha is the mean Lyman alpha transmitted flux,
# Here n = 2 => tau_2 = tau_alpha
tau_n = np.zeros((n_transitions_max, len(wavelength)))
if redshift <= 4:
tau_a = 0.00211 * np.power(1. + redshift, 3.7)
tau_n[2, :] = 0.00211 * np.power(1. + z_n[2, :], 3.7)
elif redshift > 4:
tau_a = 0.00058 * np.power(1. + redshift, 4.5)
tau_n[2, :] = 0.00058 * np.power(1. + z_n[2, :], 4.5)
# Then, tau_n is the mean optical depth value for transitions
# n = 3 - 9 -> 1
for n in range(3, n_transitions_max):
if n <= 5:
w = np.where(z_n[n, :] < 3)
tau_n[n, w] = (tau_a * fact[n] *
np.power(0.25 * (1. + z_n[n, w]), (1. / 3.)))
w = np.where(z_n[n, :] >= 3)
tau_n[n, w] = (tau_a * fact[n] *
np.power(0.25 * (1. + z_n[n, w]), (1. / 6.)))
elif 5 < n <= 9:
tau_n[n, :] = (tau_a * fact[n] *
np.power(0.25 * (1. + z_n[n, :]), (1. / 3.)))
else:
tau_n[n, :] = (tau_n[9, :] * 720. /
(float(n) * (float(n*n - 1.))))
for n in range(2, n_transitions_max):
w = np.where(z_n[n, :] >= redshift)
tau_n[n, w] = 0.
z_l = wavelength / lambda_limit - 1.
tau_l_igm = np.zeros_like(wavelength)
w = np.where(z_l < redshift)
tau_l_igm[w] = (0.805 * np.power(1. + z_l[w], 3) *
(1. / (1. + z_l[w]) - 1. / (1. + redshift)))
term1 = gamma - np.exp(-1.)
n = np.arange(n_transitions_low - 1)
term2 = np.sum(np.power(-1., n) / (factorial(n) * (2*n - 1)))
term3 = ((1.+redshift) * np.power(wavelength[w]/lambda_limit, 1.5) -
np.power(wavelength[w]/lambda_limit, 2.5))
term4 = np.sum(np.array(
[((2.*np.power(-1., n) / (factorial(n) * ((6*n - 5)*(2*n - 1)))) *
((1.+redshift) ** (2.5-(3 * n)) *
(wavelength[w]/lambda_limit) ** (3*n) -
(wavelength[w]/lambda_limit) ** 2.5))
for n in np.arange(1, n_transitions_low)]), axis=0)
tau_l_lls = np.zeros_like(wavelength)
w = np.where(z_l < redshift)
tau_l_lls[w] = n0 * ((term1 - term2) * term3 - term4)
tau_taun = np.sum(tau_n[2:n_transitions_max, :], axis=0.)
lambda_min_igm = (1+redshift)*70.
weight = np.ones_like(wavelength)
w = np.where(wavelength < lambda_min_igm)
weight[w] = np.power(wavelength[w]/lambda_min_igm, 2.)
# Another weight using erf function can be used.
# However, you would need to add: from scipy.special import erf
# weight[w] = 0.5*(1.+erf(0.05*(wavelength[w]-lambda_min_igm)))
tau = tau_taun + tau_l_igm + tau_l_lls
igm_transmission = np.exp(-tau) * weight
return igm_transmission
class Redshifting(CreationModule):
"""Redshift a SED
This module redshift a rest-frame SED. If the SED is already redshifted, an
exception is raised.
"""
parameter_list = OrderedDict([
("redshift", (
"float",
"Redshift to apply to the galaxy. Leave empty to use the redshifts"
"from the input file.",
None
))
])
def _init_code(self):
"""Compute the age of the Universe at a given redshift
"""
self.redshift = float(self.parameters["redshift"])
self.universe_age = cosmology.age(self.redshift).value * 1000.
if self.redshift == 0.:
self.luminosity_distance = 10. * parsec
else:
self.luminosity_distance = (
cosmology.luminosity_distance(self.redshift).value * 1e6 *
parsec)
# We do not define the values of the IGM attenuation component yet.
# This is because we need the wavelength grid for that first. This
# will be assigned on the first call.
self.igm_attenuation = None
def process(self, sed):
"""Redshift the SED
Parameters
----------
sed: pcigale.sed.SED object
"""
# If the SED is already redshifted, raise an error.
if 'redshift' in sed.info.keys() and sed.info['redshift'] > 0.:
raise Exception("The SED is already redshifted <z={}>."
.format(sed.info['redshift']))
# Raise an error when applying a negative redshift. This module is
# not for blue-shifting.
if self.redshift < 0:
raise Exception("The redshift provided is negative <{}>."
.format(self.redshift))
# We redshift directly the SED wavelength grid
sed.wavelength_grid *= 1. + self.redshift
# We modify each luminosity contribution to keep energy constant
sed.luminosities /= 1. + self.redshift
sed.add_info("redshift", self.redshift)
sed.add_info("universe.luminosity_distance", self.luminosity_distance)
sed.add_info("universe.age", self.universe_age)
if self.igm_attenuation is None:
self.igm_attenuation = igm_transmission(sed.wavelength_grid,
self.redshift) - 1.
sed.add_contribution('igm', sed.wavelength_grid,
self.igm_attenuation * sed.luminosity)
sed.add_module(self.name, self.parameters)
# CreationModule to be returned by get_module
Module = Redshifting
Supports Markdown
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