dustatt_powerlaw.py 7.88 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
60
61
62
63
        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
    """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

100
101
    This module computes the attenuation using a power law 
    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.
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,
Yannick Roehlly's avatar
Yannick Roehlly committed
193
194
                                     uv_bump_wavelength, uv_bump_width,
                                     uv_bump_amplitude)
195
  
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 **
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)
209
210
211
212
            sed.add_info("attenuation.Av.stellar"+
                         contrib[contrib.find('_'):], av[age])
            sed.add_info("attenuation.stellar"+
                         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.
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