dustatt_calzleit.py 10.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

Yannick Roehlly's avatar
Yannick Roehlly committed
7 8 9 10 11 12 13 14 15
"""
Calzetti et al. (2000) and Leitherer et al. (2002) attenuation module
=====================================================================

This module implements the Calzetti et al. (2000) and  Leitherer et al. (2002)
attenuation formulae, adding an UV-bump and a power law.

"""

16
from collections import OrderedDict
17

18
import numpy as np
19

20
from . import SedModule
21 22 23 24 25 26


def k_calzetti2000(wavelength):
    """Compute the Calzetti et al. (2000) A(λ)/E(B-V)∗

    Given a wavelength grid, this function computes the selective attenuation
27 28 29
    A(λ)/E(B-V)∗ using the formula from Calzetti at al. (2000). This formula
    is given for wavelengths between 120 nm and 2200 nm, but this function
    makes the computation outside.
30 31 32

    Parameters
    ----------
33
    wavelength: array of floats
34 35 36 37 38 39 40 41 42 43 44
        Wavelength grid in nm.

    Returns
    -------
    a numpy array of floats

    """
    wavelength = np.array(wavelength)
    result = np.zeros(len(wavelength))

    # Attenuation between 120 nm and 630 nm
45
    mask = (wavelength < 630)
Médéric Boquien's avatar
Médéric Boquien committed
46 47 48
    result[mask] = 2.659 * (-2.156 + 1.509e3 / wavelength[mask] -
                            0.198e6 / wavelength[mask] ** 2 +
                            0.011e9 / wavelength[mask] ** 3) + 4.05
49 50

    # Attenuation between 630 nm and 2200 nm
51
    mask = (wavelength >= 630)
52 53 54 55 56 57 58 59 60
    result[mask] = 2.659 * (-1.857 + 1.040e3 / wavelength[mask]) + 4.05

    return result


def k_leitherer2002(wavelength):
    """Compute the Leitherer et al. (2002) A(λ)/E(B-V)∗

    Given a wavelength grid, this function computes the selective attenuation
61 62 63
    A(λ)/E(B-V)∗ using the formula from Leitherer at al. (2002). This formula
    is given for wavelengths between 91.2 nm and 180 nm, but this function
    makes the computation outside.
64 65 66

    Parameters
    ----------
67
    wavelength: array of floats
68 69 70 71 72 73 74 75
        Wavelength grid in nm.

    Returns
    -------
    a numpy array of floats

    """
    wavelength = np.array(wavelength)
Médéric Boquien's avatar
Médéric Boquien committed
76 77 78
    result = (5.472 + 0.671e3 / wavelength -
              9.218e3 / wavelength ** 2 +
              2.620e6 / wavelength ** 3)
79 80 81 82

    return result


83 84 85 86 87
def uv_bump(wavelength, central_wave, gamma, ebump):
    """Compute the Lorentzian-like Drude profile.

    Parameters
    ----------
88
    wavelength: array of floats
89
        Wavelength grid in nm.
90
    central_wave: float
91
        Central wavelength of the bump in nm.
92
    gamma: float
93
        Width (FWHM) of the bump in nm.
94
    ebump: float
95 96 97 98 99 100 101 102
        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
103 104
            ((wavelength ** 2 - central_wave ** 2) ** 2 +
             wavelength ** 2 * gamma ** 2))
105 106 107 108 109 110 111


def power_law(wavelength, delta):
    """Power law 'centered' on 550 nm..

    Parameters
    ----------
112
    wavelength: array of floats
113
        The wavelength grid in nm.
114
    delta: float
115 116 117 118 119 120 121 122 123 124 125 126 127
        The slope of the power law.

    Returns
    -------
    array of floats

    """
    return (wavelength / 550) ** delta


def a_vs_ebv(wavelength, bump_wave, bump_width, bump_ampl, power_slope):
    """Compute the complete attenuation curve A(λ)/E(B-V)*

128 129 130 131 132
    The Leitherer et al. (2002) formula is used between 91.2 nm and 150 nm, and
    the Calzetti et al. (2000) formula is used after 150 (we do an
    extrapolation after 2200 nm). When the attenuation becomes negative, it is
    kept to 0. This continuum is multiplied by the power law and then the UV
    bump is added.
133 134 135

    Parameters
    ----------
136
    wavelength: array of floats
137
        The wavelength grid (in nm) to compute the attenuation curve on.
138
    bump_wave: float
139
        Central wavelength (in nm) of the UV bump.
140
    bump_width: float
141
        Width (FWHM, in nm) of the UV bump.
142
    bump_ampl: float
143
        Amplitude of the UV bump.
144
    power_slope: float
145 146 147 148
        Slope of the power law.

    Returns
    -------
149
    attenuation: array of floats
150 151 152 153 154 155
        The A(λ)/E(B-V)* attenuation at each wavelength of the grid.

    """
    attenuation = np.zeros(len(wavelength))

    # Leitherer et al.
156
    mask = (wavelength > 91.2) & (wavelength < 150)
157 158 159 160 161 162 163 164
    attenuation[mask] = k_leitherer2002(wavelength[mask])
    # Calzetti et al.
    mask = (wavelength >= 150)
    attenuation[mask] = k_calzetti2000(wavelength[mask])
    # We set attenuation to 0 where it becomes negative
    mask = (attenuation < 0)
    attenuation[mask] = 0
    # Power law
165
    attenuation *= power_law(wavelength, power_slope)
166

167 168 169 170 171 172 173 174 175 176
    # As the powerlaw slope changes E(B-V), we correct this so that the curve
    # always has the same E(B-V) as the starburst curve. This ensures that the
    # E(B-V) requested by the user is the actual E(B-V) of the curve.
    wl_BV = np.array([440., 550.])
    EBV_calz = ((k_calzetti2000(wl_BV) * power_law(wl_BV, 0.)) +
                uv_bump(wl_BV, bump_wave, bump_width, bump_ampl))
    EBV = ((k_calzetti2000(wl_BV) * power_law(wl_BV, power_slope)) +
           uv_bump(wl_BV, bump_wave, bump_width, bump_ampl))
    attenuation *= (EBV_calz[1]-EBV_calz[0]) / (EBV[1]-EBV[0])

177 178 179 180 181
    # UV bump. It is added after the renormalization as the bump strength
    # should correspond to the requested E(B-V) and should therefore not be
    # changed by the renormalization.
    attenuation += uv_bump(wavelength, bump_wave, bump_width, bump_ampl)

182 183 184
    return attenuation


185
class CalzLeit(SedModule):
Yannick Roehlly's avatar
Yannick Roehlly committed
186 187
    """Calzetti + Leitherer attenuation module

188 189 190 191
    This module computes the dust attenuation using the formulae from
    Calzetti et al. (2000) and Leitherer et al. (2002). Note that both the
    stars and the gas are attenuated with the same curve as opposed to Calzetti
    et al. (2000) where the gas is attenuated with a Milky Way curve.
192

Yannick Roehlly's avatar
Yannick Roehlly committed
193 194
    The attenuation can be computed on the whole spectrum or on a specific
    contribution and is added to the SED as a negative contribution.
195 196 197

    """

198
    parameter_list = OrderedDict([
199
        ("E_BVs_young", (
200
            "cigale_list(minvalue=0.)",
201
            "E(B-V)*, the colour excess of the stellar continuum light for "
202
            "the young population.",
203
            0.3
204 205
        )),
        ("E_BVs_old_factor", (
206
            "cigale_list(minvalue=0., maxvalue=1.)",
207 208
            "Reduction factor for the E(B-V)* of the old population compared "
            "to the young one (<1).",
209
            1.0
210 211
        )),
        ("uv_bump_wavelength", (
212
            "cigale_list(minvalue=0.)",
213 214
            "Central wavelength of the UV bump in nm.",
            217.5
215 216
        )),
        ("uv_bump_width", (
217
            "cigale_list()",
218
            "Width (FWHM) of the UV bump in nm.",
219
            35.
220 221
        )),
        ("uv_bump_amplitude", (
222
            "cigale_list(minvalue=0.)",
223 224
            "Amplitude of the UV bump. For the Milky Way: 3.",
            0.
225 226
        )),
        ("powerlaw_slope", (
227
            "cigale_list()",
228
            "Slope delta of the power law modifying the attenuation curve.",
229
            0.
230 231
        )),
        ("filters", (
232
            "string()",
233 234 235
            "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).",
236
            "B_B90 & V_B90 & FUV"
237 238 239
        ))
    ])

240 241
    def _init_code(self):
        """Get the filters from the database"""
242 243 244 245 246 247 248 249 250
        self.ebvs = {}
        self.ebvs['young'] = float(self.parameters["E_BVs_young"])
        self.ebvs_old_factor = float(self.parameters["E_BVs_old_factor"])
        self.ebvs['old'] = self.ebvs_old_factor * self.ebvs['young']
        self.uv_bump_wavelength = float(self.parameters["uv_bump_wavelength"])
        self.uv_bump_width = float(self.parameters["uv_bump_width"])
        self.uv_bump_amplitude = float(self.parameters["uv_bump_amplitude"])
        self.powerlaw_slope = float(self.parameters["powerlaw_slope"])

251 252
        self.filter_list = [item.strip() for item in
                            self.parameters["filters"].split("&")]
253 254
        # We cannot compute the attenuation until we know the wavelengths. Yet,
        # we reserve the object.
255
        self.curve = {}
256 257

    def process(self, sed):
258 259 260 261
        """Add the CCM dust attenuation to the SED.

        Parameters
        ----------
262
        sed: pcigale.sed.SED object
263 264 265 266

        """
        wavelength = sed.wavelength_grid

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

270
        # Compute attenuation curve
271 272 273 274 275 276
        if len(self.curve) == 0:
            sel_att = a_vs_ebv(wavelength, self.uv_bump_wavelength,
                               self.uv_bump_width, self.uv_bump_amplitude,
                               self.powerlaw_slope)
            for age in ['young', 'old']:
                self.curve[age] = 10 ** (-.4 * self.ebvs[age] * sel_att)
277

278
        attenuation_total = 0.
279 280 281
        contribs = [contrib for contrib in sed.contribution_names if
                    'absorption' not in contrib]
        for contrib in contribs:
282 283
            age = contrib.split('.')[-1].split('_')[-1]
            luminosity = sed.get_lumin_contribution(contrib)
284
            attenuation_spectrum = luminosity * (self.curve[age] - 1.)
285 286
            # We integrate the amount of luminosity attenuated (-1 because the
            # spectrum is negative).
287
            attenuation = -1. * np.trapz(attenuation_spectrum, wavelength)
288 289 290
            attenuation_total += attenuation

            sed.add_module(self.name, self.parameters)
291
            sed.add_info("attenuation.E_BVs." + contrib, self.ebvs[age])
292
            sed.add_info("attenuation." + contrib, attenuation, True)
293 294 295 296
            sed.add_contribution("attenuation." + contrib, wavelength,
                                 attenuation_spectrum)

        # Total attenuation
297 298 299 300 301 302
        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)
303

304
        # Fλ fluxes (only from continuum) in each filter after attenuation.
Médéric Boquien's avatar
Médéric Boquien committed
305
        flux_att = {filt: sed.compute_fnu(filt) for filt in self.filter_list}
306 307

        # Attenuation in each filter
308 309 310
        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
311

312 313 314 315 316
        sed.add_info('attenuation.ebvs_old_factor', self.ebvs_old_factor)
        sed.add_info('attenuation.uv_bump_wavelength', self.uv_bump_wavelength)
        sed.add_info('attenuation.uv_bump_width', self.uv_bump_width)
        sed.add_info('attenuation.uv_bump_amplitude', self.uv_bump_amplitude)
        sed.add_info('attenuation.powerlaw_slope', self.powerlaw_slope)
317

318

319
# SedModule to be returned by get_module
Yannick Roehlly's avatar
Yannick Roehlly committed
320
Module = CalzLeit