redshifting.py 7.78 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
23
import numpy as np
from scipy.constants import parsec
from scipy.misc import factorial
24
from astropy.cosmology import WMAP7 as cosmology
25

Yannick Roehlly's avatar
Yannick Roehlly committed
26
from ..creation_modules import CreationModule
27

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

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

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

91
92
93
    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
94

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

    tau_l_igm = np.zeros_like(wavelength)
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
    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)

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

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

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

131
    igm_transmission = np.exp(-tau_taun-tau_l_igm-tau_l_lls) * weight
132
133

    return igm_transmission
Yannick Roehlly's avatar
Yannick Roehlly committed
134
135
136
137
138
139
140
141
142
143


class Redshifting(CreationModule):
    """Redshift a SED

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

    """

144
    parameter_list = dict([
Yannick Roehlly's avatar
Yannick Roehlly committed
145
146
        ("redshift", (
            "float",
147
            "Redshift to apply to the galaxy. Leave empty to use the redshifts"
Médéric Boquien's avatar
Médéric Boquien committed
148
            " from the input file.",
Yannick Roehlly's avatar
Yannick Roehlly committed
149
150
151
152
            None
        ))
    ])

153
154
155
156
157
    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.
158
159
160
161
162
163
        if self.redshift == 0.:
            self.luminosity_distance = 10. * parsec
        else:
            self.luminosity_distance = (
                cosmology.luminosity_distance(self.redshift).value * 1e6 *
                parsec)
164
165
166
        # 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.
167
        self.igm_attenuation = {}
168

Yannick Roehlly's avatar
Yannick Roehlly committed
169
170
171
172
173
    def process(self, sed):
        """Redshift the SED

        Parameters
        ----------
174
        sed: pcigale.sed.SED object
Yannick Roehlly's avatar
Yannick Roehlly committed
175
176

        """
177
178
        redshift = self.redshift

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

        # Raise an error when applying a negative redshift. This module is
        # not for blue-shifting.
187
        if redshift < 0:
188
            raise Exception("The redshift provided is negative <{}>."
189
                            .format(redshift))
Yannick Roehlly's avatar
Yannick Roehlly committed
190

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

195
196
            # We modify each luminosity contribution to keep energy constant
            sed.luminosities /= 1. + redshift
197
            sed.luminosity /= 1. + redshift
Yannick Roehlly's avatar
Yannick Roehlly committed
198

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

        # 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
209
        if key not in self.igm_attenuation:
210
211
            self.igm_attenuation[key] = igm_transmission(sed.wavelength_grid,
                                                         redshift) - 1.
212
        sed.add_contribution('igm', sed.wavelength_grid,
213
                             self.igm_attenuation[key] * sed.luminosity)
Yannick Roehlly's avatar
Yannick Roehlly committed
214
215
216
217
        sed.add_module(self.name, self.parameters)

# CreationModule to be returned by get_module
Module = Redshifting