dustatt_powerlaw.py 7.81 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
        )),
        ("uv_bump_wavelength", (
121
122
123
            "float",
            "Central wavelength of the UV bump in nm.",
            217.5
124
125
        )),
        ("uv_bump_width", (
126
127
128
            "float",
            "Width (FWHM) of the UV bump in nm.",
            None
129
130
        )),
        ("uv_bump_amplitude", (
131
132
133
            "float",
            "Amplitude of the UV bump in nm.",
            None
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
164
    def _init_code(self):
        """Get the filters from the database"""
165
166
        filter_list = [item.strip() for item in
                       self.parameters["filters"].split("&")]
167
        self.filters = {}
168
169
170
        with Database() as base:
            for filter_name in filter_list:
                self.filters[filter_name] = base.get_filter(filter_name)
171

172
    def process(self, sed):
173
174
175
176
177
178
179
        """Add the CCM dust attenuation to the SED.

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

        """
180
        av = {}
181
        wavelength = sed.wavelength_grid
182
183
        av['young'] = float(self.parameters["Av_young"])
        av['old'] = float(self.parameters["Av_old_factor"] * av['young'])
184
185
186
187
        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"])
188
        filters = self.filters
189

190
        # Fλ fluxes (only from continuum)) in each filter before attenuation.
191
        flux_noatt = {}
192
        for filter_name, filter_ in filters.items():
193
            flux_noatt[filter_name] = sed.compute_fnu(
194
                filter_.trans_table,
195
196
                filter_.effective_wavelength,
                add_line_fluxes=False)
197
198

        # Compute attenuation curve
Yannick Roehlly's avatar
Yannick Roehlly committed
199
200
201
        sel_attenuation = alambda_av(wavelength, powerlaw_slope,
                                     uv_bump_wavelength, uv_bump_width,
                                     uv_bump_amplitude)
202

203
204
205
206
207
208
        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 **
                                    (av[age] * sel_attenuation / -2.5))
209
            attenuation_spectrum = attenuated_luminosity - luminosity
210
211
212
213
214
215
216
            # 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)
            sed.add_info("attenuation.Av." + contrib, av[age])
217
            sed.add_info("attenuation." + contrib, attenuation, True)
218
219
220
221
            sed.add_contribution("attenuation." + contrib, wavelength,
                                 attenuation_spectrum)

        # Total attenuation
222
        sed.add_info("dust.luminosity", attenuation_total, True)
223

224
        # Fλ fluxes (only in continuum) in each filter after attenuation.
225
        flux_att = {}
226
227
228
        for filter_name, filter_ in filters.items():
            flux_att[filter_name] = sed.compute_fnu(
                filter_.trans_table,
229
230
                filter_.effective_wavelength,
                add_line_fluxes=False)
231
232

        # Attenuation in each filter
233
        for filter_name in filters:
234
            sed.add_info("attenuation." + filter_name,
Yannick Roehlly's avatar
Yannick Roehlly committed
235
236
                         -2.5 * np.log10(flux_att[filter_name] /
                                         flux_noatt[filter_name]))
Yannick Roehlly's avatar
Yannick Roehlly committed
237
238
239

# CreationModule to be returned by get_module
Module = PowerLawAtt