dustatt_calzleit.py 12.6 KB
Newer Older
1
# -*- coding: utf-8 -*-
2
3
4
# Copyright (C) 2013 Centre de données Astrophysiques de Marseille
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
# Author: Yannick Roehlly <yannick.roehlly@oamp.fr>
5
6

import numpy as np
7
from collections import OrderedDict
8
9
10
11
12
13
14
15
from . import common
from ...data import Database


def k_calzetti2000(wavelength):
    """Compute the Calzetti et al. (2000) A(λ)/E(B-V)∗

    Given a wavelength grid, this function computes the selective attenuation
16
17
18
    A(λ)/E(B-V)∗ using the formula from Calzetti at al. (2000). This formula
    is given for wavelengths between 120 nm and 2200 nm, but this function
    makes the computation outside.
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

    Parameters
    ----------
    wavelength : array of floats
        Wavelength grid in nm.

    Returns
    -------
    a numpy array of floats

    """
    wavelength = np.array(wavelength)
    result = np.zeros(len(wavelength))

    # Attenuation between 120 nm and 630 nm
34
    mask = (wavelength < 630)
35
36
37
38
39
    result[mask] = 2.659 * (-2.156 + 1.509e3 / wavelength[mask]
                            - 0.198e6 / wavelength[mask] ** 2
                            + 0.011e9 / wavelength[mask] ** 3) + 4.05

    # Attenuation between 630 nm and 2200 nm
40
    mask = (wavelength >= 630)
41
42
43
44
45
46
47
48
49
    result[mask] = 2.659 * (-1.857 + 1.040e3 / wavelength[mask]) + 4.05

    return result


def k_leitherer2002(wavelength):
    """Compute the Leitherer et al. (2002) A(λ)/E(B-V)∗

    Given a wavelength grid, this function computes the selective attenuation
50
51
52
    A(λ)/E(B-V)∗ using the formula from Leitherer at al. (2002). This formula
    is given for wavelengths between 91.2 nm and 180 nm, but this function
    makes the computation outside.
53
54
55
56
57
58
59
60
61
62
63
64

    Parameters
    ----------
    wavelength : array of floats
        Wavelength grid in nm.

    Returns
    -------
    a numpy array of floats

    """
    wavelength = np.array(wavelength)
65
66
67
    result = (5.472 + 0.671e3 / wavelength
              - 9.218e3 / wavelength ** 2
              + 2.620e6 / wavelength ** 3)
68
69
70
71

    return result


72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
def uv_bump(wavelength, central_wave, gamma, ebump):
    """Compute the Lorentzian-like Drude profile.

    Parameters
    ----------
    wavelength : array of floats
        Wavelength grid in nm.
    central_wave : float
        Central wavelength of the bump in nm.
    gamma : float
        Width (FWHM) of the bump in nm.
    ebump : float
        Amplitude of the bump.

    Returns
    -------
    a numpy array of floats

    """
    return (ebump * wavelength ** 2 * gamma ** 2 /
            ((wavelength ** 2 - central_wave ** 2) ** 2
             + wavelength ** 2 * gamma ** 2))


def power_law(wavelength, delta):
    """Power law 'centered' on 550 nm..

    Parameters
    ----------
    wavelength : array of floats
        The wavelength grid in nm.
    delta : float
        The slope of the power law.

    Returns
    -------
    array of floats

    """
    return (wavelength / 550) ** delta


def a_vs_ebv(wavelength, bump_wave, bump_width, bump_ampl, power_slope):
    """Compute the complete attenuation curve A(λ)/E(B-V)*

    The Leitherer et al. (2002) formula is used between 91.2 and 150 nm and
    the Calzetti et al. (2000) formula is used after 150 (we do an
    extrapolation after 2200 nm). When the attenuation becomes negative, it is
    kept to 0. This continuum is multiplied by the power law and then the UV
    bump is added.

    Parameters
    ----------
    wavelength : array of floats
        The wavelength grid (in nm) to compute the attenuation curve on.
    bump_wave : float
        Central wavelength (in nm) of the UV bump.
    bump_width : float
        Width (FWHM, in nm) of the UV bump.
    bump_ampl : float
        Amplitude of the UV bump.
    power_slope : float
        Slope of the power law.

    Returns
    -------
    attenuation : array of floats
        The A(λ)/E(B-V)* attenuation at each wavelength of the grid.

    """
    attenuation = np.zeros(len(wavelength))

    # Leitherer et al.
    mask = (wavelength >= 91.2) & (wavelength < 150)
    attenuation[mask] = k_leitherer2002(wavelength[mask])
    # Calzetti et al.
    mask = (wavelength >= 150)
    attenuation[mask] = k_calzetti2000(wavelength[mask])
    # We set attenuation to 0 where it becomes negative
    mask = (attenuation < 0)
    attenuation[mask] = 0
    # Power law
    attenuation = attenuation * power_law(wavelength, power_slope)
    # UV bump
    attenuation = attenuation + uv_bump(wavelength, bump_wave,
                                        bump_width, bump_ampl)

    return attenuation


162
163
164
165
166
167
168
169
170
class Module(common.SEDCreationModule):
    """Add CCM dust attenuation based on the Calzetti formula

    If a contribution name is given in the parameter list, the attenuation
    will be applied only to the flux of this contribution; else, it will be
    applied to the whole spectrum.

    """

171
172
    parameter_list = OrderedDict([
        ("E_BVs_young", (
173
174
175
176
            "float",
            "E(B-V)*, the colour excess of the stellar continuum light for "
            "the young population.",
            None
177
178
        )),
        ("E_BVs_old_factor", (
179
180
181
182
            "float",
            "Reduction factor for the E(B-V)* of the old population compared "
            "to the young one (<1).",
            None
183
184
        )),
        ("young_contribution_name", (
185
186
187
            "string",
            "Name of the contribution containing the spectrum of the "
            "young population.",
188
            "ssp_young"
189
190
        )),
        ("old_contribution_name", (
191
192
193
194
            "string",
            "Name of the contribution containing the spectrum of the "
            "old population. If it is set to 'None', only one population "
            "is considered.",
195
            "ssp_old"
196
197
        )),
        ("uv_bump_wavelength", (
198
199
200
            "float",
            "Central wavelength of the UV bump in nm.",
            217.5
201
202
        )),
        ("uv_bump_width", (
203
204
205
            "float",
            "Width (FWHM) of the UV bump in nm.",
            None
206
207
        )),
        ("uv_bump_amplitude", (
208
209
210
            "float",
            "Amplitude of the UV bump in nm.",
            None
211
212
        )),
        ("powerlaw_slope", (
213
214
215
            "float",
            "Slope delta of the power law modifying the attenuation curve.",
            None
216
217
        )),
        ("filters", (
218
219
220
221
222
            "string",
            "Filters for which the attenuation will be computed and added to "
            "the SED information dictionary. You can give several filter "
            "names separated by a & (don't use commas).",
            "V_B90 & FUV"
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
        ))
    ])

    out_parameter_list = OrderedDict([
        ("E_BVs_young", "E(B-V)*, the colour excess of the stellar continuum "
                        "light for the young population."),
        ("E_BVs_old", "E(B-V)*, the colour excess of the stellar "
                      "continuum light for the old population."),
        ("attenuation_young", "Amount of luminosity attenuated from the "
                              "young population in W."),
        ("E_BVs_old_factor", "Ratio of the old population E(B-V)* to the "
                             "young one."),
        ("attenuation_old", "Amount of luminosity attenuated from the "
                            "old population in W."),
        ("attenuation", "Total amount of luminosity attenuated in W."),
        ("uv_bump_wavelength", "Central wavelength of UV bump in nm."),
        ("uv_bump_width", "Width of the UV bump in nm."),
        ("uv_bump_amplitude", "Amplitude of the UV bump in nm."),
        ("powerlaw_slope", "Slope of the power law."),
242
        ("FILTER_attenuation", "Attenuation in the FILTER filter.")
243
    ])
244

245
246
    def _init_code(self):
        """Get the filters from the database"""
247
248
        filter_list = [item.strip() for item in
                       self.parameters["filters"].split("&")]
249
250
        self.filters = {}
        base = Database()
251
        for filter_name in filter_list:
252
253
254
255
            self.filters[filter_name] = base.get_filter(filter_name)
        base.close()

    def process(self, sed):
256
257
258
259
260
261
262
263
264
        """Add the CCM dust attenuation to the SED.

        Parameters
        ----------
        sed : pcigale.sed.SED object

        """

        wavelength = sed.wavelength_grid
265
266
267
268
269
270
271
272
273
        ebvs_young = float(self.parameters["E_BVs_young"])
        ebvs_old = float(self.parameters["E_BVs_old_factor"]) * ebvs_young
        young_contrib = self.parameters["young_contribution_name"]
        old_contrib = self.parameters["old_contribution_name"]
        uv_bump_wavelength = float(self.parameters["uv_bump_wavelength"])
        uv_bump_width = float(self.parameters["uv_bump_wavelength"])
        uv_bump_amplitude = float(self.parameters["uv_bump_amplitude"])
        powerlaw_slope = float(self.parameters["powerlaw_slope"])
        filters = self.filters
274

275
        # Fλ fluxes (only from continuum) in each filter before attenuation.
276
        flux_noatt = {}
277
        for filter_name, filter_ in filters.items():
278
            flux_noatt[filter_name] = sed.compute_fnu(
279
                filter_.trans_table,
280
281
                filter_.effective_wavelength,
                add_line_fluxes=False)
282

283
284
285
286
        # Compute attenuation curve
        sel_attenuation = a_vs_ebv(wavelength, uv_bump_wavelength,
                                   uv_bump_width, uv_bump_amplitude,
                                   powerlaw_slope)
287
288
289

        # Young population attenuation
        luminosity = sed.get_lumin_contribution(young_contrib)
290
291
292
        attenuated_luminosity = luminosity * 10 ** (ebvs_young
                                                    * sel_attenuation / -2.5)
        attenuated_luminosity[wavelength < 91.2] = 0
293
294
295
296
297
        attenuation_spectrum = attenuated_luminosity - luminosity
        # We integrate the amount of luminosity attenuated (-1 because the
        # spectrum is negative).
        attenuation_young = -1 * np.trapz(attenuation_spectrum, wavelength)

298
299
300
301
302
        sed.add_module(self.name, self.parameters)
        sed.add_info("E_BVs_young" + self.postfix, ebvs_young)
        sed.add_info("attenuation_young" + self.postfix, attenuation_young)
        sed.add_contribution("attenuation_young" + self.postfix,
                             wavelength, attenuation_spectrum)
303
304
305
306

        # Old population (if any) attenuation
        if old_contrib:
            luminosity = sed.get_lumin_contribution(old_contrib)
307
308
309
310
            attenuated_luminosity = (luminosity *
                                     10 ** (ebvs_old
                                            * sel_attenuation / -2.5))
            attenuated_luminosity[wavelength < 91.2] = 0
311
312
313
            attenuation_spectrum = attenuated_luminosity - luminosity
            attenuation_old = -1 * np.trapz(attenuation_spectrum, wavelength)

314
315
            sed.add_info("E_BVs_old" + self.postfix, ebvs_old)
            sed.add_info("E_BVs_old_factor" + self.postfix,
316
                         self.parameters["E_BVs_old_factor"])
317
318
            sed.add_info("attenuation_old" + self.postfix, attenuation_old)
            sed.add_contribution("attenuation_old" + self.postfix,
319
320
321
322
                                 wavelength, attenuation_spectrum)
        else:
            attenuation_old = 0

323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
        # Attenuation of spectral lines
        if young_contrib in sed.lines:
            line_attenuation = a_vs_ebv(sed.lines[young_contrib][0],
                                        uv_bump_wavelength, uv_bump_width,
                                        uv_bump_amplitude, powerlaw_slope)
            sed.lines[young_contrib][1] *= 10 ** (ebvs_young *
                                                  line_attenuation / -2.5)
        if old_contrib in sed.lines:
            line_attenuation = a_vs_ebv(sed.lines[old_contrib][0],
                                        uv_bump_wavelength, uv_bump_width,
                                        uv_bump_amplitude, powerlaw_slope)
            sed.lines[old_contrib][1] *= 10 ** (ebvs_old *
                                                line_attenuation / -2.5)

        # Total attenuation (we don't take into account the energy attenuated
        # in the spectral lines)
339
        sed.add_info("attenuation" + self.postfix,
340
341
                     attenuation_young + attenuation_old)

342
        # Fλ fluxes (only from continuum) in each filter after attenuation.
343
        flux_att = {}
344
345
346
        for filter_name, filter_ in filters.items():
            flux_att[filter_name] = sed.compute_fnu(
                filter_.trans_table,
347
348
                filter_.effective_wavelength,
                add_line_fluxes=False)
349
350

        # Attenuation in each filter
351
        for filter_name in filters:
352
            sed.add_info(filter_name + "_attenuation" + self.postfix,
353
354
                         -2.5 * np.log(flux_att[filter_name] /
                                       flux_noatt[filter_name]))