nebular.py 8.47 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 98
        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'])
        self.emission = bool(self.parameters["emission"])
99

100
        if self.fesc < 0. or self.fesc > 1:
101 102
            raise Exception("Escape fraction must be between 0 and 1")

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

107
        if self.fesc + self.fdust > 1:
108 109
            raise Exception("Escape fraction+f_dust>1")

110 111
        if self.emission:
            with Database() as db:
112
                metallicities = db.get_nebular_continuum_parameters()['metallicity']
113
                self.lines_template = {m: db.get_nebular_lines(m, self.logU)
114
                                    for m in metallicities}
115
                self.cont_template = {m: db.get_nebular_continuum(m, self.logU)
116 117 118 119 120 121
                                    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}
122 123 124 125

            for lines in self.lines_template.values():
                new_wave = np.array([])
                for line_wave in lines.wave:
126
                    width = line_wave * self.lines_width * 1e3 / cst.c
127 128 129 130 131 132 133
                    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):
134
                    width = line_wave * self.lines_width * 1e3 / cst.c
135 136 137 138 139 140 141 142 143
                    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))
144
            alpha_B = 2.58e-19  # Ferland 1980, m³ s¯¹
145 146 147 148
            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))

149
            self.corr = k
150 151 152
        self.idx_Ly_break = None
        self.absorbed_old = None
        self.absorbed_young = None
153 154 155 156 157 158

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

        Parameters
        ----------
159 160
        sed: pcigale.sed.SED object
        parameters: dictionary containing the parameters
161 162

        """
163 164 165 166 167
        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)

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

175
        sed.add_module(self.name, self.parameters)
176 177
        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
178
        sed.add_info('dust.luminosity', (sed.info['stellar.lum_ly_young'] +
179 180 181 182 183 184 185 186 187 188
                     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']
189 190 191 192 193
            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]
194

195 196
            sed.add_info('nebular.lines_width', self.lines_width)
            sed.add_info('nebular.logU', self.logU)
197

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

204
            sed.add_contribution('nebular.lines_old', lines.wave,
205
                                 lines.ratio * NLy_old * self.corr)
206
            sed.add_contribution('nebular.lines_young', lines.wave,
207
                                 lines.ratio * NLy_young * self.corr)
208 209

            sed.add_contribution('nebular.continuum_old', cont.wave,
210
                                 cont.lumin * NLy_old * self.corr)
211
            sed.add_contribution('nebular.continuum_young', cont.wave,
212
                                 cont.lumin * NLy_young * self.corr)
213

214

215
# SedModule to be returned by get_module
216
Module = NebularEmission