dustatt_calzleit.py 9.87 KB
Newer Older
1
# -*- coding: utf-8 -*-
2
# Copyright (C) 2013 Centre de données Astrophysiques de Marseille
3
# Copyright (C) 2014 Laboratoire d'Astrophysique de Marseille
4
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
5
# Author: Yannick Roehlly, Denis Burgarella
6

Yannick Roehlly's avatar
Yannick Roehlly committed
7
8
9
10
11
12
13
14
15
"""
Calzetti et al. (2000) and Leitherer et al. (2002) attenuation module
=====================================================================

This module implements the Calzetti et al. (2000) and  Leitherer et al. (2002)
attenuation formulae, adding an UV-bump and a power law.

"""

16
from collections import OrderedDict
17

18
import numpy as np
19

20
from . import CreationModule
21
22
23
24
25
26


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

    Given a wavelength grid, this function computes the selective attenuation
27
28
29
    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.
30
31
32

    Parameters
    ----------
33
    wavelength: array of floats
34
35
36
37
38
39
40
41
42
43
44
        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
45
    mask = (wavelength < 630)
Médéric Boquien's avatar
Médéric Boquien committed
46
47
48
    result[mask] = 2.659 * (-2.156 + 1.509e3 / wavelength[mask] -
                            0.198e6 / wavelength[mask] ** 2 +
                            0.011e9 / wavelength[mask] ** 3) + 4.05
49
50

    # Attenuation between 630 nm and 2200 nm
51
    mask = (wavelength >= 630)
52
53
54
55
56
57
58
59
60
    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
61
62
63
    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.
64
65
66

    Parameters
    ----------
67
    wavelength: array of floats
68
69
70
71
72
73
74
75
        Wavelength grid in nm.

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

    """
    wavelength = np.array(wavelength)
Médéric Boquien's avatar
Médéric Boquien committed
76
77
78
    result = (5.472 + 0.671e3 / wavelength -
              9.218e3 / wavelength ** 2 +
              2.620e6 / wavelength ** 3)
79
80
81
82

    return result


83
84
85
86
87
def uv_bump(wavelength, central_wave, gamma, ebump):
    """Compute the Lorentzian-like Drude profile.

    Parameters
    ----------
88
    wavelength: array of floats
89
        Wavelength grid in nm.
90
    central_wave: float
91
        Central wavelength of the bump in nm.
92
    gamma: float
93
        Width (FWHM) of the bump in nm.
94
    ebump: float
95
96
97
98
99
100
101
102
        Amplitude of the bump.

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

    """
    return (ebump * wavelength ** 2 * gamma ** 2 /
Médéric Boquien's avatar
Médéric Boquien committed
103
104
            ((wavelength ** 2 - central_wave ** 2) ** 2 +
             wavelength ** 2 * gamma ** 2))
105
106
107
108
109
110
111


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

    Parameters
    ----------
112
    wavelength: array of floats
113
        The wavelength grid in nm.
114
    delta: float
115
116
117
118
119
120
121
122
123
124
125
126
127
        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)*

128
129
130
131
132
    The Leitherer et al. (2002) formula is used bellow 150 nm (even if it is
    defined only after 91.2 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.
133
134
135

    Parameters
    ----------
136
    wavelength: array of floats
137
        The wavelength grid (in nm) to compute the attenuation curve on.
138
    bump_wave: float
139
        Central wavelength (in nm) of the UV bump.
140
    bump_width: float
141
        Width (FWHM, in nm) of the UV bump.
142
    bump_ampl: float
143
        Amplitude of the UV bump.
144
    power_slope: float
145
146
147
148
        Slope of the power law.

    Returns
    -------
149
    attenuation: array of floats
150
151
152
153
154
155
        The A(λ)/E(B-V)* attenuation at each wavelength of the grid.

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

    # Leitherer et al.
156
    mask = (wavelength < 150)
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
    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


Yannick Roehlly's avatar
Yannick Roehlly committed
173
class CalzLeit(CreationModule):
Yannick Roehlly's avatar
Yannick Roehlly committed
174
175
    """Calzetti + Leitherer attenuation module

176
    This module computes the dust attenuation using the
Yannick Roehlly's avatar
Yannick Roehlly committed
177
    formulae from Calzetti et al. (2000) and Leitherer et al. (2002).
178

Yannick Roehlly's avatar
Yannick Roehlly committed
179
180
    The attenuation can be computed on the whole spectrum or on a specific
    contribution and is added to the SED as a negative contribution.
181
182
183

    """

184
    parameter_list = OrderedDict([
185
        ("E_BVs_young", (
186
187
188
            "float",
            "E(B-V)*, the colour excess of the stellar continuum light for "
            "the young population.",
189
            0.3
190
191
        )),
        ("E_BVs_old_factor", (
192
193
194
            "float",
            "Reduction factor for the E(B-V)* of the old population compared "
            "to the young one (<1).",
195
            0.44
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
            "float",
            "Width (FWHM) of the UV bump in nm.",
205
            35.
206
207
        )),
        ("uv_bump_amplitude", (
208
            "float",
209
210
            "Amplitude of the UV bump. For the Milky Way: 3.",
            0.
211
212
        )),
        ("powerlaw_slope", (
213
214
            "float",
            "Slope delta of the power law modifying the attenuation curve.",
215
            0.
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
    def _init_code(self):
        """Get the filters from the database"""
228
229
        self.filter_list = [item.strip() for item in
                            self.parameters["filters"].split("&")]
230
231
232
        # We cannot compute the attenuation until we know the wavelengths. Yet,
        # we reserve the object.
        self.sel_attenuation = None
233
234

    def process(self, sed):
235
236
237
238
        """Add the CCM dust attenuation to the SED.

        Parameters
        ----------
239
        sed: pcigale.sed.SED object
240
241

        """
242
        ebvs = {}
243
        wavelength = sed.wavelength_grid
244
245
        ebvs['young'] = float(self.parameters["E_BVs_young"])
        ebvs['old'] = float(self.parameters["E_BVs_old_factor"]) * ebvs['young']
246
        uv_bump_wavelength = float(self.parameters["uv_bump_wavelength"])
247
        uv_bump_width = float(self.parameters["uv_bump_width"])
248
249
        uv_bump_amplitude = float(self.parameters["uv_bump_amplitude"])
        powerlaw_slope = float(self.parameters["powerlaw_slope"])
250

251
        # Fλ fluxes (only from continuum) in each filter before attenuation.
Médéric Boquien's avatar
Médéric Boquien committed
252
        flux_noatt = {filt: sed.compute_fnu(filt) for filt in self.filter_list}
253

254
        # Compute attenuation curve
255
256
        if self.sel_attenuation is None:
            self.sel_attenuation = a_vs_ebv(wavelength, uv_bump_wavelength,
Médéric Boquien's avatar
Médéric Boquien committed
257
258
                                            uv_bump_width, uv_bump_amplitude,
                                            powerlaw_slope)
259

260
261
262
263
264
        attenuation_total = 0.
        for contrib in list(sed.contribution_names):
            age = contrib.split('.')[-1].split('_')[-1]
            luminosity = sed.get_lumin_contribution(contrib)
            attenuated_luminosity = (luminosity * 10 **
Médéric Boquien's avatar
Médéric Boquien committed
265
                                     (ebvs[age] * self.sel_attenuation / -2.5))
266
            attenuation_spectrum = attenuated_luminosity - luminosity
267
268
269
270
271
272
273
            # We integrate the amount of luminosity attenuated (-1 because the
            # spectrum is negative).
            attenuation = -1 * np.trapz(attenuation_spectrum, wavelength)
            attenuation_total += attenuation

            sed.add_module(self.name, self.parameters)
            sed.add_info("attenuation.E_BVs." + contrib, ebvs[age])
274
            sed.add_info("attenuation." + contrib, attenuation, True)
275
276
277
278
            sed.add_contribution("attenuation." + contrib, wavelength,
                                 attenuation_spectrum)

        # Total attenuation
279
280
281
282
283
284
        if 'dust.luminosity' in sed.info:
            sed.add_info("dust.luminosity",
                         sed.info["dust.luminosity"]+attenuation_total, True,
                         True)
        else:
            sed.add_info("dust.luminosity", attenuation_total, True)
285

286
        # Fλ fluxes (only from continuum) in each filter after attenuation.
Médéric Boquien's avatar
Médéric Boquien committed
287
        flux_att = {filt: sed.compute_fnu(filt) for filt in self.filter_list}
288
289

        # Attenuation in each filter
290
291
292
        for filt in self.filter_list:
            sed.add_info("attenuation." + filt,
                         -2.5 * np.log10(flux_att[filt] / flux_noatt[filt]))
Yannick Roehlly's avatar
Yannick Roehlly committed
293

294
295
296
297
298
299
300
        sed.add_info('attenuation.ebvs_main', ebvs['old'])
        sed.add_info('attenuation.ebvs_young', ebvs['young'])
        sed.add_info('attenuation.uv_bump_wavelength', uv_bump_wavelength)
        sed.add_info('attenuation.uv_bump_width', uv_bump_width)
        sed.add_info('attenuation.uv_bump_amplitude', uv_bump_amplitude)
        sed.add_info('attenuation.powerlaw_slope', powerlaw_slope)

Yannick Roehlly's avatar
Yannick Roehlly committed
301
302
# CreationModule to be returned by get_module
Module = CalzLeit