nebular.py 8.61 KB
Newer Older
1
2
3
4
5
# -*- coding: utf-8 -*-
# Copyright (C) 2014 University of Cambridge
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
# Author: Médéric Boquien <mboquien@ast.cam.ac.uk>

6
from collections import OrderedDict
7
from copy import deepcopy
8

9
10
import numpy as np
import scipy.constants as cst
11
12

from pcigale.data import Database
13
from . import SedModule
14

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
default_lines = ['Ly-alpha',
                 'CII-133.5',
                 'SiIV-139.7',
                 'CIV-154.9',
                 'HeII-164.0',
                 'OIII-166.5',
                 'CIII-190.9',
                 'CII-232.6',
                 'MgII-279.8',
                 'OII-372.7',
                 'H-10',
                 'H-9',
                 'NeIII-386.9',
                 'HeI-388.9',
                 'H-epsilon',
                 'SII-407.0',
                 'H-delta',
                 'H-gamma',
                 'H-beta',
                 'OIII-495.9',
                 'OIII-500.7',
                 'OI-630.0',
                 'NII-654.8',
                 'H-alpha',
                 'NII-658.4',
                 'SII-671.6',
                 'SII-673.1'
                 ]
43

44
class NebularEmission(SedModule):
45
    """
46
47
    Module computing the nebular emission from the ultraviolet to the
    near-infrared. It includes both the nebular lines and the nubular
48
49
    continuum (optional). It takes into account the escape fraction and the
    absorption by dust.
50
51
52
53

    Given the number of Lyman continuum photons, we compute the Hβ line
    luminosity. We then compute the other lines using the
    metallicity-dependent templates that provide the ratio between individual
54
55
    lines and Hβ. The nebular continuum is scaled directly from the number of
    ionizing photons.
56
57
58

    """

59
    parameter_list = OrderedDict([
60
        ('logU', (
61
62
63
64
            'cigale_list(options=-4.0 & -3.9 & -3.8 & -3.7 & -3.6 & -3.5 & '
            '-3.4 & -3.3 & -3.2 & -3.1 & -3.0 & -2.9 & -2.8 & -2.7 & -2.6 & '
            '-2.5 & -2.4 & -2.3 & -2.2 & -2.1 & -2.0 & -1.9 & -1.8 & -1.7 & '
            '-1.6 & -1.5 & -1.4 & -1.3 & -1.2 & -1.1 & -1.0)',
65
            "Ionisation parameter",
66
67
            -2.
        )),
68
        ('f_esc', (
69
            'cigale_list(minvalue=0., maxvalue=1.)',
70
71
72
73
            "Fraction of Lyman continuum photons escaping the galaxy",
            0.
        )),
        ('f_dust', (
74
            'cigale_list(minvalue=0., maxvalue=1.)',
75
76
77
            "Fraction of Lyman continuum photons absorbed by dust",
            0.
        )),
78
        ('lines_width', (
79
            'cigale_list(minvalue=0.)',
80
            "Line width in km/s",
81
            300.
82
83
        )),
        ('emission', (
84
            'boolean()',
85
86
            "Include nebular emission.",
            True
87
88
89
90
91
        ))
    ])

    def _init_code(self):
        """Get the nebular emission lines out of the database and resample
92
           them to see the line profile. Compute scaling coefficients.
93
        """
94
95
96
97
        self.logU = float(self.parameters['logU'])
        self.fesc = float(self.parameters['f_esc'])
        self.fdust = float(self.parameters['f_dust'])
        self.lines_width = float(self.parameters['lines_width'])
98
99
100
101
        if type(self.parameters["emission"]) is str:
            self.emission = self.parameters["emission"].lower() == 'true'
        else:
            self.emission = bool(self.parameters["emission"])
102

103
        if self.fesc < 0. or self.fesc > 1:
104
105
            raise Exception("Escape fraction must be between 0 and 1")

106
        if self.fdust < 0 or self.fdust > 1:
107
108
109
            raise Exception("Fraction of lyman photons absorbed by dust must "
                            "be between 0 and 1")

110
        if self.fesc + self.fdust > 1:
111
112
            raise Exception("Escape fraction+f_dust>1")

113
114
        if self.emission:
            with Database() as db:
115
                metallicities = db.get_nebular_continuum_parameters()['metallicity']
116
                self.lines_template = {m: db.get_nebular_lines(m, self.logU)
117
                                    for m in metallicities}
118
                self.cont_template = {m: db.get_nebular_continuum(m, self.logU)
119
120
121
122
123
124
                                    for m in metallicities}

            self.linesdict = {m: dict(zip(self.lines_template[m].name,
                                          zip(self.lines_template[m].wave,
                                              self.lines_template[m].ratio)))
                              for m in metallicities}
125
126
127
128

            for lines in self.lines_template.values():
                new_wave = np.array([])
                for line_wave in lines.wave:
129
                    width = line_wave * self.lines_width * 1e3 / cst.c
130
131
132
133
134
135
136
                    new_wave = np.concatenate((new_wave,
                                            np.linspace(line_wave - 3. * width,
                                                        line_wave + 3. * width,
                                                        9)))
                new_wave.sort()
                new_flux = np.zeros_like(new_wave)
                for line_flux, line_wave in zip(lines.ratio, lines.wave):
137
                    width = line_wave * self.lines_width * 1e3 / cst.c
138
139
140
141
142
143
144
145
146
                    new_flux += (line_flux * np.exp(- 4. * np.log(2.) *
                                (new_wave - line_wave) ** 2. / (width * width)) /
                                (width * np.sqrt(np.pi / np.log(2.)) / 2.))
                lines.wave = new_wave
                lines.ratio = new_flux

            # To take into acount the escape fraction and the fraction of Lyman
            # continuum photons absorbed by dust we correct by a factor
            # k=(1-fesc-fdust)/(1+(α1/αβ)*(fesc+fdust))
147
            alpha_B = 2.58e-19  # Ferland 1980, m³ s¯¹
148
149
150
151
            alpha_1 = 1.54e-19  # αA-αB, Ferland 1980, m³ s¯¹
            k = (1. - self.fesc - self.fdust) / (1. + alpha_1 / alpha_B * (
                self.fesc + self.fdust))

152
            self.corr = k
153
154
155
        self.idx_Ly_break = None
        self.absorbed_old = None
        self.absorbed_young = None
156
157
158
159
160
161

    def process(self, sed):
        """Add the nebular emission lines

        Parameters
        ----------
162
163
        sed: pcigale.sed.SED object
        parameters: dictionary containing the parameters
164
165

        """
166
167
168
169
170
        if self.idx_Ly_break is None:
            self.idx_Ly_break = np.searchsorted(sed.wavelength_grid, 91.2)
            self.absorbed_old = np.zeros(sed.wavelength_grid.size)
            self.absorbed_young = np.zeros(sed.wavelength_grid.size)

171
        self.absorbed_old[:self.idx_Ly_break] = -(
172
173
            sed.get_lumin_contribution('stellar.old')[:self.idx_Ly_break] *
            (1. - self.fesc))
174
        self.absorbed_young[:self.idx_Ly_break] = -(
175
176
            sed.get_lumin_contribution('stellar.young')[:self.idx_Ly_break] *
            (1. - self.fesc))
177

178
        sed.add_module(self.name, self.parameters)
179
180
        sed.add_info('nebular.f_esc', self.fesc)
        sed.add_info('nebular.f_dust', self.fdust)
Médéric Boquien's avatar
Médéric Boquien committed
181
        sed.add_info('dust.luminosity', (sed.info['stellar.lum_ly_young'] +
182
183
184
185
186
187
188
189
190
191
                     sed.info['stellar.lum_ly_old']) * self.fdust, True)

        sed.add_contribution('nebular.absorption_old', sed.wavelength_grid,
                             self.absorbed_old)
        sed.add_contribution('nebular.absorption_young', sed.wavelength_grid,
                             self.absorbed_young)

        if self.emission:
            NLy_old = sed.info['stellar.n_ly_old']
            NLy_young = sed.info['stellar.n_ly_young']
192
193
194
195
196
            NLy_tot = NLy_old + NLy_young
            metallicity = sed.info['stellar.metallicity']
            lines = self.lines_template[metallicity]
            linesdict = self.linesdict[metallicity]
            cont = self.cont_template[metallicity]
197

198
199
            sed.add_info('nebular.lines_width', self.lines_width)
            sed.add_info('nebular.logU', self.logU)
200

201
202
203
204
205
206
            for line in default_lines:
                wave, ratio = linesdict[line]
                sed.lines[line] = (wave,
                                   ratio * NLy_old * self.corr,
                                   ratio * NLy_young * self.corr)

207
            sed.add_contribution('nebular.lines_old', lines.wave,
208
                                 lines.ratio * NLy_old * self.corr)
209
            sed.add_contribution('nebular.lines_young', lines.wave,
210
                                 lines.ratio * NLy_young * self.corr)
211
212

            sed.add_contribution('nebular.continuum_old', cont.wave,
213
                                 cont.lumin * NLy_old * self.corr)
214
            sed.add_contribution('nebular.continuum_young', cont.wave,
215
                                 cont.lumin * NLy_young * self.corr)
216

217

218
# SedModule to be returned by get_module
219
Module = NebularEmission