dustatt_powerlaw.py 10.3 KB
Newer Older
1
# -*- coding: utf-8 -*-
2
3
# Copyright (C) 2013 Centre de données Astrophysiques de Marseille
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
Yannick Roehlly's avatar
Yannick Roehlly committed
4
# Author: Yannick Roehlly
5

6
"""
Yannick Roehlly's avatar
Yannick Roehlly committed
7
8
9
10
Charlot and Fall (2000) power law attenuation module
====================================================

This module implements the attenuation based on a power law as defined
11
12
13
14
15
in Charlot and Fall (2000) with a UV bump added.

"""

import numpy as np
16
from collections import OrderedDict
17
from . import CreationModule
18
from ..data import Database
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63


def power_law(wavelength, delta):
    """Compute the power law (λ / λv)^δ

    Parameters
    ----------
    wavelength : array of float
        Wavelength grid in nm.
    delta : float
        Power law slope.

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

    """
    wave = np.array(wavelength)
    return (wave / 550) ** delta


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))


Yannick Roehlly's avatar
Yannick Roehlly committed
64
def alambda_av(wavelength, delta, bump_wave, bump_width, bump_ampl):
65
66
67
68
69
70
71
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
    """Compute the complete attenuation curve A(λ)/Av

    The continuum is a power law (λ / λv) ** δ to with is added a UV bump.

    Parameters
    ----------
    wavelength : array of floats
        The wavelength grid (in nm) to compute the attenuation curve on.
    delta : float
        Slope of the power law.
    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.

    Returns
    -------
    attenuation : array of floats
        The A(λ)/Av attenuation at each wavelength of the grid.

    """
    wave = np.array(wavelength)

    attenuation = power_law(wave, delta)
    attenuation = attenuation + uv_bump(wavelength, bump_wave,
                                        bump_width, bump_ampl)

    return attenuation


Yannick Roehlly's avatar
Yannick Roehlly committed
97
class PowerLawAtt(CreationModule):
Yannick Roehlly's avatar
Yannick Roehlly committed
98
99
100
101
102
103
104
105
    """Power law attenuation module

    This module computes the Cardelli, Clayton and Mathis attenuation using
    a power law as defined in Charlot and Fall (2000).

    The attenuation can be computed on the whole spectrum or on a specific
    contribution and is added to the SED as a negative contribution.

106
107
    """

108
109
    parameter_list = OrderedDict([
        ("Av_young", (
110
111
112
            "float",
            "V-band attenuation of the young population.",
            None
113
114
        )),
        ("Av_old_factor", (
115
116
117
118
            "float",
            "Reduction factor for the V-band attenuation of the old "
            "population compared to the young one (<1).",
            None
119
120
        )),
        ("young_contribution_name", (
121
122
123
            "string",
            "Name of the contribution containing the spectrum of the "
            "young population.",
124
            "ssp_young"
125
126
        )),
        ("old_contribution_name", (
127
128
129
130
            "string",
            "Name of the contribution containing the spectrum of the "
            "old population. If it is set to 'None', only one population "
            "is considered.",
131
            "ssp_old"
132
133
        )),
        ("uv_bump_wavelength", (
134
135
136
            "float",
            "Central wavelength of the UV bump in nm.",
            217.5
137
138
        )),
        ("uv_bump_width", (
139
140
141
            "float",
            "Width (FWHM) of the UV bump in nm.",
            None
142
143
        )),
        ("uv_bump_amplitude", (
144
145
146
            "float",
            "Amplitude of the UV bump in nm.",
            None
147
148
        )),
        ("powerlaw_slope", (
149
150
151
            "float",
            "Slope delta of the power law continuum.",
            -0.7
152
153
        )),
        ("filters", (
154
155
156
157
158
            "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"
159
160
161
162
        ))
    ])

    out_parameter_list = OrderedDict([
163
164
165
166
167
168
169
170
171
172
173
174
175
176
        ("Av_young", "V-band attenuation of the young population."),
        ("Av_old", "V-band attenuation of the old population."),
        ("attenuation_young", "Amount of luminosity attenuated from the "
                              "young population in W."),
        ("Av_old_factor", "Reduction factor for the V-band attenuation "
                          "of  the old population compared to the young "
                          "one (<1)."),
        ("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."),
177
        ("FILTER_attenuation", "Attenuation in the FILTER filter.")
178
    ])
179

180
181
    def _init_code(self):
        """Get the filters from the database"""
182
183
        filter_list = [item.strip() for item in
                       self.parameters["filters"].split("&")]
184
185
        self.filters = {}
        base = Database()
186
        for filter_name in filter_list:
187
188
189
            self.filters[filter_name] = base.get_filter(filter_name)
        base.close()

190
    def process(self, sed):
191
192
193
194
195
196
197
198
199
        """Add the CCM dust attenuation to the SED.

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

        """

        wavelength = sed.wavelength_grid
200
201
202
203
204
205
206
207
        av_young = float(self.parameters["Av_young"])
        av_old = float(self.parameters["Av_old_factor"] * av_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"])
208
        filters = self.filters
209

210
        # Fλ fluxes (only from continuum)) in each filter before attenuation.
211
        flux_noatt = {}
212
        for filter_name, filter_ in filters.items():
213
            flux_noatt[filter_name] = sed.compute_fnu(
214
                filter_.trans_table,
215
216
                filter_.effective_wavelength,
                add_line_fluxes=False)
217
218

        # Compute attenuation curve
Yannick Roehlly's avatar
Yannick Roehlly committed
219
220
221
        sel_attenuation = alambda_av(wavelength, powerlaw_slope,
                                     uv_bump_wavelength, uv_bump_width,
                                     uv_bump_amplitude)
222
223
224
225
226
227
228
229
230
231

        # Young population attenuation
        luminosity = sed.get_lumin_contribution(young_contrib)
        attenuated_luminosity = luminosity * 10 ** (av_young
                                                    * sel_attenuation / -2.5)
        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)

232
233
234
235
236
        sed.add_module(self.name, self.parameters)
        sed.add_info("Av_young" + self.postfix, av_young)
        sed.add_info("attenuation_young" + self.postfix, attenuation_young)
        sed.add_contribution("attenuation_young" + self.postfix, wavelength,
                             attenuation_spectrum)
237
238
239
240
241
242
243
244
245
246

        # Old population (if any) attenuation
        if old_contrib:
            luminosity = sed.get_lumin_contribution(old_contrib)
            attenuated_luminosity = (luminosity *
                                     10 ** (av_old
                                            * sel_attenuation / -2.5))
            attenuation_spectrum = attenuated_luminosity - luminosity
            attenuation_old = -1 * np.trapz(attenuation_spectrum, wavelength)

247
248
            sed.add_info("Av_old" + self.postfix, av_old)
            sed.add_info("Av_old_factor" + self.postfix,
249
                         self.parameters["Av_old_factor"])
250
            sed.add_info("attenuation_old" + self.postfix, attenuation_old)
Yannick Roehlly's avatar
Yannick Roehlly committed
251
            sed.add_contribution("attenuation_old" + self.postfix,
252
253
254
255
                                 wavelength, attenuation_spectrum)
        else:
            attenuation_old = 0

256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
        # Attenuation of spectral lines
        if young_contrib in sed.lines:
            line_attenuation = alambda_av(sed.lines[young_contrib][0],
                                          powerlaw_slope, uv_bump_wavelength,
                                          uv_bump_width, uv_bump_amplitude)
            sed.lines[young_contrib][1] *= 10 ** (av_young *
                                                  line_attenuation / -2.5)
        if old_contrib in sed.lines:
            line_attenuation = alambda_av(sed.lines[old_contrib][0],
                                          powerlaw_slope, uv_bump_wavelength,
                                          uv_bump_width, uv_bump_amplitude)
            sed.lines[old_contrib][1] *= 10 ** (av_old *
                                                line_attenuation / -2.5)

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

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

        # Attenuation in each filter
284
        for filter_name in filters:
285
            sed.add_info(filter_name + "_attenuation" + self.postfix,
286
287
                         -2.5 * np.log(flux_att[filter_name] /
                                       flux_noatt[filter_name]))
Yannick Roehlly's avatar
Yannick Roehlly committed
288
289
290

# CreationModule to be returned by get_module
Module = PowerLawAtt