dustatt_powerlaw.py 9.95 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
"""
7
8
9
10
11
12
This module implements the attenuation base on a power law as defined
in Charlot and Fall (2000) with a UV bump added.

"""

import numpy as np
13
from collections import OrderedDict
14
from . import common
15
from ..data import Database
16
17
18
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


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
61
def alambda_av(wavelength, delta, bump_wave, bump_width, bump_ampl):
62
63
64
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
97
    """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


class Module(common.SEDCreationModule):
    """Add CCM dust attenuation based on Charlot and Fall (2000) power law.
    """

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

    out_parameter_list = OrderedDict([
153
154
155
156
157
158
159
160
161
162
163
164
165
166
        ("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."),
167
        ("FILTER_attenuation", "Attenuation in the FILTER filter.")
168
    ])
169

170
171
    def _init_code(self):
        """Get the filters from the database"""
172
173
        filter_list = [item.strip() for item in
                       self.parameters["filters"].split("&")]
174
175
        self.filters = {}
        base = Database()
176
        for filter_name in filter_list:
177
178
179
            self.filters[filter_name] = base.get_filter(filter_name)
        base.close()

180
    def process(self, sed):
181
182
183
184
185
186
187
188
189
        """Add the CCM dust attenuation to the SED.

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

        """

        wavelength = sed.wavelength_grid
190
191
192
193
194
195
196
197
        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"])
198
        filters = self.filters
199

200
        # Fλ fluxes (only from continuum)) in each filter before attenuation.
201
        flux_noatt = {}
202
        for filter_name, filter_ in filters.items():
203
            flux_noatt[filter_name] = sed.compute_fnu(
204
                filter_.trans_table,
205
206
                filter_.effective_wavelength,
                add_line_fluxes=False)
207
208

        # Compute attenuation curve
Yannick Roehlly's avatar
Yannick Roehlly committed
209
210
211
        sel_attenuation = alambda_av(wavelength, powerlaw_slope,
                                     uv_bump_wavelength, uv_bump_width,
                                     uv_bump_amplitude)
212
213
214
215
216
217
218
219
220
221

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

222
223
224
225
226
        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)
227
228
229
230
231
232
233
234
235
236

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

237
238
            sed.add_info("Av_old" + self.postfix, av_old)
            sed.add_info("Av_old_factor" + self.postfix,
239
                         self.parameters["Av_old_factor"])
240
            sed.add_info("attenuation_old" + self.postfix, attenuation_old)
Yannick Roehlly's avatar
Yannick Roehlly committed
241
            sed.add_contribution("attenuation_old" + self.postfix,
242
243
244
245
                                 wavelength, attenuation_spectrum)
        else:
            attenuation_old = 0

246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
        # 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)
262
        sed.add_info("attenuation" + self.postfix,
263
264
                     attenuation_young + attenuation_old)

265
        # Fλ fluxes (only in continuum) in each filter after attenuation.
266
        flux_att = {}
267
268
269
        for filter_name, filter_ in filters.items():
            flux_att[filter_name] = sed.compute_fnu(
                filter_.trans_table,
270
271
                filter_.effective_wavelength,
                add_line_fluxes=False)
272
273

        # Attenuation in each filter
274
        for filter_name in filters:
275
            sed.add_info(filter_name + "_attenuation" + self.postfix,
276
277
                         -2.5 * np.log(flux_att[filter_name] /
                                       flux_noatt[filter_name]))