dustatt_powerlaw.py

# -*- coding: utf-8 -*-
# Copyright (C) 2013 Centre de données Astrophysiques de Marseille
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
# Author: Yannick Roehlly

"""
Charlot and Fall (2000) power law attenuation module
====================================================

This module implements the attenuation based on a power law as defined
in Charlot and Fall (2000) with a UV bump added.

"""

import numpy as np
from collections import OrderedDict
from . import CreationModule
from ..data import Database


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))


def alambda_av(wavelength, delta, bump_wave, bump_width, bump_ampl):
    """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


class PowerLawAtt(CreationModule):
    """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.

    """

    parameter_list = OrderedDict([
        ("Av_young", (
            "float",
            "V-band attenuation of the young population.",
            None
        )),
        ("Av_old_factor", (
            "float",
            "Reduction factor for the V-band attenuation of the old "
            "population compared to the young one (<1).",
            None
        )),
        ("young_contribution_name", (
            "string",
            "Name of the contribution containing the spectrum of the "
            "young population.",
            "ssp_young"
        )),
        ("old_contribution_name", (
            "string",
            "Name of the contribution containing the spectrum of the "
            "old population. If it is set to 'None', only one population "
            "is considered.",
            "ssp_old"
        )),
        ("uv_bump_wavelength", (
            "float",
            "Central wavelength of the UV bump in nm.",
            217.5
        )),
        ("uv_bump_width", (
            "float",
            "Width (FWHM) of the UV bump in nm.",
            None
        )),
        ("uv_bump_amplitude", (
            "float",
            "Amplitude of the UV bump in nm.",
            None
        )),
        ("powerlaw_slope", (
            "float",
            "Slope delta of the power law continuum.",
            -0.7
        )),
        ("filters", (
            "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"
        ))
    ])

    out_parameter_list = OrderedDict([
        ("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."),
        ("FILTER_attenuation", "Attenuation in the FILTER filter.")
    ])

    def _init_code(self):
        """Get the filters from the database"""
        filter_list = [item.strip() for item in
                       self.parameters["filters"].split("&")]
        self.filters = {}
        base = Database()
        for filter_name in filter_list:
            self.filters[filter_name] = base.get_filter(filter_name)
        base.close()

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

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

        """

        wavelength = sed.wavelength_grid
        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"])
        filters = self.filters

        # Fλ fluxes (only from continuum)) in each filter before attenuation.
        flux_noatt = {}
        for filter_name, filter_ in filters.items():
            flux_noatt[filter_name] = sed.compute_fnu(
                filter_.trans_table,
                filter_.effective_wavelength,
                add_line_fluxes=False)

        # Compute attenuation curve
        sel_attenuation = alambda_av(wavelength, powerlaw_slope,
                                     uv_bump_wavelength, uv_bump_width,
                                     uv_bump_amplitude)

        # 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)

        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)

        # 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)

            sed.add_info("Av_old" + self.postfix, av_old)
            sed.add_info("Av_old_factor" + self.postfix,
                         self.parameters["Av_old_factor"])
            sed.add_info("attenuation_old" + self.postfix, attenuation_old)
            sed.add_contribution("attenuation_old" + self.postfix,
                                 wavelength, attenuation_spectrum)
        else:
            attenuation_old = 0

        # 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)
        sed.add_info("attenuation" + self.postfix,
                     attenuation_young + attenuation_old)

        # Fλ fluxes (only in continuum) in each filter after attenuation.
        flux_att = {}
        for filter_name, filter_ in filters.items():
            flux_att[filter_name] = sed.compute_fnu(
                filter_.trans_table,
                filter_.effective_wavelength,
                add_line_fluxes=False)

        # Attenuation in each filter
        for filter_name in filters:
            sed.add_info(filter_name + "_attenuation" + self.postfix,
                         -2.5 * np.log(flux_att[filter_name] /
                                       flux_noatt[filter_name]))

# CreationModule to be returned by get_module
Module = PowerLawAtt