__init__.py 13.4 KB
Newer Older
Yannick Roehlly's avatar
Yannick Roehlly committed
1
# -*- coding: utf-8 -*-
2
3
# Copyright (C) 2012, 2013 Centre de données Astrophysiques de Marseille
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
Yannick Roehlly's avatar
Yannick Roehlly committed
4
# Author: Yannick Roehlly
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
- info: a dictionary containing various information about the SED.
Yannick Roehlly's avatar
Yannick Roehlly committed
27

28
- mass_proportional_info: the set of keys in the info dictionary whose value
Yannick Roehlly's avatar
Yannick Roehlly committed
29
  is proportional to the galaxy mass.
Yannick Roehlly's avatar
Yannick Roehlly committed
30

Yannick Roehlly's avatar
Yannick Roehlly committed
31
"""
Yannick Roehlly's avatar
Yannick Roehlly committed
32

Yannick Roehlly's avatar
Yannick Roehlly committed
33
import numpy as np
34
35
from scipy.constants import c, parsec

Yannick Roehlly's avatar
Yannick Roehlly committed
36
from . import utils
37
from .io.vo import save_sed_to_vo
38
39
from ..data import Database

Yannick Roehlly's avatar
Yannick Roehlly committed
40

41

42
43
44
45
46
# 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

47

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

52
53
54
55
    # We declare the filters cache here as to be efficient it needs to be
    # shared between different objects.
    cache_filters = {}

Yannick Roehlly's avatar
Yannick Roehlly committed
56
57
58
59
60
    def __init__(self, sfh=None):
        """Create a new SED

        Parameters
        ----------
61
        sfh: (numpy.array, numpy.array)
Yannick Roehlly's avatar
Yannick Roehlly committed
62
63
64
65
66
67
            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
68
69
70
        self.modules = []
        self.wavelength_grid = None
        self.contribution_names = []
71
        self.luminosity = None
Yannick Roehlly's avatar
Yannick Roehlly committed
72
        self.luminosities = None
73
        self.info = dict()
74
        self.mass_proportional_info = set()
Yannick Roehlly's avatar
Yannick Roehlly committed
75

Yannick Roehlly's avatar
Yannick Roehlly committed
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):
Yannick Roehlly's avatar
Yannick Roehlly committed
87
88
89
90

        # The SFH can be set multiple times. Maybe it's better to make is
        # settable only once and then provide an update_sfh method for when
        # it's needed.
Yannick Roehlly's avatar
Yannick Roehlly committed
91
92
        self._sfh = value

93
94
95
        if value:
            sfh_time, sfh_sfr = value
            self._sfh = value
Yannick Roehlly's avatar
Yannick Roehlly committed
96
97
98
99
100
101
            self.add_info("sfh.sfr", sfh_sfr[-1], True, force=True)
            self.add_info("sfh.sfr10Myrs", np.mean(sfh_sfr[-10:]), True,
                          force=True)
            self.add_info("sfh.sfr100Myrs", np.mean(sfh_sfr[-100:]), True,
                          force=True)
            self.add_info("sfh.age", sfh_time[-1], False, force=True)
102

103
104
105
    @property
    def fnu(self):
        """Total Fν flux density of the SED
Yannick Roehlly's avatar
Yannick Roehlly committed
106

107
108
        Return the total Fν density vector, i.e the total luminosity converted
        to Fν flux in mJy.
Yannick Roehlly's avatar
Yannick Roehlly committed
109
110
        """

Yannick Roehlly's avatar
Yannick Roehlly committed
111
        # Fλ flux density in W/m²/nm
112
113
114
115
116
117
        if 'universe.luminosity_distance' in self.info:
            f_lambda = utils.luminosity_to_flux(self.luminosity,
                                                self.info
                                                ['universe.luminosity_distance'])
        else:
            f_lambda = utils.luminosity_to_flux(self.luminosity, 10. * parsec)
Yannick Roehlly's avatar
Yannick Roehlly committed
118

119
120
        # Fν flux density in mJy
        f_nu = utils.lambda_flambda_to_fnu(self.wavelength_grid, f_lambda)
Yannick Roehlly's avatar
Yannick Roehlly committed
121

122
        return f_nu
Yannick Roehlly's avatar
Yannick Roehlly committed
123

124
    def add_info(self, key, value, mass_proportional=False, force=False):
Yannick Roehlly's avatar
Yannick Roehlly committed
125
126
127
        """
        Add a key / value to the information dictionary

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

Yannick Roehlly's avatar
Yannick Roehlly committed
132
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
133
        ----------
134
        key: any immutable
Yannick Roehlly's avatar
Yannick Roehlly committed
135
           The key used to retrieve the information.
136
        value: anything
Yannick Roehlly's avatar
Yannick Roehlly committed
137
           The information.
138
        mass_proportional: boolean
139
140
           If True, the added variable is set as proportional to the
           mass.
141
        force: boolean
142
           If false (default), adding a key that already exists in the info
143
144
           dictionary will raise an error. If true, doing this will update
           the associated value.
Yannick Roehlly's avatar
Yannick Roehlly committed
145
146

        """
147
        if (key not in self.info) or force:
Yannick Roehlly's avatar
Yannick Roehlly committed
148
            self.info[key] = value
149
            if mass_proportional:
150
                self.mass_proportional_info.add(key)
Yannick Roehlly's avatar
Yannick Roehlly committed
151
        else:
152
            raise KeyError("The information %s is already present "
Yannick Roehlly's avatar
Yannick Roehlly committed
153
154
                           "in the SED. " % key)

155
156
157
    def add_module(self, module_name, module_conf):
        """Add a new module information to the SED.

Yannick Roehlly's avatar
Yannick Roehlly committed
158
        Parameters
159
        ----------
160
        module_name: string
161
162
            Name of the module. This name can be suffixed with anything
            using a dot.
163
        module_conf: dictionary
Yannick Roehlly's avatar
Yannick Roehlly committed
164
            Dictionary containing the module parameters.
165

Yannick Roehlly's avatar
Yannick Roehlly committed
166
        TODO: Complete the parameter dictionary with the default values from
167
168
              the module if they are not present.

Yannick Roehlly's avatar
Yannick Roehlly committed
169
        """
170
        self.modules.append((module_name, module_conf))
Yannick Roehlly's avatar
Yannick Roehlly committed
171

172
173
174
175
    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
176
177
178
179
180
181
182
183

        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
184
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
185
        ----------
186
        contribution_name: string
Yannick Roehlly's avatar
Yannick Roehlly committed
187
188
189
190
            Name of the contribution added. This name is used to retrieve the
            luminosity contribution and allows one module to add more than
            one contribution.

191
        results_wavelengths: array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
192
193
            The vector of the wavelengths of the module results (in nm).

194
        results_lumin: array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
195
196
197
198
199
200
201
            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.
202
        if self.luminosity is None:
203
204
205
            self.wavelength_grid = results_wavelengths.copy()
            self.luminosity = results_lumin.copy()
            self.luminosities = results_lumin.copy()
Yannick Roehlly's avatar
Yannick Roehlly committed
206
        else:
207
208
            # If the added luminosity contribution changes the SED wavelength
            # grid, we interpolate everything on a common wavelength grid.
209
210
211
212
            if (results_wavelengths.size != self.wavelength_grid.size or
                    not np.all(results_wavelengths == self.wavelength_grid)):
                # Interpolate each luminosity component to the new wavelength
                # grid setting everything outside the wavelength domain to 0.
213
214
                self.wavelength_grid, self.luminosities = \
                    utils.interpolate_lumin(self.wavelength_grid,
215
                                            self.luminosities,
216
217
218
                                            results_wavelengths,
                                            results_lumin)

219
                self.luminosity = self.luminosities.sum(0)
220
221
222
            else:
                self.luminosities = np.vstack((self.luminosities,
                                               results_lumin))
223
                self.luminosity += results_lumin
Yannick Roehlly's avatar
Yannick Roehlly committed
224
225
226
227
228
229
230

    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
231
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
232
        ----------
233
        name: string
Yannick Roehlly's avatar
Yannick Roehlly committed
234
235
236
237
            Name of the contribution

        Returns
        -------
238
        luminosities: array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
239
240
241
242
243
            Vector of the luminosity density contribution based on the SED
            wavelength grid.

        """
        # Find the index of the _last_ name element
Médéric Boquien's avatar
Médéric Boquien committed
244
245
        idx = (len(self.contribution_names) - 1 -
               self.contribution_names[::-1].index(name))
Yannick Roehlly's avatar
Yannick Roehlly committed
246
        return self.luminosities[idx]
Yannick Roehlly's avatar
Yannick Roehlly committed
247

248
    def compute_fnu(self, filter_name):
Yannick Roehlly's avatar
Yannick Roehlly committed
249
        """
250
        Compute the Fν flux density in a given filter
Yannick Roehlly's avatar
Yannick Roehlly committed
251
252
253
254

        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
255
        by the luminosity distance of the source.
Yannick Roehlly's avatar
Yannick Roehlly committed
256

Yannick Roehlly's avatar
Yannick Roehlly committed
257
        Fλ is in W/m²/nm. At redshift 0, the flux is computed at 10 pc. Then,
Yannick Roehlly's avatar
Yannick Roehlly committed
258
259
260
261
262
263
264
265
266
        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
267
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
268
        ----------
269
270
271
        filter_name: string
            Name of the filter to integrate into. It must be presnt in the
            database.
Yannick Roehlly's avatar
Yannick Roehlly committed
272
273
274

        Return
        ------
275
        fnu: float
Yannick Roehlly's avatar
Yannick Roehlly committed
276
277
278
            The integrated Fν density in mJy.
        """

279
        wavelength = self.wavelength_grid
Yannick Roehlly's avatar
Yannick Roehlly committed
280

281
282
283
284
285
286
287
        # First we try to fetch the filter's wavelength, transmission and
        # effective wavelength from the cache. The two keys are the size of the
        # spectrum wavelength grid and the name of the filter. The first key is
        # necessary because different spectra may have different sampling. To
        # avoid having the resample the filter every time on the optimal grid
        # (spectrum+filter), we store the resampled filter. That way we only
        # have to resample to spectrum.
Médéric Boquien's avatar
Médéric Boquien committed
288
289
290
        if 'universe.redshift' in self.info:
            key = (wavelength.size, filter_name,
                   self.info['universe.redshift'])
291
            dist = self.info['universe.luminosity_distance']
292
293
        else:
            key = (wavelength.size, filter_name, 0.)
294
295
            dist = 10. * parsec

296
        if key in self.cache_filters:
297
            wavelength_r, transmission_r, lambda_eff = self.cache_filters[key]
Yannick Roehlly's avatar
Yannick Roehlly committed
298
        else:
299
300
301
302
303
304
305
306
307
308
309
310
            with Database() as db:
                filter_ = db.get_filter(filter_name)
            trans_table = filter_.trans_table
            lambda_eff = filter_.effective_wavelength
            lambda_min = filter_.trans_table[0][0]
            lambda_max = filter_.trans_table[0][-1]

            # Test if the filter covers all the spectrum extent. If not then
            # the flux is not defined
            if ((wavelength[0] > lambda_min) or (wavelength[-1] < lambda_max)):
                return -99.

Yannick Roehlly's avatar
Yannick Roehlly committed
311
312
313
            # 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
314
            # filter one.
Médéric Boquien's avatar
Médéric Boquien committed
315
316
            w = np.where((wavelength >= lambda_min) &
                         (wavelength <= lambda_max))
317
            wavelength_r = utils.best_grid(wavelength[w], trans_table[0])
Médéric Boquien's avatar
Médéric Boquien committed
318
319
            transmission_r = np.interp(wavelength_r, trans_table[0],
                                       trans_table[1])
Yannick Roehlly's avatar
Yannick Roehlly committed
320

Médéric Boquien's avatar
Médéric Boquien committed
321
322
            self.cache_filters[key] = (wavelength_r, transmission_r,
                                       lambda_eff)
323

324
        l_lambda_r = np.interp(wavelength_r, wavelength, self.luminosity)
Yannick Roehlly's avatar
Yannick Roehlly committed
325

326
        f_lambda = utils.luminosity_to_flux(
327
            utils.flux_trapz(transmission_r * l_lambda_r, wavelength_r, key),
328
            dist)
Yannick Roehlly's avatar
Yannick Roehlly committed
329

330
331
        # Return Fν in mJy. The 1e-9 factor is because λ is in nm and 1e29 for
        # convert from W/m²/Hz to mJy.
332
        return 1e-9 / c * 1e29 * lambda_eff * lambda_eff * f_lambda
333
334
335
336
337
338
339

    def to_votable(self, filename, mass=1.):
        """
        Save the SED to a VO-table file

        Parameters
        ----------
340
        filename: string
341
            Name of the VO-table file
342
        mass: float
343
344
345
346
347
            Galaxy mass in solar mass. When need, the saved data will be
            multiplied by this mass.

        """
        save_sed_to_vo(self, filename, mass)
348
349

    def copy(self):
350
351
352
353
354
355
        """
        Create a new copy of the object. This is done manually rather than
        using copy.deepcopy() for speed reasons. As we know the structure of
        the object, we can do a better job.

        """
356
        sed = SED()
357
        if self._sfh is not None:
358
            sed._sfh = (self._sfh[0], self._sfh[1])
359
360
361
362
363
364
365
        sed.modules = self.modules[:]
        if self.wavelength_grid is not None:
            sed.wavelength_grid = self.wavelength_grid.copy()
            sed.luminosity = self.luminosity.copy()
            sed.luminosities = self.luminosities.copy()
        sed.contribution_names = self.contribution_names[:]
        sed.info = self.info.copy()
366
        sed.mass_proportional_info = self.mass_proportional_info.copy()
367
368

        return sed