dustatt_powerlaw.py 10.3 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
        )),
        ("young_contribution_name", (
121 122 123
            "string",
            "Name of the contribution containing the spectrum of the "
            "young population.",
124
            "ssp_young"
125 126
        )),
        ("old_contribution_name", (
127 128 129 130
            "string",
            "Name of the contribution containing the spectrum of the "
            "old population. If it is set to 'None', only one population "
            "is considered.",
131
            "ssp_old"
132 133
        )),
        ("uv_bump_wavelength", (
134 135 136
            "float",
            "Central wavelength of the UV bump in nm.",
            217.5
137 138
        )),
        ("uv_bump_width", (
139 140 141
            "float",
            "Width (FWHM) of the UV bump in nm.",
            None
142 143
        )),
        ("uv_bump_amplitude", (
144 145 146
            "float",
            "Amplitude of the UV bump in nm.",
            None
147 148
        )),
        ("powerlaw_slope", (
149 150 151
            "float",
            "Slope delta of the power law continuum.",
            -0.7
152 153
        )),
        ("filters", (
154 155 156 157 158
            "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"
159 160 161 162
        ))
    ])

    out_parameter_list = OrderedDict([
163 164 165 166 167 168 169 170 171 172 173 174 175 176
        ("Av_young", "V-band attenuation of the young population."),
        ("Av_old", "V-band attenuation of the old population."),
        ("attenuation_young", "Amount of luminosity attenuated from the "
                              "young population in W."),
        ("Av_old_factor", "Reduction factor for the V-band attenuation "
                          "of  the old population compared to the young "
                          "one (<1)."),
        ("attenuation_old", "Amount of luminosity attenuated from the "
                            "old population in W."),
        ("attenuation", "Total amount of luminosity attenuated in W."),
        ("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."),
177
        ("FILTER_attenuation", "Attenuation in the FILTER filter.")
178
    ])
179

180 181
    def _init_code(self):
        """Get the filters from the database"""
182 183
        filter_list = [item.strip() for item in
                       self.parameters["filters"].split("&")]
184 185
        self.filters = {}
        base = Database()
186
        for filter_name in filter_list:
187 188 189
            self.filters[filter_name] = base.get_filter(filter_name)
        base.close()

190
    def process(self, sed):
191 192 193 194 195 196 197 198 199
        """Add the CCM dust attenuation to the SED.

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

        """

        wavelength = sed.wavelength_grid
200 201 202 203 204 205 206 207
        av_young = float(self.parameters["Av_young"])
        av_old = float(self.parameters["Av_old_factor"] * av_young)
        young_contrib = self.parameters["young_contribution_name"]
        old_contrib = self.parameters["old_contribution_name"]
        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"])
208
        filters = self.filters
209

210
        # Fλ fluxes (only from continuum)) in each filter before attenuation.
211
        flux_noatt = {}
212
        for filter_name, filter_ in filters.items():
213
            flux_noatt[filter_name] = sed.compute_fnu(
214
                filter_.trans_table,
215 216
                filter_.effective_wavelength,
                add_line_fluxes=False)
217 218

        # Compute attenuation curve
Yannick Roehlly's avatar
Yannick Roehlly committed
219 220 221
        sel_attenuation = alambda_av(wavelength, powerlaw_slope,
                                     uv_bump_wavelength, uv_bump_width,
                                     uv_bump_amplitude)
222 223 224 225 226 227 228 229 230 231

        # Young population attenuation
        luminosity = sed.get_lumin_contribution(young_contrib)
        attenuated_luminosity = luminosity * 10 ** (av_young
                                                    * sel_attenuation / -2.5)
        attenuation_spectrum = attenuated_luminosity - luminosity
        # We integrate the amount of luminosity attenuated (-1 because the
        # spectrum is negative).
        attenuation_young = -1 * np.trapz(attenuation_spectrum, wavelength)

232 233 234 235 236
        sed.add_module(self.name, self.parameters)
        sed.add_info("Av_young" + self.postfix, av_young)
        sed.add_info("attenuation_young" + self.postfix, attenuation_young)
        sed.add_contribution("attenuation_young" + self.postfix, wavelength,
                             attenuation_spectrum)
237 238 239 240 241 242 243 244 245 246

        # Old population (if any) attenuation
        if old_contrib:
            luminosity = sed.get_lumin_contribution(old_contrib)
            attenuated_luminosity = (luminosity *
                                     10 ** (av_old
                                            * sel_attenuation / -2.5))
            attenuation_spectrum = attenuated_luminosity - luminosity
            attenuation_old = -1 * np.trapz(attenuation_spectrum, wavelength)

247 248
            sed.add_info("Av_old" + self.postfix, av_old)
            sed.add_info("Av_old_factor" + self.postfix,
249
                         self.parameters["Av_old_factor"])
250
            sed.add_info("attenuation_old" + self.postfix, attenuation_old)
Yannick Roehlly's avatar
Yannick Roehlly committed
251
            sed.add_contribution("attenuation_old" + self.postfix,
252 253 254 255
                                 wavelength, attenuation_spectrum)
        else:
            attenuation_old = 0

256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271
        # Attenuation of spectral lines
        if young_contrib in sed.lines:
            line_attenuation = alambda_av(sed.lines[young_contrib][0],
                                          powerlaw_slope, uv_bump_wavelength,
                                          uv_bump_width, uv_bump_amplitude)
            sed.lines[young_contrib][1] *= 10 ** (av_young *
                                                  line_attenuation / -2.5)
        if old_contrib in sed.lines:
            line_attenuation = alambda_av(sed.lines[old_contrib][0],
                                          powerlaw_slope, uv_bump_wavelength,
                                          uv_bump_width, uv_bump_amplitude)
            sed.lines[old_contrib][1] *= 10 ** (av_old *
                                                line_attenuation / -2.5)

        # Total attenuation (we don't take into account the energy attenuated
        # in the spectral lines)
272
        sed.add_info("attenuation" + self.postfix,
273 274
                     attenuation_young + attenuation_old)

275
        # Fλ fluxes (only in continuum) in each filter after attenuation.
276
        flux_att = {}
277 278 279
        for filter_name, filter_ in filters.items():
            flux_att[filter_name] = sed.compute_fnu(
                filter_.trans_table,
280 281
                filter_.effective_wavelength,
                add_line_fluxes=False)
282 283

        # Attenuation in each filter
284
        for filter_name in filters:
285
            sed.add_info(filter_name + "_attenuation" + self.postfix,
286 287
                         -2.5 * np.log(flux_att[filter_name] /
                                       flux_noatt[filter_name]))
Yannick Roehlly's avatar
Yannick Roehlly committed
288 289 290

# CreationModule to be returned by get_module
Module = PowerLawAtt