casey2012.py 4.41 KB
Newer Older
1
# -*- coding: utf-8 -*-
2
# Copyright (C) 2013 Centre de données Astrophysiques de Marseille
Yannick Roehlly's avatar
Yannick Roehlly committed
3
# Copyright (C) 2013 Institute of Astronomy, University of Cambridge
4
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
Yannick Roehlly's avatar
Yannick Roehlly committed
5
# Author: Médéric Boquien
6

Yannick Roehlly's avatar
Yannick Roehlly committed
7 8 9 10 11 12 13 14
"""
Casey (2012) IR models module
=============================

This module implements the Casey (2012) infra-red models.

"""

15 16
import numpy as np
import scipy.constants as cst
17
from collections import OrderedDict
18
from . import CreationModule
19 20


Yannick Roehlly's avatar
Yannick Roehlly committed
21
class Casey2012(CreationModule):
Yannick Roehlly's avatar
Yannick Roehlly committed
22
    """Casey (2012) templates IR re-emission
23 24 25 26 27 28 29

    Given an amount of attenuation (e.g. resulting from the action of a dust
    attenuation module) this module normalises the Casey (2012) template
    corresponding to a given α to this amount of energy and add it to the SED.

    """

30 31
    parameter_list = OrderedDict([
        ("temperature", (
32 33 34
            "float",
            "Temperature of the dust in K.",
            None
35 36
        )),
        ("beta", (
37 38 39
            "float",
            "Emissivity index of the dust.",
            None
40 41
        )),
        ("alpha", (
42 43 44
            "float",
            "Mid-infrared powerlaw slope.",
            None
45
        )),
46 47 48 49 50 51 52
        ('attenuation_value_keys', (
            'string',
            "Keys of the SED information dictionary where the module will "
            "look for the attenuation (in W) to re-emit. You can give several "
            "keys separated with a & (don't use commas), a re-emission "
            "contribution will be added for each key.",
            "attenuation"
53 54 55 56 57 58 59 60
        ))
    ])

    out_parameter_list = OrderedDict([
        ("temperature", "Temperature of the dust in K."),
        ("beta", "Emissivity index of the dust."),
        ("alpha", "Mid-infrared powerlaw slope.")
    ])
61 62 63 64


    def _init_code(self):
        """Build the model for a given set of parameters."""
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 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111
        T = float(self.parameters["temperature"])
        beta = float(self.parameters["beta"])
        alpha = float(self.parameters["alpha"])

        # We define various constants necessary to compute the model. For
        # consistency, we define the speed of light in nm s¯¹ rather than in
        # m s¯¹.
        c = cst.c * 1e9
        b1 = 26.68
        b2 = 6.246
        b3 = 1.905e-4
        b4 = 7.243e-5
        lambda_c = 0.75e3 / ((b1 + b2 * alpha) ** -2. + (b3 + b4 * alpha) * T)
        lambda_0 = 200e3
        Npl = ((1. - np.exp(-(lambda_0 / lambda_c) ** beta)) * (c / lambda_c)
              ** 3. / (np.exp(cst.h * c / (lambda_c * cst.k * T)) - 1.))

        self.wave = np.logspace(3., 6., 1000.)
        conv = c / (self.wave * self.wave)

        self.lumin_blackbody = conv * (1. - np.exp(-(lambda_0 / self.wave)
                              ** beta)) * (c / self.wave) ** 3. / (np.exp(
                              cst.h * c / (self.wave * cst.k * T)) - 1.)
        self.lumin_powerlaw = (conv * Npl * (self.wave / lambda_c) ** alpha *
                        np.exp(-(self.wave / lambda_c) ** 2.))

        # TODO, save the right normalisation factor to retrieve the dust mass
        norm = np.trapz(self.lumin_powerlaw + self.lumin_blackbody,
                        x=self.wave)
        self.lumin_powerlaw /= norm
        self.lumin_blackbody /= norm
        self.lumin = self.lumin_powerlaw + self.lumin_blackbody


    def process(self, sed):
        """Add the IR re-emission contributions.

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

        """

        # Base name for adding information to the SED.
        name = self.name or 'casey2012'

112 113 114 115
        attenuation_value_keys = [
            item.strip() for item in
            self.parameters["attenuation_value_keys"].split("&")]

116
        sed.add_module(name, self.parameters)
117 118 119
        sed.add_info("temperature" + self.postfix, self.parameters["temperature"])
        sed.add_info("alpha" + self.postfix, self.parameters["alpha"])
        sed.add_info("beta" + self.postfix, self.parameters["beta"])
120

121
        for attenuation in attenuation_value_keys:
122 123 124 125 126 127 128 129 130 131
            sed.add_contribution(
                name + '_powerlaw_' + attenuation,
                self.wave,
                sed.info[attenuation] * self.lumin_powerlaw
            )
            sed.add_contribution(
                name + '_blackbody_' + attenuation,
                self.wave,
                sed.info[attenuation] * self.lumin_blackbody
            )
Yannick Roehlly's avatar
Yannick Roehlly committed
132 133 134

# CreationModule to be returned by get_module
Module = Casey2012