__init__.py 14.2 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

10
- sfh: the Star Formation History of the galaxy.
Yannick Roehlly's avatar
Yannick Roehlly committed
11

Yannick Roehlly's avatar
Yannick Roehlly committed
12
13
- 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
14

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

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

Yannick Roehlly's avatar
Yannick Roehlly committed
20
21
22
23
- 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
24

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

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

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

Yannick Roehlly's avatar
Yannick Roehlly committed
32
import numpy as np
33
from numpy.core.multiarray import interp # Compiled version
34
from scipy.constants import parsec
35

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

Yannick Roehlly's avatar
Yannick Roehlly committed
41

Yannick Roehlly's avatar
Yannick Roehlly committed
42
43
class SED(object):
    """Spectral Energy Distribution with associated information
Yannick Roehlly's avatar
Yannick Roehlly committed
44
45
    """

46
47
48
49
    # 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
50
51
52
53
54
    def __init__(self, sfh=None):
        """Create a new SED

        Parameters
        ----------
55
        sfh: (numpy.array, numpy.array)
Yannick Roehlly's avatar
Yannick Roehlly committed
56
57
58
59
60
61
            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
62
63
64
        self.modules = []
        self.wavelength_grid = None
        self.contribution_names = []
65
        self.luminosity = None
Yannick Roehlly's avatar
Yannick Roehlly committed
66
        self.luminosities = None
67
        self.lines = dict()
68
        self.info = dict()
69
        self.mass_proportional_info = set()
Yannick Roehlly's avatar
Yannick Roehlly committed
70

Yannick Roehlly's avatar
Yannick Roehlly committed
71
72
73
74
75
76
77
78
79
80
81
    @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
82
83
84
85

        # 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
86
87
        self._sfh = value

88
89
        if value is not None:
            sfh_sfr = value
90
            self._sfh = value
Yannick Roehlly's avatar
Yannick Roehlly committed
91
92
93
94
95
            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)
96
            self.add_info("sfh.age", sfh_sfr.size, False, force=True)
97

98
99
100
    @property
    def fnu(self):
        """Total Fν flux density of the SED
Yannick Roehlly's avatar
Yannick Roehlly committed
101

102
103
        Return the total Fν density vector, i.e the total luminosity converted
        to Fν flux in mJy.
Yannick Roehlly's avatar
Yannick Roehlly committed
104
105
        """

Yannick Roehlly's avatar
Yannick Roehlly committed
106
        # Fλ flux density in W/m²/nm
107
108
109
110
111
112
        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
113

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

117
        return f_nu
Yannick Roehlly's avatar
Yannick Roehlly committed
118

119
    def add_info(self, key, value, mass_proportional=False, force=False):
Yannick Roehlly's avatar
Yannick Roehlly committed
120
121
122
        """
        Add a key / value to the information dictionary

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

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

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

150
151
152
    def add_module(self, module_name, module_conf):
        """Add a new module information to the SED.

Yannick Roehlly's avatar
Yannick Roehlly committed
153
        Parameters
154
        ----------
155
        module_name: string
156
157
            Name of the module. This name can be suffixed with anything
            using a dot.
158
        module_conf: dictionary
Yannick Roehlly's avatar
Yannick Roehlly committed
159
            Dictionary containing the module parameters.
160

Yannick Roehlly's avatar
Yannick Roehlly committed
161
        TODO: Complete the parameter dictionary with the default values from
162
163
              the module if they are not present.

Yannick Roehlly's avatar
Yannick Roehlly committed
164
        """
165
        self.modules.append((module_name, module_conf))
Yannick Roehlly's avatar
Yannick Roehlly committed
166

167
168
169
170
    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
171
172
173
174
175
176
177
178

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

186
        results_wavelengths: array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
187
188
            The vector of the wavelengths of the module results (in nm).

189
        results_lumin: array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
190
191
192
193
194
195
196
            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.
197
        if self.luminosity is None:
198
199
200
            self.wavelength_grid = results_wavelengths.copy()
            self.luminosity = results_lumin.copy()
            self.luminosities = results_lumin.copy()
Yannick Roehlly's avatar
Yannick Roehlly committed
201
        else:
202
203
            # If the added luminosity contribution changes the SED wavelength
            # grid, we interpolate everything on a common wavelength grid.
204
205
206
207
            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.
208
209
                self.wavelength_grid, self.luminosities = \
                    utils.interpolate_lumin(self.wavelength_grid,
210
                                            self.luminosities,
211
212
213
                                            results_wavelengths,
                                            results_lumin)

214
                self.luminosity = self.luminosities.sum(0)
215
216
217
            else:
                self.luminosities = np.vstack((self.luminosities,
                                               results_lumin))
218
                self.luminosity += results_lumin
Yannick Roehlly's avatar
Yannick Roehlly committed
219
220
221
222
223
224
225

    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
226
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
227
        ----------
228
        name: string
Yannick Roehlly's avatar
Yannick Roehlly committed
229
230
231
232
            Name of the contribution

        Returns
        -------
233
        luminosities: array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
234
235
236
237
238
            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
239
240
        idx = (len(self.contribution_names) - 1 -
               self.contribution_names[::-1].index(name))
Yannick Roehlly's avatar
Yannick Roehlly committed
241
        return self.luminosities[idx]
Yannick Roehlly's avatar
Yannick Roehlly committed
242

243
    def compute_fnu(self, filter_name):
Yannick Roehlly's avatar
Yannick Roehlly committed
244
        """
245
        Compute the Fν flux density in a given filter
Yannick Roehlly's avatar
Yannick Roehlly committed
246

247
248
249
250
        The filters are stored in the database in such a way that after
        integration and conversion from luminosity to flux we directly get the
        latter in units of mJy. If the SED spectrum does not cover all the
        filter response table, NaN is returned.
Yannick Roehlly's avatar
Yannick Roehlly committed
251

Yannick Roehlly's avatar
Yannick Roehlly committed
252
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
253
        ----------
254
255
256
        filter_name: string
            Name of the filter to integrate into. It must be presnt in the
            database.
Yannick Roehlly's avatar
Yannick Roehlly committed
257
258
259

        Return
        ------
260
        fnu: float
Yannick Roehlly's avatar
Yannick Roehlly committed
261
262
263
            The integrated Fν density in mJy.
        """

264
        wavelength = self.wavelength_grid
Yannick Roehlly's avatar
Yannick Roehlly committed
265

266
        # First we try to fetch the filter's wavelength, transmission and
267
        # pivot wavelength from the cache. The two keys are the size of the
268
269
270
271
272
        # 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
273
        if 'universe.redshift' in self.info:
274
275
276
277
278
279
280
            if 'nebular.lines_width' in self.info:
                key = (wavelength.size, filter_name,
                       self.info['nebular.lines_width'],
                       self.info['universe.redshift'])
            else:
                key = (wavelength.size, filter_name,
                       self.info['universe.redshift'])
281
            dist = self.info['universe.luminosity_distance']
282
        else:
283
284
285
286
287
            if 'nebular.lines_width' in self.info:
                key = (wavelength.size, filter_name,
                       self.info['nebular.lines_width'], 0.)
            else:
                key = (wavelength.size, filter_name, 0.)
288
289
            dist = 10. * parsec

290
291
292
293
294
295
296
        if filter_name.startswith('line.'):
            lum = 0
            for name in filter_name.split('+'):
                line = self.lines[name.split('.', maxsplit=1)[1]]
                lum += line[1] + line[2]  # Young and old components
            return utils.luminosity_to_flux(lum, dist)

297
        if key in self.cache_filters:
298
            wavelength_r, transmission_r, lambda_piv = self.cache_filters[key]
Yannick Roehlly's avatar
Yannick Roehlly committed
299
        else:
300
301
302
            with Database() as db:
                filter_ = db.get_filter(filter_name)
            trans_table = filter_.trans_table
303
            lambda_piv = filter_.pivot_wavelength
304
305
306
307
308
309
            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)):
310
                return np.nan
311

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

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

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

327
328
329
330
331
332
        # We compute directly Fν from ∫T×Fλ×dλ/∫T×c/λ²×dλ. The filter bandpass
        # in the database is already normalised so that we do not need to
        # compute the denominator (it is a constant that does not depend on the
        # spectrum) and the normalisation is such that the results we obtain
        # are directly in mJy.
        f_nu = utils.luminosity_to_flux(
333
            utils.flux_trapz(transmission_r * l_lambda_r, wavelength_r, key),
334
            dist)
Yannick Roehlly's avatar
Yannick Roehlly committed
335

336
        return f_nu
337
338
339
340
341
342
343

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

        Parameters
        ----------
344
        filename: string
345
            Name of the VO-table file
346
        mass: float
347
348
349
350
351
            Galaxy mass in solar mass. When need, the saved data will be
            multiplied by this mass.

        """
        save_sed_to_vo(self, filename, mass)
352

353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
    def to_fits(self, prefix, mass=1.):
        """
        Save the SED to FITS files.

        Parameters
        ----------
        prefix: string
            Prefix of the fits file containing the path and the id of the model
        mass: float
            Galaxy mass in solar masses. When needed, the data will be scaled
            to this mass

        """
        save_sed_to_fits(self, prefix, mass)

368
    def copy(self):
369
370
371
372
373
374
        """
        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.

        """
375
        sed = SED()
376
        if self._sfh is not None:
377
            sed._sfh = self._sfh
378
379
380
381
382
383
        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[:]
384
        sed.lines = self.lines.copy()
385
        sed.info = self.info.copy()
386
        sed.mass_proportional_info = self.mass_proportional_info.copy()
387
388

        return sed