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


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
67
68
69
70
    """Compute the complete attenuation curve A(λ)/Av

    The continuum is a power law (λ / λv) ** δ to with is added a UV bump.

    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
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
        """Add the CCM dust attenuation to the SED.

        Parameters
        ----------
177
        sed: pcigale.sed.SED object
178
179

        """
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
                filter_.effective_wavelength)
196
197

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

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

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

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

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

# CreationModule to be returned by get_module
Module = PowerLawAtt