schreiber2016.py 4.13 KB
Newer Older
Laure Ciesla's avatar
Laure Ciesla committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
# -*- coding: utf-8 -*-
# Copyright (C) 2013 Centre de données Astrophysiques de Marseille
# Copyright (C) 2013 Institute of Astronomy, University of Cambridge
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
# Author: Laure Ciesla

"""
Schreiber (2016) IR models module
=====================================

This module implements the Schreiber et al. (2016) infra-red models.

"""

from collections import OrderedDict
import numpy as np
from pcigale.data import Database
18
from . import SedModule
Laure Ciesla's avatar
Laure Ciesla committed
19 20


21
class Schreiber2016(SedModule):
Laure Ciesla's avatar
Laure Ciesla committed
22 23 24 25 26 27 28 29 30 31 32
    """Schreiber et al. (2016) templates IR re-emission module

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

    """

    parameter_list = OrderedDict([
        ('tdust', (
33 34 35 36 37
            'cigale_list(options=15. & 16. & 17. & 18. & 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.)',
Laure Ciesla's avatar
Laure Ciesla committed
38 39 40 41 42
            "Dust temperature. "
            "Between 15 and 60K, with 1K step.",
            20.
        )),
        ('fpah', (
43
            'cigale_list(minvalue=0., maxvalue=1.)',
Laure Ciesla's avatar
Laure Ciesla committed
44 45
            "Mass fraction of PAH. "
            "Between 0 and 1.",
46
            0.05
Laure Ciesla's avatar
Laure Ciesla committed
47 48 49 50 51 52 53 54 55 56 57
        ))
    ])

    out_parameter_list = OrderedDict([
        ('tdust', 'Dust temperature'),
        ('fpah', 'Mass fraction of PAH')
    ])

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

58 59
        self.tdust = float(self.parameters["tdust"])
        self.fpah = float(self.parameters["fpah"])
Laure Ciesla's avatar
Laure Ciesla committed
60
        with Database() as database:
61 62
            self.model_dust = database.get_schreiber2016(0, self.tdust)
            self.model_pah = database.get_schreiber2016(1, self.tdust)
Laure Ciesla's avatar
Laure Ciesla committed
63 64 65 66 67 68 69 70

        # The models in memory are in W/nm/kg. 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
        # mass in W kg¯¹, The gamma parameter does not affect the fact that it
        # is for 1 kg because it represents a mass fraction of each component.

71 72
        self.emissivity = np.trapz((1. - self.fpah) * self.model_dust.lumin +
                                   self.fpah * self.model_pah.lumin,
Laure Ciesla's avatar
Laure Ciesla committed
73
                                   x=self.model_dust.wave)
74

Laure Ciesla's avatar
Laure Ciesla committed
75 76
        # We want to be able to display the respective contributions of both
        # components, therefore we keep they separately.
77 78
        self.model_dust.lumin *= (1. - self.fpah) / self.emissivity
        self.model_pah.lumin *= self.fpah / self.emissivity
Laure Ciesla's avatar
Laure Ciesla committed
79 80 81 82 83 84 85 86 87 88 89 90 91 92 93

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

        Parameters
        ----------
        sed: pcigale.sed.SED object
        parameters: dictionary containing the parameters

        """
        if 'dust.luminosity' not in sed.info.keys():
            sed.add_info('dust.luminosity', 1., True)
        luminosity = sed.info['dust.luminosity']

        sed.add_module(self.name, self.parameters)
94 95
        sed.add_info('dust.tdust', self.tdust)
        sed.add_info('dust.fpah', self.fpah)
Laure Ciesla's avatar
Laure Ciesla committed
96 97 98 99 100 101 102 103 104 105 106
        # To compute the dust mass we simply divide the luminosity by the
        # emissivity and then by the expected MH/Mdust as the emissivity was
        # computed for 1 kg of H. Note that we take 100 here but it should vary
        # with the exact model. Fix that later. Maybe directly in the database.
        sed.add_info('dust.mass', luminosity / self.emissivity, True)

        sed.add_contribution('dust.dust_continuum', self.model_dust.wave,
                             luminosity * self.model_dust.lumin)
        sed.add_contribution('dust.pah', self.model_pah.wave,
                             luminosity * self.model_pah.lumin)

107

108
# SedModule to be returned by get_module
Laure Ciesla's avatar
Laure Ciesla committed
109
Module = Schreiber2016