redshifting.py 7.74 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, Guang Yang
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
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

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
49
    n_transitions_max = 31
50
51
52
53
54
55
56
57
    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))
58
        z_n[n, :] = (wavelength / lambda_n[n]) - 1.
Médéric Boquien's avatar
Médéric Boquien committed
59
60

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

92
    for n in range(2, n_transitions_max):
Guang's avatar
Guang committed
93
        # If z_n>=redshift or z_n<0, the photon cannot be absorbed by Lyman n->1
Guang's avatar
Guang committed
94
        w = (z_n[n, :] >= redshift) | (z_n[n, :] < 0)
95
        tau_n[n, w] = 0.
Médéric Boquien's avatar
Médéric Boquien committed
96

97
    z_l = wavelength / lambda_limit - 1.
Guang's avatar
Guang committed
98
    w = 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)

Guang's avatar
Guang committed
122
    # Reset for short wavelength (z_l<0)
Guang's avatar
Guang committed
123
    w = z_l<0
Guang's avatar
Guang committed
124
125
126
127
128
129
    # Get the normalization factor at z_l=0
    tau_norm_l_igm = np.interp(0, z_l, tau_l_igm)
    tau_norm_l_lls = np.interp(0, z_l, tau_l_lls)
    # Calculate tau_l_igm & tau_l_lls, assuming cross section ~ lambda^2.75 (O'Meara et al. 2013)
    tau_l_igm[w] = tau_norm_l_igm * (z_l[w]+1)**2.75
    tau_l_lls[w] = tau_norm_l_lls * (z_l[w]+1)**2.75
130

Guang's avatar
Guang committed
131
    tau_taun = np.sum(tau_n[2:n_transitions_max, :], axis=0)
132

Guang's avatar
Guang committed
133
    igm_transmission = np.exp(-tau_taun-tau_l_igm-tau_l_lls)
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