dustatt_powerlaw.py 7.75 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 . import CreationModule
18
19
20
21
22
23
24


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

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


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

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

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

    Returns
    -------
83
    attenuation: array of floats
84
85
86
87
88
89
90
91
92
93
94
95
        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
96
class PowerLawAtt(CreationModule):
Yannick Roehlly's avatar
Yannick Roehlly committed
97
98
    """Power law attenuation module

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

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

105
106
    """

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

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

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

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

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

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

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

        # Compute attenuation curve
190
191
        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
192
193
194
                                              uv_bump_wavelength,
                                              uv_bump_width, uv_bump_amplitude)

195
196
197
198
199
        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
200
                                     (av[age] * self.sel_attenuation / -2.5))
201
            attenuation_spectrum = attenuated_luminosity - luminosity
202
203
204
205
206
207
            # 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)
208
209
            sed.add_info("attenuation.Av." + contrib, av[age])
            sed.add_info("attenuation." + contrib, attenuation, True)
210
211
212
            sed.add_contribution("attenuation." + contrib, wavelength,
                                 attenuation_spectrum)

213
        # Bump and slope of the dust attenuation
214
215
        sed.add_info("attenuation.uv_bump_amplitude", uv_bump_amplitude)
        sed.add_info("attenuation.powerlaw_slope", powerlaw_slope)
216

217
        # Total attenuation
218
219
220
221
222
223
        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)
224

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

        # Attenuation in each filter
229
230
231
        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
232
233
234

# CreationModule to be returned by get_module
Module = PowerLawAtt