dustatt_calzleit.py 12.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# -*- coding: utf-8 -*-
"""
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>

"""


import numpy as np
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
20
21
22
    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.
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

    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
38
    mask = (wavelength < 630)
39
40
41
42
43
    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
44
    mask = (wavelength >= 630)
45
46
47
48
49
50
51
52
53
    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
54
55
56
    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.
57
58
59
60
61
62
63
64
65
66
67
68

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

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

    """
    wavelength = np.array(wavelength)
69
70
71
    result = (5.472 + 0.671e3 / wavelength
              - 9.218e3 / wavelength ** 2
              + 2.620e6 / wavelength ** 3)
72
73
74
75

    return result


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
162
163
164
165
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


166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
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.

    """

    parameter_list = {
        "E_BVs_young": (
            "float",
            "E(B-V)*, the colour excess of the stellar continuum light for "
            "the young population.",
            None
        ),
        "E_BVs_old_factor": (
            "float",
            "Reduction factor for the E(B-V)* of the old population compared "
            "to the young one (<1).",
            None
        ),
        "young_contribution_name": (
            "string",
            "Name of the contribution containing the spectrum of the "
            "young population.",
Yannick Roehlly's avatar
Yannick Roehlly committed
192
            "m2005_young"
193
194
195
196
197
198
        ),
        "old_contribution_name": (
            "string",
            "Name of the contribution containing the spectrum of the "
            "old population. If it is set to 'None', only one population "
            "is considered.",
Yannick Roehlly's avatar
Yannick Roehlly committed
199
            "m2005_old"
200
        ),
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
        "uv_bump_wavelength": (
            "float",
            "Central wavelength of the UV bump in nm.",
            217.5
        ),
        "uv_bump_width": (
            "float",
            "Width (FWHM) of the UV bump in nm.",
            None
        ),
        "uv_bump_amplitude": (
            "float",
            "Amplitude of the UV bump in nm.",
            None
        ),
        "powerlaw_slope": (
            "float",
            "Slope delta of the power law modifying the attenuation curve.",
            None
        ),
221
222
223
224
225
226
227
228
229
230
        "filters": (
            "list of strings",
            "List of the filters for which the attenuation will be computed.",
            ['V_B90', 'FUV']
        )
    }

    out_parameter_list = {
        "NAME_E_BVs_young": "E(B-V)*, the colour excess of the stellar "
                            "continuum light for the young population.",
231
232
        "NAME_E_BVs_old": "E(B-V)*, the colour excess of the stellar "
                          "continuum light for the old population.",
233
234
235
236
237
238
239
        "NAME_attenuation_young": "Amount of luminosity attenuated from the "
                                  "young population in W.",
        "NAME_E_BVs_old_factor": "Ratio of the old population E(B-V)* to the "
                                 "young one.",
        "NAME_attenuation_old": "Amount of luminosity attenuated from the "
                                "old population in W.",
        "NAME_attenuation": "Total amount of luminosity attenuated in W.",
240
241
242
243
        "NAME_uv_bump_wavelength": "Central wavelength of UV bump in nm.",
        "NAME_uv_bump_width": "Width of the UV bump in nm.",
        "NAME_uv_bump_amplitude": "Amplitude of the UV bump in nm.",
        "NAME_powerlaw_slope": "Slope of the power law.",
244
245
246
        "NAME_FILTER": "Attenuation in the FILTER filter.",
    }

247
248
249
250
251
252
253
254
255
    def _init_code(self):
        """Get the filters from the database"""
        self.filters = {}
        base = Database()
        for filter_name in self.parameters["filters"]:
            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

        """

        # Base name for adding information to the SED.
265
        name = self.name or 'dustatt_calzleit_'
266
267

        wavelength = sed.wavelength_grid
268
269
270
271
272
273
274
275
276
        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
277

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

286
287
288
289
        # Compute attenuation curve
        sel_attenuation = a_vs_ebv(wavelength, uv_bump_wavelength,
                                   uv_bump_width, uv_bump_amplitude,
                                   powerlaw_slope)
290
291
292

        # Young population attenuation
        luminosity = sed.get_lumin_contribution(young_contrib)
293
294
295
        attenuated_luminosity = luminosity * 10 ** (ebvs_young
                                                    * sel_attenuation / -2.5)
        attenuated_luminosity[wavelength < 91.2] = 0
296
297
298
299
300
        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)

301
        sed.add_module(name, self.parameters)
302
303
304
305
306
307
308
        sed.add_info(name + "_E_BVs_young", ebvs_young)
        sed.add_info(name + "_attenuation_young", attenuation_young)
        sed.add_contribution(name + "_young", wavelength, attenuation_spectrum)

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

            sed.add_info(name + "_E_BVs_old", ebvs_old)
            sed.add_info(name + "_E_BVs_old_factor",
318
                         self.parameters["E_BVs_old_factor"])
319
320
321
322
323
324
            sed.add_info(name + "_attenuation_old", attenuation_old)
            sed.add_contribution(name + "_old",
                                 wavelength, attenuation_spectrum)
        else:
            attenuation_old = 0

325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
        # 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)
341
342
343
        sed.add_info(name + "_attenuation",
                     attenuation_young + attenuation_old)

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

        # Attenuation in each filter
353
        for filter_name in filters:
354
355
356
            sed.add_info(name + "_" + filter_name,
                         -2.5 * np.log(flux_att[filter_name] /
                                       flux_noatt[filter_name]))