redshifting.py 7.64 KB
Newer Older
Yannick Roehlly's avatar
Yannick Roehlly committed
1
# -*- coding: utf-8 -*-
2
# Copyright (C) 2014 Yannick Roehlly, Médéric Boquien, Denis Burgarella
Yannick Roehlly's avatar
Yannick Roehlly committed
3
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
4
# Author: Yannick Roehlly, Médéric Boquien, Denis Burgarella
Yannick Roehlly's avatar
Yannick Roehlly committed
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

"""
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.

"""

21 22
from collections import OrderedDict

23 24
import numpy as np
from scipy.constants import parsec
25
from scipy.special import factorial
26
from ..utils.cosmology import age, luminosity_distance
27

28
from . import SedModule
29

Médéric Boquien's avatar
Médéric Boquien committed
30

31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
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
50
    n_transitions_max = 31
51 52 53 54 55 56 57 58
    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))
59
        z_n[n, :] = (wavelength / lambda_n[n]) - 1.
Médéric Boquien's avatar
Médéric Boquien committed
60 61

    # From Table 1 in Meiksin (2006), only n >= 3 are relevant.
62
    # fact has a length equal to n_transitions_low.
63 64 65 66 67 68 69
    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:
Médéric Boquien's avatar
Médéric Boquien committed
70
        tau_a = 0.00211 * np.power(1. + redshift,  3.7)
71
        tau_n[2, :] = 0.00211 * np.power(1. + z_n[2, :], 3.7)
72
    elif redshift > 4:
Médéric Boquien's avatar
Médéric Boquien committed
73
        tau_a = 0.00058 * np.power(1. + redshift,  4.5)
74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92
        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.))))

93 94 95
    for n in range(2, n_transitions_max):
        w = np.where(z_n[n, :] >= redshift)
        tau_n[n, w] = 0.
Médéric Boquien's avatar
Médéric Boquien committed
96

97 98
    z_l = wavelength / lambda_limit - 1.
    w = np.where(z_l < redshift)
99 100

    tau_l_igm = np.zeros_like(wavelength)
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121
    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)
    tau_l_lls[w] = n0 * ((term1 - term2) * term3 - term4)

122
    tau_taun = np.sum(tau_n[2:n_transitions_max, :], axis=0)
Médéric Boquien's avatar
Médéric Boquien committed
123

124 125
    lambda_min_igm = (1+redshift)*70.
    w = np.where(wavelength < lambda_min_igm)
126 127

    weight = np.ones_like(wavelength)
128 129
    weight[w] = np.power(wavelength[w]/lambda_min_igm, 2.)
    # Another weight using erf function can be used.
130
    # However, you would need to add: from scipy.special import erf
Médéric Boquien's avatar
Médéric Boquien committed
131
    # weight[w] = 0.5*(1.+erf(0.05*(wavelength[w]-lambda_min_igm)))
132

133
    igm_transmission = np.exp(-tau_taun-tau_l_igm-tau_l_lls) * weight
134 135

    return igm_transmission
Yannick Roehlly's avatar
Yannick Roehlly committed
136 137


138
class Redshifting(SedModule):
Yannick Roehlly's avatar
Yannick Roehlly committed
139 140 141 142 143 144 145
    """Redshift a SED

    This module redshift a rest-frame SED. If the SED is already redshifted, an
    exception is raised.

    """

146
    parameter_list = OrderedDict([
Yannick Roehlly's avatar
Yannick Roehlly committed
147
        ("redshift", (
148
            "cigale_list(minvalue=0.)",
149 150
            "Redshift of the objects. Leave empty to use the redshifts from the"
            " input file.",
Yannick Roehlly's avatar
Yannick Roehlly committed
151 152 153 154
            None
        ))
    ])

155 156 157 158
    def _init_code(self):
        """Compute the age of the Universe at a given redshift
        """
        self.redshift = float(self.parameters["redshift"])
159 160 161 162

        # Raise an error when applying a negative redshift. This module is
        # not for blue-shifting.
        if self.redshift < 0.:
163 164
            raise Exception(f"The redshift provided is negative "
                            f"({self.redshift}).")
165

166 167
        self.universe_age = age(self.redshift)
        self.luminosity_distance = luminosity_distance(self.redshift)
168 169 170
        # 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.
171
        self.igm_attenuation = {}
172

Yannick Roehlly's avatar
Yannick Roehlly committed
173 174 175 176 177
    def process(self, sed):
        """Redshift the SED

        Parameters
        ----------
178
        sed: pcigale.sed.SED object
Yannick Roehlly's avatar
Yannick Roehlly committed
179 180

        """
181 182
        redshift = self.redshift

Yannick Roehlly's avatar
Yannick Roehlly committed
183
        # If the SED is already redshifted, raise an error.
Médéric Boquien's avatar
Médéric Boquien committed
184 185
        if ('universe.redshift' in sed.info and
            sed.info['universe.redshift'] > 0.):
186 187
            raise Exception(f"The SED is already redshifted (z="
                            f"{sed.info['universe.redshift']}).")
Yannick Roehlly's avatar
Yannick Roehlly committed
188

189 190 191
        if redshift > 0.:
            # We redshift directly the SED wavelength grid
            sed.wavelength_grid *= 1. + redshift
Yannick Roehlly's avatar
Yannick Roehlly committed
192

193
            # We modify each luminosity contribution to keep energy constant
194 195
            sed.luminosities *= 1. / (1. + redshift)
            sed.luminosity *= 1. / (1. + redshift)
Yannick Roehlly's avatar
Yannick Roehlly committed
196

Médéric Boquien's avatar
Médéric Boquien committed
197
        sed.add_info("universe.redshift", redshift)
198
        sed.add_info("universe.luminosity_distance", self.luminosity_distance)
199
        sed.add_info("universe.age", self.universe_age)
200 201 202 203 204 205 206

        # We identify the right grid from the length of the wavelength array.
        # It is not completely foolproof but it is good enough. We need to do
        # that in case two models do not have the same wavelength sampling.
        # This is the case for instance if some but not all models have an AGN
        # fraction of 0.
        key = sed.wavelength_grid.size
207
        if key not in self.igm_attenuation:
208 209
            self.igm_attenuation[key] = igm_transmission(sed.wavelength_grid,
                                                         redshift) - 1.
210
        sed.add_contribution('igm', sed.wavelength_grid,
211
                             self.igm_attenuation[key] * sed.luminosity)
Yannick Roehlly's avatar
Yannick Roehlly committed
212 213
        sed.add_module(self.name, self.parameters)

214

215
# SedModule to be returned by get_module
Yannick Roehlly's avatar
Yannick Roehlly committed
216
Module = Redshifting