dustatt_powerlaw.py 7.9 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

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

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

"""

import numpy as np
17
from collections import OrderedDict
18
from . import CreationModule
19
20
21
22
23
24
25


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

    Parameters
    ----------
26
    wavelength: array of float
27
        Wavelength grid in nm.
28
    delta: float
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
        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
    ----------
45
    wavelength: array of floats
46
        Wavelength grid in nm.
47
    central_wave: float
48
        Central wavelength of the bump in nm.
49
    gamma: float
50
        Width (FWHM) of the bump in nm.
51
    ebump: float
52
53
54
55
56
57
58
59
        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
60
61
            ((wavelength ** 2 - central_wave ** 2) ** 2 +
             wavelength ** 2 * gamma ** 2))
62
63


Yannick Roehlly's avatar
Yannick Roehlly committed
64
def alambda_av(wavelength, delta, bump_wave, bump_width, bump_ampl):
65
66
    """Compute the complete attenuation curve A(λ)/Av

67
    The continuum is a power law (λ / λv) ** δ to which is added a UV bump.
68
69
70

    Parameters
    ----------
71
    wavelength: array of floats
72
        The wavelength grid (in nm) to compute the attenuation curve on.
73
    delta: float
74
        Slope of the power law.
75
    bump_wave: float
76
        Central wavelength (in nm) of the UV bump.
77
    bump_width: float
78
        Width (FWHM, in nm) of the UV bump.
79
    bump_ampl: float
80
81
82
83
        Amplitude of the UV bump.

    Returns
    -------
84
    attenuation: array of floats
85
86
87
88
89
90
91
92
93
94
95
96
        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
    """Power law attenuation module

Médéric Boquien's avatar
Médéric Boquien committed
100
    This module computes the attenuation using a power law
101
    as defined in Charlot and Fall (2000).
Yannick Roehlly's avatar
Yannick Roehlly committed
102
103
104
105

    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
            "float",
            "V-band attenuation of the young population.",
112
            1.
113
114
        )),
        ("Av_old_factor", (
115
116
117
            "float",
            "Reduction factor for the V-band attenuation of the old "
            "population compared to the young one (<1).",
118
            0.44
119
120
        )),
        ("uv_bump_wavelength", (
121
122
123
            "float",
            "Central wavelength of the UV bump in nm.",
            217.5
124
125
        )),
        ("uv_bump_width", (
126
127
            "float",
            "Width (FWHM) of the UV bump in nm.",
128
            35.
129
130
        )),
        ("uv_bump_amplitude", (
131
            "float",
132
133
            "Amplitude of the UV bump. For the Milky Way: 3.",
            0.
134
135
        )),
        ("powerlaw_slope", (
136
137
138
            "float",
            "Slope delta of the power law continuum.",
            -0.7
139
140
        )),
        ("filters", (
141
142
143
144
145
            "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"
146
147
148
149
        ))
    ])

    out_parameter_list = OrderedDict([
150
        ("Av", "V-band attenuation."),
151
152
153
        ("Av_old_factor", "Reduction factor for the V-band attenuation "
                          "of  the old population compared to the young "
                          "one (<1)."),
154
155
        ("attenuation", "Amount of luminosity attenuated in W for each "
                        "component and in total."),
156
157
158
159
        ("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."),
160
        ("FILTER_attenuation", "Attenuation in the FILTER filter.")
161
    ])
162

163
    def _init_code(self):
164
165
        self.filter_list = [item.strip() for item in
                            self.parameters["filters"].split("&")]
166
167
168
        # We cannot compute the attenuation until we know the wavelengths. Yet,
        # we reserve the object.
        self.sel_attenuation = None
169

170
    def process(self, sed):
171
        """Add the dust attenuation to the SED.
172
173
174

        Parameters
        ----------
175
        sed: pcigale.sed.SED object
176
177

        """
178
        av = {}
179
        wavelength = sed.wavelength_grid
180
181
        av['young'] = float(self.parameters["Av_young"])
        av['old'] = float(self.parameters["Av_old_factor"] * av['young'])
182
        uv_bump_wavelength = float(self.parameters["uv_bump_wavelength"])
183
        uv_bump_width = float(self.parameters["uv_bump_width"])
184
185
        uv_bump_amplitude = float(self.parameters["uv_bump_amplitude"])
        powerlaw_slope = float(self.parameters["powerlaw_slope"])
186

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

        # Compute attenuation curve
191
192
        if self.sel_attenuation is None:
            self.sel_attenuation = alambda_av(wavelength, powerlaw_slope,
Médéric Boquien's avatar
Médéric Boquien committed
193
194
195
                                              uv_bump_wavelength,
                                              uv_bump_width, uv_bump_amplitude)

196
197
198
199
200
        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
201
                                     (av[age] * self.sel_attenuation / -2.5))
202
            attenuation_spectrum = attenuated_luminosity - luminosity
203
204
205
206
207
208
            # 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)
Médéric Boquien's avatar
Médéric Boquien committed
209
            sed.add_info("attenuation.Av.stellar" +
210
                         contrib[contrib.find('_'):], av[age])
Médéric Boquien's avatar
Médéric Boquien committed
211
            sed.add_info("attenuation.stellar" +
212
                         contrib[contrib.find('_'):], attenuation, True)
213
214
215
            sed.add_contribution("attenuation." + contrib, wavelength,
                                 attenuation_spectrum)

216
        # Bump and slope of the dust attenuation
217
218
        sed.add_info("attenuation.uv_bump_amplitude", uv_bump_amplitude)
        sed.add_info("attenuation.powerlaw_slope", powerlaw_slope)
219

220
        # Total attenuation
221
222
223
224
225
226
        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)
227

228
        # Fλ fluxes (only in continuum) in each filter after attenuation.
Médéric Boquien's avatar
Médéric Boquien committed
229
        flux_att = {filt: sed.compute_fnu(filt) for filt in self.filter_list}
230
231

        # Attenuation in each filter
232
233
234
        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
235
236
237

# CreationModule to be returned by get_module
Module = PowerLawAtt