dustatt_calzleit.py 11.1 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 SedModule
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 between 91.2 nm 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.
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 > 91.2) & (wavelength < 150)
157
158
159
160
161
162
163
164
    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
165
    attenuation *= power_law(wavelength, power_slope)
166

167
168
169
170
171
172
173
174
175
176
    # As the powerlaw slope changes E(B-V), we correct this so that the curve
    # always has the same E(B-V) as the starburst curve. This ensures that the
    # E(B-V) requested by the user is the actual E(B-V) of the curve.
    wl_BV = np.array([440., 550.])
    EBV_calz = ((k_calzetti2000(wl_BV) * power_law(wl_BV, 0.)) +
                uv_bump(wl_BV, bump_wave, bump_width, bump_ampl))
    EBV = ((k_calzetti2000(wl_BV) * power_law(wl_BV, power_slope)) +
           uv_bump(wl_BV, bump_wave, bump_width, bump_ampl))
    attenuation *= (EBV_calz[1]-EBV_calz[0]) / (EBV[1]-EBV[0])

177
178
179
180
181
    # UV bump. It is added after the renormalization as the bump strength
    # should correspond to the requested E(B-V) and should therefore not be
    # changed by the renormalization.
    attenuation += uv_bump(wavelength, bump_wave, bump_width, bump_ampl)

182
183
184
    return attenuation


185
class CalzLeit(SedModule):
Yannick Roehlly's avatar
Yannick Roehlly committed
186
187
    """Calzetti + Leitherer attenuation module

188
189
190
191
    This module computes the dust attenuation using the formulae from
    Calzetti et al. (2000) and Leitherer et al. (2002). Note that both the
    stars and the gas are attenuated with the same curve as opposed to Calzetti
    et al. (2000) where the gas is attenuated with a Milky Way curve.
192

Yannick Roehlly's avatar
Yannick Roehlly committed
193
194
    The attenuation can be computed on the whole spectrum or on a specific
    contribution and is added to the SED as a negative contribution.
195
196
197

    """

198
    parameter_list = OrderedDict([
199
        ("E_BVs_young", (
200
            "cigale_list(minvalue=0.)",
201
            "E(B-V)*, the colour excess of the stellar continuum light for "
202
            "the young population.",
203
            0.3
204
205
        )),
        ("E_BVs_old_factor", (
206
            "cigale_list(minvalue=0., maxvalue=1.)",
207
208
            "Reduction factor for the E(B-V)* of the old population compared "
            "to the young one (<1).",
209
            1.0
210
211
        )),
        ("uv_bump_wavelength", (
212
            "cigale_list(minvalue=0.)",
213
214
            "Central wavelength of the UV bump in nm.",
            217.5
215
216
        )),
        ("uv_bump_width", (
217
            "cigale_list()",
218
            "Width (FWHM) of the UV bump in nm.",
219
            35.
220
221
        )),
        ("uv_bump_amplitude", (
222
            "cigale_list(minvalue=0.)",
223
224
            "Amplitude of the UV bump. For the Milky Way: 3.",
            0.
225
226
        )),
        ("powerlaw_slope", (
227
            "cigale_list()",
228
            "Slope delta of the power law modifying the attenuation curve.",
229
            0.
230
231
        )),
        ("filters", (
232
            "string()",
233
234
235
            "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).",
236
            "B_B90 & V_B90 & FUV"
237
238
239
        ))
    ])

240
241
    def _init_code(self):
        """Get the filters from the database"""
242
243
244
245
246
247
248
249
250
        self.ebvs = {}
        self.ebvs['young'] = float(self.parameters["E_BVs_young"])
        self.ebvs_old_factor = float(self.parameters["E_BVs_old_factor"])
        self.ebvs['old'] = self.ebvs_old_factor * self.ebvs['young']
        self.uv_bump_wavelength = float(self.parameters["uv_bump_wavelength"])
        self.uv_bump_width = float(self.parameters["uv_bump_width"])
        self.uv_bump_amplitude = float(self.parameters["uv_bump_amplitude"])
        self.powerlaw_slope = float(self.parameters["powerlaw_slope"])

251
252
        self.filter_list = [item.strip() for item in
                            self.parameters["filters"].split("&")]
253
254
255
        # We cannot compute the attenuation until we know the wavelengths. Yet,
        # we reserve the object.
        self.sel_attenuation = None
256
257

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

        Parameters
        ----------
262
        sed: pcigale.sed.SED object
263
264
265
266

        """
        wavelength = sed.wavelength_grid

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

270
        # Compute attenuation curve
271
        if self.sel_attenuation is None:
272
273
274
275
276
            self.sel_attenuation = a_vs_ebv(wavelength,
                                            self.uv_bump_wavelength,
                                            self.uv_bump_width,
                                            self.uv_bump_amplitude,
                                            self.powerlaw_slope)
277

278
        attenuation_total = 0.
279
280
281
        contribs = [contrib for contrib in sed.contribution_names if
                    'absorption' not in contrib]
        for contrib in contribs:
282
283
            age = contrib.split('.')[-1].split('_')[-1]
            luminosity = sed.get_lumin_contribution(contrib)
284
285
            attenuated_luminosity = (luminosity * 10. ** (self.ebvs[age] *
                                     self.sel_attenuation / -2.5))
286
            attenuation_spectrum = attenuated_luminosity - luminosity
287
288
            # We integrate the amount of luminosity attenuated (-1 because the
            # spectrum is negative).
289
            attenuation = -1. * np.trapz(attenuation_spectrum, wavelength)
290
291
292
            attenuation_total += attenuation

            sed.add_module(self.name, self.parameters)
293
            sed.add_info("attenuation.E_BVs." + contrib, self.ebvs[age])
294
            sed.add_info("attenuation." + contrib, attenuation, True)
295
296
297
298
            sed.add_contribution("attenuation." + contrib, wavelength,
                                 attenuation_spectrum)

        # Total attenuation
299
300
301
302
303
304
        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)
305

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

        # Attenuation in each filter
310
311
312
        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
313

314
315
316
317
318
        sed.add_info('attenuation.ebvs_old_factor', self.ebvs_old_factor)
        sed.add_info('attenuation.uv_bump_wavelength', self.uv_bump_wavelength)
        sed.add_info('attenuation.uv_bump_width', self.uv_bump_width)
        sed.add_info('attenuation.uv_bump_amplitude', self.uv_bump_amplitude)
        sed.add_info('attenuation.powerlaw_slope', self.powerlaw_slope)
319

320

321
# SedModule to be returned by get_module
Yannick Roehlly's avatar
Yannick Roehlly committed
322
Module = CalzLeit