dl2007.py 5 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
"""
Draine and Li (2007) IR models module
=====================================

This module implements the Draine and Li (2007) infra-red models.

"""

15
from collections import OrderedDict
16

17
import numpy as np
18

19
from pcigale.data import Database
20
from . import SedModule
21
22


23
class DL2007(SedModule):
Yannick Roehlly's avatar
Yannick Roehlly committed
24
    """Draine and Li (2007) templates IR re-emission module
25
26
27
28
29
30
31
32
33
34

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

    Information added to the SED: NAME_alpha.

    """

35
    parameter_list = OrderedDict([
36
        ('qpah', (
37
38
            'cigale_list(options=0.47 & 1.12 & 1.77 & 2.50 & 3.19 & 3.90 & '
            '4.58)',
39
            "Mass fraction of PAH. Possible values are: 0.47, 1.12, 1.77, "
40
            "2.50, 3.19, 3.90, 4.58.",
41
            2.50
42
43
        )),
        ('umin', (
44
45
46
            'cigale_list(options=0.10 & 0.15 & 0.20 & 0.30 & 0.40 & 0.50 & '
            '0.70 & 0.80 & 1.00 & 1.20 & 1.50 & 2.00 & 2.50 & 3.00 & 4.00 & '
            '5.00 & 7.00 & 8.00 & 10.0 & 12.0 & 15.0 & 20.0 & 25.0)',
47
            "Minimum radiation field. Possible values are: 0.10, 0.15, 0.20, "
Médéric Boquien's avatar
Médéric Boquien committed
48
49
            "0.30, 0.40, 0.50, 0.70, 0.80, 1.00, 1.20, 1.50, 2.00, 2.50, "
            "3.00, 4.00, 5.00, 7.00, 8.00, 10.0, 12.0, 15.0, 20.0, 25.0.",
50
            1.0
51
52
        )),
        ('umax', (
53
            'cigale_list(options=1e3 & 1e4 & 1e5 & 1e6)',
54
55
            "Maximum radiation field. Possible values are: 1e3, 1e4, 1e5, "
            "1e6.",
56
            1e6
57
58
        )),
        ('gamma', (
59
            'cigale_list(minvalue=0., maxvalue=1.)',
60
61
            "Fraction illuminated from Umin to Umax. Possible values between "
            "0 and 1.",
62
            0.1
63
64
        ))
    ])
65
66
67
68

    def _init_code(self):
        """Get the model out of the database"""

69
70
71
72
        self.qpah = float(self.parameters["qpah"])
        self.umin = float(self.parameters["umin"])
        self.umax = float(self.parameters["umax"])
        self.gamma = float(self.parameters["gamma"])
73

74
75
76
77
78
        # We also compute <U>
        self.umean = (1. - self.gamma) * self.umin + \
                     self.gamma * np.log(self.umax / self.umin) / \
                     (1. / self.umin  - 1. / self.umax)

79
        with Database() as database:
80
81
82
83
            self.model_minmin = database.get_dl2007(self.qpah, self.umin,
                                                    self.umin)
            self.model_minmax = database.get_dl2007(self.qpah, self.umin,
                                                    self.umax)
84
85
86
87
88

        # The models in memory are in W/nm for 1 kg of dust. At the same time
        # we need to normalize them to 1 W here to easily scale them from the
        # power absorbed in the UV-optical. If we want to retrieve the dust
        # mass at a later point, we have to save their "emissivity" per unit
89
90
91
        # mass in W (kg of dust)¯¹, The gamma parameter does not affect the
        # fact that it is for 1 kg because it represents a mass fraction of
        # each component.
92
93
        self.emissivity = np.trapz((1.-self.gamma) * self.model_minmin.lumin +
                                   self.gamma * self.model_minmax.lumin,
94
                                   x=self.model_minmin.wave)
95

Yannick Roehlly's avatar
Yannick Roehlly committed
96
        # We want to be able to display the respective contributions of both
97
        # components, therefore we keep they separately.
98
99
        self.model_minmin.lumin *= (1. - self.gamma) / self.emissivity
        self.model_minmax.lumin *= self.gamma / self.emissivity
100
101
102
103
104
105

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

        Parameters
        ----------
106
107
        sed: pcigale.sed.SED object
        parameters: dictionary containing the parameters
108
109

        """
110
        if 'dust.luminosity' not in sed.info:
111
            sed.add_info('dust.luminosity', 1., True, unit='W')
112
        luminosity = sed.info['dust.luminosity']
113
114

        sed.add_module(self.name, self.parameters)
115
116
117
        sed.add_info('dust.qpah', self.qpah)
        sed.add_info('dust.umin', self.umin)
        sed.add_info('dust.umax', self.umax)
118
        sed.add_info('dust.umean', self.umean)
119
        sed.add_info('dust.gamma', self.gamma)
120
121
        # To compute the dust mass we simply divide the luminosity in W by the
        # emissivity in W/kg of dust.
122
123
        sed.add_info('dust.mass', luminosity / self.emissivity, True,
                     unit='solMass')
124
125
126
127
128

        sed.add_contribution('dust.Umin_Umin', self.model_minmin.wave,
                             luminosity * self.model_minmin.lumin)
        sed.add_contribution('dust.Umin_Umax', self.model_minmax.wave,
                             luminosity * self.model_minmax.lumin)
Yannick Roehlly's avatar
Yannick Roehlly committed
129

130

131
# SedModule to be returned by get_module
Yannick Roehlly's avatar
Yannick Roehlly committed
132
Module = DL2007