__init__.py 15 KB
Newer Older
Yannick Roehlly's avatar
Yannick Roehlly committed
1
# -*- coding: utf-8 -*-
2 3 4
# Copyright (C) 2012, 2013 Centre de données Astrophysiques de Marseille
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
# Author: Yannick Roehlly <yannick.roehlly@oamp.fr>
Yannick Roehlly's avatar
Yannick Roehlly committed
5

6
"""
Yannick Roehlly's avatar
Yannick Roehlly committed
7 8
This class represents a Spectral Energy Distribution (SED) as used by pcigale.
Such SED is characterised by:
Yannick Roehlly's avatar
Yannick Roehlly committed
9

Yannick Roehlly's avatar
Yannick Roehlly committed
10 11
- sfh: a tuple (time [Myr], Star Formation Rate [Msun/yr]) representing the
  Star Formation History of the galaxy.
Yannick Roehlly's avatar
Yannick Roehlly committed
12

Yannick Roehlly's avatar
Yannick Roehlly committed
13 14
- modules: a list of tuples (module name, parameter dictionary) containing all
  the pcigale modules the SED 'went through'.
Yannick Roehlly's avatar
Yannick Roehlly committed
15

Yannick Roehlly's avatar
Yannick Roehlly committed
16
- wavelength_grid: the grid of wavelengths [nm] used for the luminosities.
Yannick Roehlly's avatar
Yannick Roehlly committed
17

Yannick Roehlly's avatar
Yannick Roehlly committed
18 19
- contribution_names: the list of the names of the luminosity contributions
  making part of the SED.
Yannick Roehlly's avatar
Yannick Roehlly committed
20

Yannick Roehlly's avatar
Yannick Roehlly committed
21 22 23 24
- luminosities: a two axis numpy array containing all the luminosity density
  [W/nm] contributions to the SED. The index in the first axis corresponds to
  the contribution (in the contribution_names list) and the index of the
  second axis corresponds to the wavelength in the wavelength grid.
Yannick Roehlly's avatar
Yannick Roehlly committed
25

Yannick Roehlly's avatar
Yannick Roehlly committed
26 27
- lines: a dictionary containing the emission lines associated with the SED.
  A dictionary is used to allow the storage of various sets of lines. The
28 29 30
  lines are stored in a three axis numpy array: axis 0 is the central
  wavelength [nm], axis 1 is the line luminosity [W] and axis 2 is the line
  width [km.s-1].
31

Yannick Roehlly's avatar
Yannick Roehlly committed
32
- info: a dictionary containing various information about the SED.
Yannick Roehlly's avatar
Yannick Roehlly committed
33

Yannick Roehlly's avatar
Yannick Roehlly committed
34
- mass_proportional_info: the list of keys in the info dictionary whose value
Yannick Roehlly's avatar
Yannick Roehlly committed
35
  is proportional to the galaxy mass.
Yannick Roehlly's avatar
Yannick Roehlly committed
36

Yannick Roehlly's avatar
Yannick Roehlly committed
37
"""
Yannick Roehlly's avatar
Yannick Roehlly committed
38

Yannick Roehlly's avatar
Yannick Roehlly committed
39 40 41 42
import numpy as np
from . import utils
from scipy.constants import c
from scipy.interpolate import interp1d
Yannick Roehlly's avatar
Yannick Roehlly committed
43

44 45 46 47 48
# Time lapse used to compute the average star formation rate. We use a
# constant to keep it easily changeable for advanced user while limiting the
# number of parameters. The value is in Myr.
AV_LAPSE = 100

49

Yannick Roehlly's avatar
Yannick Roehlly committed
50 51
class SED(object):
    """Spectral Energy Distribution with associated information
Yannick Roehlly's avatar
Yannick Roehlly committed
52 53
    """

Yannick Roehlly's avatar
Yannick Roehlly committed
54 55 56 57 58 59 60 61 62 63 64 65
    def __init__(self, sfh=None):
        """Create a new SED

        Parameters
        ----------
        sfh : (numpy.array, numpy.array)
            Star Formation History: tuple of two numpy array, the first is the
            time in Myr and the second is the Star Formation Rate in Msun/yr.
            If no SFH is given, it's set to None.

        """
        self.sfh = sfh
Yannick Roehlly's avatar
Yannick Roehlly committed
66 67 68
        self.modules = []
        self.wavelength_grid = None
        self.contribution_names = []
Yannick Roehlly's avatar
Yannick Roehlly committed
69
        self.luminosities = None
70
        self.lines = {}
Yannick Roehlly's avatar
Yannick Roehlly committed
71
        self.info = {}
72
        self.mass_proportional_info = []
Yannick Roehlly's avatar
Yannick Roehlly committed
73

Yannick Roehlly's avatar
Yannick Roehlly committed
74 75 76 77 78 79 80 81 82 83 84 85 86
    @property
    def sfh(self):
        """Return a copy of the star formation history
        """
        if self._sfh is None:
            return None
        else:
            return np.copy(self._sfh)

    @sfh.setter
    def sfh(self, value):
        self._sfh = value

87 88 89 90 91 92 93 94 95
        if value:
            sfh_time, sfh_sfr = value
            sfh_age = np.max(sfh_time) - sfh_time
            self._sfh = value
            self.add_info("sfr", sfh_sfr[-1], True, True)
            self.add_info("average_sfr", np.mean(sfh_sfr[sfh_age <= AV_LAPSE]),
                          True, True)
            self.add_info("age", np.max(sfh_time), True, True)

Yannick Roehlly's avatar
Yannick Roehlly committed
96 97
    @property
    def wavelength_grid(self):
Yannick Roehlly's avatar
Yannick Roehlly committed
98
        """ Return a copy of the wavelength grid
Yannick Roehlly's avatar
Yannick Roehlly committed
99 100 101 102 103 104 105 106 107 108 109
        """
        if self._wavelength_grid is None:
            return None
        else:
            return np.copy(self._wavelength_grid)

    @wavelength_grid.setter
    def wavelength_grid(self, value):
        self._wavelength_grid = value

    @property
Yannick Roehlly's avatar
Yannick Roehlly committed
110 111
    def luminosities(self):
        """ Return a copy of the luminosity contributions
Yannick Roehlly's avatar
Yannick Roehlly committed
112
        """
Yannick Roehlly's avatar
Yannick Roehlly committed
113
        if self._luminosities is None:
Yannick Roehlly's avatar
Yannick Roehlly committed
114 115
            return None
        else:
Yannick Roehlly's avatar
Yannick Roehlly committed
116
            return np.copy(self._luminosities)
Yannick Roehlly's avatar
Yannick Roehlly committed
117

Yannick Roehlly's avatar
Yannick Roehlly committed
118 119 120
    @luminosities.setter
    def luminosities(self, value):
        self._luminosities = value
Yannick Roehlly's avatar
Yannick Roehlly committed
121 122 123

    @property
    def luminosity(self):
Yannick Roehlly's avatar
Yannick Roehlly committed
124 125
        """Total luminosity of the SED

Yannick Roehlly's avatar
Yannick Roehlly committed
126 127 128
        Return the total luminosity density vector, i.e. the sum of all the
        contributions in W/nm.
        """
Yannick Roehlly's avatar
Yannick Roehlly committed
129 130 131 132
        if self._luminosities is None:
            return None
        else:
            return self._luminosities.sum(0)
Yannick Roehlly's avatar
Yannick Roehlly committed
133

Yannick Roehlly's avatar
Yannick Roehlly committed
134
    def lambda_fnu(self, redshift=0, apply_redshift=False):
Yannick Roehlly's avatar
Yannick Roehlly committed
135 136 137 138
        """
        Return the (redshifted if asked) total Fν flux density vs wavelength
        spectrum of the SED.

Yannick Roehlly's avatar
Yannick Roehlly committed
139
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
140 141
        ----------
        redshift : float, default = 0
Yannick Roehlly's avatar
Yannick Roehlly committed
142 143
            If 0 (the default), the flux at 10 pc is computed.
        apply_redshift : boolean, default = None
144 145 146
            If true, the spectrum will be redshifted before computing the
            flux. The default is False because we generally use a specific
            module to apply the redshift.
Yannick Roehlly's avatar
Yannick Roehlly committed
147 148 149 150 151 152 153

        Returns
        -------
        wavelength, f_nu : tuple of array of floats
            The wavelength is in nm and the Fν is in mJy.

        """
Yannick Roehlly's avatar
Yannick Roehlly committed
154 155
        wavelength = self.wavelength_grid
        luminosity = self.luminosity
Yannick Roehlly's avatar
Yannick Roehlly committed
156

Yannick Roehlly's avatar
Yannick Roehlly committed
157 158 159
        if apply_redshift:
            wavelength, luminosity = utils.redshift_spectrum(
                wavelength, luminosity, redshift)
Yannick Roehlly's avatar
Yannick Roehlly committed
160

Yannick Roehlly's avatar
Yannick Roehlly committed
161 162 163 164 165
        # Fλ flux density in W/m²/nm
        f_lambda = utils.luminosity_to_flux(luminosity, redshift)

        # Fν flux density in Jy
        f_nu = utils.lambda_flambda_to_fnu(wavelength, f_lambda)
Yannick Roehlly's avatar
Yannick Roehlly committed
166 167 168

        return wavelength, f_nu

169
    def add_info(self, key, value, mass_proportional=False, force=False):
Yannick Roehlly's avatar
Yannick Roehlly committed
170 171 172
        """
        Add a key / value to the information dictionary

Yannick Roehlly's avatar
Yannick Roehlly committed
173
        If the key is present in the dictionary, it will raise an exception.
Yannick Roehlly's avatar
Yannick Roehlly committed
174
        Use this method (instead of direct value assignment ) to avoid
175
        overriding an already present information.
Yannick Roehlly's avatar
Yannick Roehlly committed
176

Yannick Roehlly's avatar
Yannick Roehlly committed
177
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
178 179 180 181 182
        ----------
        key : any immutable
           The key used to retrieve the information.
        value : anything
           The information.
183 184 185
        mass_proportional : boolean
           If True, the added variable is set as proportional to the
           mass.
186
        force : boolean
187
           If false (default), adding a key that already exists in the info
188 189
           dictionary will raise an error. If true, doing this will update
           the associated value.
Yannick Roehlly's avatar
Yannick Roehlly committed
190 191

        """
192
        if (key not in self.info) or force:
Yannick Roehlly's avatar
Yannick Roehlly committed
193
            self.info[key] = value
194 195
            if mass_proportional:
                self.mass_proportional_info.append(key)
Yannick Roehlly's avatar
Yannick Roehlly committed
196
        else:
197
            raise KeyError("The information %s is already present "
Yannick Roehlly's avatar
Yannick Roehlly committed
198 199
                           "in the SED. " % key)

200 201 202
    def add_module(self, module_name, module_conf):
        """Add a new module information to the SED.

Yannick Roehlly's avatar
Yannick Roehlly committed
203
        Parameters
204 205 206 207 208
        ----------
        module_name : string
            Name of the module. This name can be suffixed with anything
            using a dot.
        module_conf : dictionary
Yannick Roehlly's avatar
Yannick Roehlly committed
209
            Dictionary containing the module parameters.
210

Yannick Roehlly's avatar
Yannick Roehlly committed
211
        TODO: Complete the parameter dictionary with the default values from
212 213
              the module if they are not present.

Yannick Roehlly's avatar
Yannick Roehlly committed
214
        """
215
        self.modules.append((module_name, module_conf))
Yannick Roehlly's avatar
Yannick Roehlly committed
216

217 218 219 220
    def add_contribution(self, contribution_name, results_wavelengths,
                         results_lumin):
        """
        Add a new luminosity contribution to the SED.
Yannick Roehlly's avatar
Yannick Roehlly committed
221 222 223 224 225 226 227 228

        The luminosity contribution of the module is added to the contribution
        table doing an interpolation between the current wavelength grid and
        the grid of the module contribution. During the interpolation,
        everything that is outside of the concerned wavelength domain has its
        luminosity set to 0. Also, the name of the contribution is added to
        the contribution names array.

Yannick Roehlly's avatar
Yannick Roehlly committed
229
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246
        ----------
        contribution_name : string
            Name of the contribution added. This name is used to retrieve the
            luminosity contribution and allows one module to add more than
            one contribution.

        results_wavelengths : array of floats
            The vector of the wavelengths of the module results (in nm).

        results_lumin : array of floats
            The vector of the Lλ luminosities (in W/nm) of the module results.

        """
        self.contribution_names.append(contribution_name)

        # If the SED luminosity table is empty, then there is nothing to
        # compute.
Yannick Roehlly's avatar
Yannick Roehlly committed
247
        if self.luminosities is None:
Yannick Roehlly's avatar
Yannick Roehlly committed
248
            self.wavelength_grid = np.copy(results_wavelengths)
Yannick Roehlly's avatar
Yannick Roehlly committed
249
            self.luminosities = np.copy(results_lumin)
Yannick Roehlly's avatar
Yannick Roehlly committed
250 251 252 253 254 255 256
        else:
            # Compute the new wavelength grid for the spectrum.
            new_wavelength_grid = utils.best_grid(results_wavelengths,
                                                  self.wavelength_grid)

            # Interpolate each luminosity component to the new wavelength grid
            # setting everything outside the wavelength domain to 0.
Yannick Roehlly's avatar
Yannick Roehlly committed
257 258 259 260 261 262 263
            new_luminosities = interp1d(self.wavelength_grid,
                                        self.luminosities,
                                        bounds_error=False,
                                        fill_value=0.)(new_wavelength_grid)

            # Interpolate the added luminosity array to the new wavelength
            # grid
264 265 266 267
            interp_lumin = interp1d(results_wavelengths,
                                    results_lumin,
                                    bounds_error=False,
                                    fill_value=0)(new_wavelength_grid)
Yannick Roehlly's avatar
Yannick Roehlly committed
268 269

            self.wavelength_grid = new_wavelength_grid
Yannick Roehlly's avatar
Yannick Roehlly committed
270
            self.luminosities = np.vstack((new_luminosities, interp_lumin))
Yannick Roehlly's avatar
Yannick Roehlly committed
271

272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297
    def add_lines(self, set_name, wavelengths, luminosities, widths):
        """Add a set of spectral lines to the SED.

        Parameters
        ----------
        set_name : string
            Name of the set of lines.
        wavelengths : list-like of floats
            The central wavelengths of the lines in nm.
        luminosities : list-like of floats
            The total luminosity in the spectral lines in W.
        widths : list-like of floats
            The widths of the lines in km.s-1.

        """
        wavelengths = np.array(wavelengths, dtype=float)
        luminosities = np.array(luminosities, dtype=float)
        widths = np.array(widths, dtype=float)

        if set_name in self.lines:
            raise KeyError("The line set {} is all ready present in the "
                           "SED.".format(set_name))
        else:
            self.lines[set_name] = np.vstack(
                (wavelengths, luminosities, widths))

Yannick Roehlly's avatar
Yannick Roehlly committed
298 299 300 301 302 303
    def get_lumin_contribution(self, name):
        """Get the luminosity vector of a given contribution

        If the name of the contribution is not unique in the SED, the flux of
        the last one is returned.

Yannick Roehlly's avatar
Yannick Roehlly committed
304
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
305 306 307 308 309 310 311 312 313 314 315 316 317 318
        ----------
        name : string
            Name of the contribution

        Returns
        -------
        luminosities : array of floats
            Vector of the luminosity density contribution based on the SED
            wavelength grid.

        """
        # Find the index of the _last_ name element
        idx = (len(self.contribution_names) - 1
               - self.contribution_names[::-1].index(name))
Yannick Roehlly's avatar
Yannick Roehlly committed
319
        return self.luminosities[idx]
Yannick Roehlly's avatar
Yannick Roehlly committed
320

321
    def compute_fnu(self, transmission, lambda_eff,
322 323
                    redshift=0, apply_redshift=False,
                    add_line_fluxes=True):
Yannick Roehlly's avatar
Yannick Roehlly committed
324 325 326 327 328 329 330 331 332 333 334 335
        """
        Compute the Fν flux density corresponding the filter which
        transmission is given.

        As the SED stores the Lλ luminosity density, we first compute the Fλ
        flux density. Fλ is the integration of the Lλ luminosity multiplied by
        the filter transmission, normalised to this transmission and corrected
        by the luminosity distance of the source. This is done by the
        pcigale.sed.utils.luminosity_to_flux function.

        Fλ = luminosity_to_flux( integ( LλT(λ)dλ ) / integ( T(λ)dλ ) )

Yannick Roehlly's avatar
Yannick Roehlly committed
336
        Fλ is in W/m²/nm. At redshift 0, the flux is computed at 10 pc. Then,
Yannick Roehlly's avatar
Yannick Roehlly committed
337 338 339 340 341 342 343 344 345
        to compute Fν, we make the approximation:

        Fν = λeff / c . λeff . Fλ

        Fν is computed in W/m²/Hz and then converted to mJy.

        If the SED spectrum does not cover all the filter response table,
        -99 is returned.

Yannick Roehlly's avatar
Yannick Roehlly committed
346
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
347 348 349 350 351 352 353 354 355
        ----------
        transmission : 2D array of floats
            A numpy 2D array containing the filter response profile
            wavelength[nm] vs transmission).

        lambda_eff : float
            Effective wavelength of the filter in nm.

        redshift : float
Yannick Roehlly's avatar
Yannick Roehlly committed
356
            The redshift of the galaxy. If 0, the flux is computed at 10 pc.
Yannick Roehlly's avatar
Yannick Roehlly committed
357

Yannick Roehlly's avatar
Yannick Roehlly committed
358
        apply_redshift : boolean
359 360 361 362
            If true, the spectrum will be redshifted before computing the
            flux. The default is False because we generally use a specific
            module to apply the redshift.

363 364 365 366
        add_line_flux : boolean
            If true (default), the flux coming from the spectral lines will be
            taken into account.

Yannick Roehlly's avatar
Yannick Roehlly committed
367 368 369 370 371 372
        Return
        ------
        fnu : float
            The integrated Fν density in mJy.
        """

373
        # Filter limits
Yannick Roehlly's avatar
Yannick Roehlly committed
374 375 376
        lambda_min = min(transmission[0])
        lambda_max = max(transmission[0])

377 378 379 380 381 382 383
        wavelength = self.wavelength_grid
        l_lambda = self.luminosity
        if apply_redshift:
            wavelength, l_lambda = utils.redshift_lambda_l_lambda(
                (wavelength, l_lambda), redshift)

        # Test if the spectrum cover all the filter extend
Yannick Roehlly's avatar
Yannick Roehlly committed
384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407
        if ((min(self.wavelength_grid) > lambda_min) or
                (max(self.wavelength_grid) < lambda_max)):
            f_nu = -99.

        else:
            # We regrid both spectrum and filter to the best wavelength grid
            # to avoid interpolating a high wavelength density curve to a low
            # density one. Also, we limit the work wavelength domain to the
            # filter one, taking care the presence of λmin and λman in the
            # used wavelength grid.
            wavelength_r = utils.best_grid(wavelength, transmission[0])
            if lambda_min not in wavelength_r:
                wavelength_r.append(lambda_min)
            if lambda_max not in wavelength_r:
                wavelength_r.append(lambda_max)
            wavelength_r.sort()
            wavelength_r = wavelength_r[wavelength_r <= lambda_max]
            wavelength_r = wavelength_r[wavelength_r >= lambda_min]

            l_lambda_r = np.interp(wavelength_r, wavelength, l_lambda)
            transmission_r = np.interp(wavelength_r, transmission[0],
                                       transmission[1])

            # TODO: Can we avoid to normalise as the filter transmission is
408
            # already normalised?
Yannick Roehlly's avatar
Yannick Roehlly committed
409 410 411 412 413 414
            f_lambda = utils.luminosity_to_flux(
                (np.trapz(transmission_r * l_lambda_r, wavelength_r) /
                 np.trapz(transmission_r, wavelength_r)),
                redshift
            )

415
            # Add the Fλ fluxes from the spectral lines.
416
            # TODO write the code
417

Yannick Roehlly's avatar
Yannick Roehlly committed
418 419 420 421 422 423 424
            # Fν in W/m²/Hz. The 1.e-9 factor is because λ is in nm.
            f_nu = lambda_eff * f_lambda * lambda_eff * 1.e-9 / c

            # Conversion from W/m²/Hz to mJy
            f_nu *= 1.e+29

        return f_nu