redshifting.py 7.77 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
21

"""
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
22
23
24
25
26

import numpy as np
from scipy.constants import parsec
from scipy.misc import factorial

Yannick Roehlly's avatar
Yannick Roehlly committed
27
from ..creation_modules import CreationModule
28
from pcigale.sed.cosmology import cosmology
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
139
140
141
142
143
144
145
146
147
148


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",
149
            "Redshift to apply to the galaxy. Leave empty to use the redshifts"
Médéric Boquien's avatar
Médéric Boquien committed
150
            " from the input file.",
Yannick Roehlly's avatar
Yannick Roehlly committed
151
152
153
154
            None
        ))
    ])

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

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

        Parameters
        ----------
176
        sed: pcigale.sed.SED object
Yannick Roehlly's avatar
Yannick Roehlly committed
177
178

        """
179
180
        redshift = self.redshift

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

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

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

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

200
        sed.add_info("redshift", redshift)
201
        sed.add_info("universe.luminosity_distance", self.luminosity_distance)
202
        sed.add_info("universe.age", self.universe_age)
203
204
205
206
207
208
209

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

# CreationModule to be returned by get_module
Module = Redshifting