casey2012.py 3.74 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
from collections import OrderedDict
16

17 18
import numpy as np
import scipy.constants as cst
19

20
from . import SedModule
21 22


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

    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.

    """

32
    parameter_list = OrderedDict([
33
        ("temperature", (
34
            "cigale_list(minvalue=0.)",
35
            "Temperature of the dust in K.",
36
            35.
37 38
        )),
        ("beta", (
39
            "cigale_list(minvalue=0.)",
40
            "Emissivity index of the dust.",
41
            1.6
42 43
        )),
        ("alpha", (
44
            "cigale_list(minvalue=0.)",
45
            "Mid-infrared powerlaw slope.",
46
            2.
47 48 49
        ))
    ])

50 51
    def _init_code(self):
        """Build the model for a given set of parameters."""
52 53
        # To compactify the following equations we only assign them to self at
        # the end of the method
54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
        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)
Médéric Boquien's avatar
Médéric Boquien committed
69
               ** 3. / (np.exp(cst.h * c / (lambda_c * cst.k * T)) - 1.))
70 71 72 73

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

Médéric Boquien's avatar
Médéric Boquien committed
74 75 76
        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.))
77
        self.lumin_powerlaw = (conv * Npl * (self.wave / lambda_c) ** alpha *
Médéric Boquien's avatar
Médéric Boquien committed
78
                               np.exp(-(self.wave / lambda_c) ** 2.))
79 80 81 82 83 84 85 86

        # 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

87 88 89 90
        self.temperature = T
        self.beta = beta
        self.alpha = alpha

91 92 93 94 95
    def process(self, sed):
        """Add the IR re-emission contributions.

        Parameters
        ----------
96
        sed: pcigale.sed.SED object
97 98

        """
99
        if 'dust.luminosity' not in sed.info:
100
            sed.add_info('dust.luminosity', 1., True)
101
        luminosity = sed.info['dust.luminosity']
102 103

        sed.add_module(self.name, self.parameters)
104 105 106 107
        sed.add_info("dust.temperature", self.temperature)
        sed.add_info("dust.beta", self.beta)
        sed.add_info("dust.alpha", self.alpha)

108 109 110 111
        sed.add_contribution('dust.powerlaw', self.wave,
                             luminosity * self.lumin_powerlaw)
        sed.add_contribution('dust.blackbody', self.wave,
                             luminosity * self.lumin_blackbody)
Yannick Roehlly's avatar
Yannick Roehlly committed
112

113

114
# SedModule to be returned by get_module
Yannick Roehlly's avatar
Yannick Roehlly committed
115
Module = Casey2012