restframe_parameters.py 8.89 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 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 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84
# -*- coding: utf-8 -*-
# Copyright (C) 2015 Laboratoire d'Astrophysique de Marseille
# Copyright (C) 2016 Universidad de Antofagasta
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
# Author: Médéric Boquien & Denis Burgarella

"""
Module that estimates other parameters, e.g., UV slope, Lick indices, etc.
==========================================================================

This module estimates additional parameters of interest close to the
observation, e.g., the ultraviolet slope (beta), the rest-frame, any type of
indices (Lick), etc.

This module has to be called right before redshifting as such it does not take
into account the IGM absorption so it should not be used to compute fluxes in
the observed frame. For that use the param_postz module.

"""

from collections import OrderedDict
from itertools import chain

import numpy as np
from scipy.constants import c, parsec

from . import SedModule
from ..sed.utils import flux_trapz


class RestframeParam(SedModule):
    """Compute miscellaneous parameters on the full SED. This is a separate
    module as computing these quantitites in other SED modules does not always
    make much sense. This module is to be called right before the redshifting
    module. This means it does not include IGM absorbtion

    """

    parameter_list = OrderedDict([
        ("beta_calz94", (
            "boolean()",
            "UV slope measured in the same way as in Calzetti et al. (1994).",
            False
        )),
        ("D4000", (
            "boolean()",
            "D4000 break using the Balogh et al. (1999) definition.",
            False
        )),
        ("IRX", (
            "boolean()",
            "IRX computed from the GALEX FUV filter and the dust luminosity.",
            False
        )),
        ("EW_lines", (
            "string()",
            "Central wavelength of the emission lines for which to compute "
            "the equivalent width. The half-bandwidth must be indicated "
            "after the '/' sign. For instance 656.3/1.0 means oth the nebular "
            "line and the continuum are integrated over 655.3-657.3 nm.",
            "500.7/1.0 & 656.3/1.0"
        )),
        ("luminosity_filters", (
            "string()",
            "Filters for which the rest-frame luminosity will be computed. "
            "You can give several filter names separated by a & (don't use "
            "commas).",
            "FUV & V_B90"
        )),
        ("colours_filters", (
            "string()",
            "Rest-frame colours to be computed. You can give several colours "
            "separated by a & (don't use commas).",
            "FUV-NUV & NUV-r_prime"
        ))
    ])

    def calz94(self, sed):
        wl = sed.wavelength_grid
        lumin = sed.luminosity

        # Attenuated (observed) UV slopes beta as defined in Calzetti et al.
        # (1994, ApJ 429, 582, Tab. 2) that excludes the 217.5 nm bump
        # wavelength range and other spectral features
85 86 87 88 89 90 91
        if 'nebular.lines_width' in sed.info:
            key = (wl.size, sed.info['nebular.lines_width'])
        else:
            key = (wl.size, )

        if key in self.w_calz94:
            w_calz94 = self.w_calz94[key]
92 93 94 95 96 97 98 99
        else:
            calz_wl = [(126.8, 128.4), (130.9, 131.6), (134.2, 137.1),
                       (140.7, 151.5), (156.2, 158.3), (167.7, 174.0),
                       (176.0, 183.3), (186.6, 189.0), (193.0, 195.0),
                       (240.0, 258.0)]
            w_calz94 = np.where(np.any([(wl >= wlseg[0]) & (wl <= wlseg[1])
                                        for wlseg in calz_wl], axis=0))

100
            self.w_calz94[key] = w_calz94
101 102 103 104 105 106 107 108 109 110 111

        # We compute the regression directly from the covariance matrix as the
        # numpy/scipy regression routines are quite slow.
        ssxm, ssxym, _, _ = np.cov(np.log10(wl[w_calz94]),
                                   np.log10(lumin[w_calz94]),
                                   bias=1).flat

        return ssxym / ssxm

    def D4000(self, sed):
        wl = sed.wavelength_grid
112
        fnu = sed.fnu
113 114 115 116 117

        # Strength of the D_4000 break using Balogh et al. (1999, ApJ 527, 54),
        # i.e., ratio of the flux in the red continuum to that in the blue
        # continuum: Blue continuum: 385.0-395.0 nm & red continuum:
        # 400.0-410.0 nm.
118 119 120 121 122 123 124 125
        if 'nebular.lines_width' in sed.info:
            key = (wl.size, sed.info['nebular.lines_width'])
        else:
            key = (wl.size, )

        if key in self.w_D4000blue:
            w_D4000blue = self.w_D4000blue[key]
            w_D4000red = self.w_D4000red[key]
126 127 128
        else:
            w_D4000blue = np.where((wl >= 385.0) & (wl <= 395.0))
            w_D4000red = np.where((wl >= 400.0) & (wl <= 410.0))
129 130
            self.w_D4000blue[key] = w_D4000blue
            self.w_D4000red[key] = w_D4000red
131

132 133
        return (np.trapz(fnu[w_D4000red], x=wl[w_D4000red]) /
                np.trapz(fnu[w_D4000blue], x=wl[w_D4000blue]))
134 135 136 137

    def EW(self, sed):
        wl = sed.wavelength_grid

138 139 140
        key = (wl.size, sed.info['nebular.lines_width'])
        if key in self.w_lines:
            w_lines = self.w_lines[key]
141 142 143 144
        else:
            w_lines = {line: np.where((wl >= line[0]-line[1]) &
                                      (wl <= line[0]+line[1]))
                       for line in self.lines}
145
            self.w_lines[key] = w_lines
146 147 148 149 150 151 152 153 154 155

        lumin_line = np.sum([sed.get_lumin_contribution(name)
                             for name in sed.contribution_names
                             if 'nebular.lines' in name], axis=0)
        lumin_cont = sed.luminosity - lumin_line

        EW = {}
        for line in self.lines:
            w_line = w_lines[line]
            wl_line = wl[w_line]
156
            key = (wl_line.size, sed.info['nebular.lines_width'], line[0], 0.)
157 158
            EW[line] = (flux_trapz(lumin_line[w_line], wl_line, key) /
                        flux_trapz(lumin_cont[w_line], wl_line, key) *
159
                        (wl_line[-1]-wl_line[0]))
160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220

        return EW

    def _init_code(self):
        # Index of the wavelengths of interest. We use a dictionary with the
        # size of the wavelength array as a key to take into account that
        # different models may have a different sampling, for instance when
        # an AGN may be present or not.
        self.w_calz94 = {}
        self.w_D4000blue = {}
        self.w_D4000red = {}
        self.w_lines = {}

        # Extract the list of lines to compute the equivalent width
        self.lines = [item.strip() for item in
                      self.parameters["EW_lines"].split("&")
                      if item.strip() != '']
        self.lines = [(float(line.split('/')[0]),
                      float(line.split('/')[1]))
                      for line in self.lines]

        # Extract the list of filters to compute the luminosities
        self.lumin_filters = [item.strip() for item in
                              self.parameters["luminosity_filters"].split("&")
                              if item.strip() != '']

        # Extract the list of rest-frame colours
        self.colours = [item.strip().split("-") for item in
                        self.parameters["colours_filters"].split("&")
                        if item.strip() != '']

        # Extract the list of rest-frame colour filters
        self.colour_filters = list(set(chain(*self.colours)))

        # Compute the list of unique filters to compute both the luminosities
        # and the colours
        self.filters = list(set(self.lumin_filters + self.colour_filters +
                                (["FUV"] if self.parameters["IRX"] is True
                                 else [])))

        # Conversion factor to go from Fnu in mJy to Lnu in W
        self.to_lumin = 1e-29 * 4. * np.pi * (10. * parsec)**2

    def process(self, sed):
        """Computes the parameters for each model.

        Parameters
        ----------
        sed: pcigale.sed.SED object

        """
        fluxes = {filt: sed.compute_fnu(filt) for filt in self.filters}

        if self.parameters['beta_calz94']:
            sed.add_info("param.beta_calz94", self.calz94(sed))
        if self.parameters['D4000']:
            sed.add_info("param.D_4000", self.D4000(sed))
        if self.parameters['IRX']:
            sed.add_info("param.IRX", np.log10(sed.info['dust.luminosity'] /
                         (fluxes['FUV'] * self.to_lumin * c / 154e-9)))

221 222
        if 'nebular.lines_young' in sed.contribution_names:
            for line, EW in self.EW(sed).items():
223
                sed.add_info(f"param.EW({line[0]}/{line[1]})", EW)
224 225

        for filt in self.lumin_filters:
Médéric Boquien's avatar
Médéric Boquien committed
226
            sed.add_info(f"param.restframe_Lnu({filt})",
227 228 229
                         fluxes[filt] * self.to_lumin,
                         True)
        for filt1, filt2 in self.colours:
230
            sed.add_info(f"param.restframe_{filt1}-{filt2}",
231 232 233 234 235
                         2.5 * np.log10(fluxes[filt2]/fluxes[filt1]))


# SedModule to be returned by get_module
Module = RestframeParam