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

Yannick Roehlly's avatar
Yannick Roehlly committed
28
- mass_proportional_info: the list 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
from collections import OrderedDict
Yannick Roehlly's avatar
Yannick Roehlly committed
35
from . import utils
36
from .io.vo import save_sed_to_vo
37
from scipy.constants import c, parsec
Yannick Roehlly's avatar
Yannick Roehlly committed
38
from scipy.interpolate import interp1d
Yannick Roehlly's avatar
Yannick Roehlly committed
39

40
41
42
43
44
# 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

45

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

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.info = OrderedDict()
68
        self.mass_proportional_info = []
Yannick Roehlly's avatar
Yannick Roehlly committed
69

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

83
84
85
86
87
88
89
        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)
90
            self.add_info("age", np.max(sfh_time), False, True)
91

Yannick Roehlly's avatar
Yannick Roehlly committed
92
93
    @property
    def wavelength_grid(self):
Yannick Roehlly's avatar
Yannick Roehlly committed
94
        """ Return a copy of the wavelength grid
Yannick Roehlly's avatar
Yannick Roehlly committed
95
96
97
98
99
100
101
102
103
104
105
        """
        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
106
107
    def luminosities(self):
        """ Return a copy of the luminosity contributions
Yannick Roehlly's avatar
Yannick Roehlly committed
108
        """
Yannick Roehlly's avatar
Yannick Roehlly committed
109
        if self._luminosities is None:
Yannick Roehlly's avatar
Yannick Roehlly committed
110
111
            return None
        else:
Yannick Roehlly's avatar
Yannick Roehlly committed
112
            return np.copy(self._luminosities)
Yannick Roehlly's avatar
Yannick Roehlly committed
113

Yannick Roehlly's avatar
Yannick Roehlly committed
114
115
116
    @luminosities.setter
    def luminosities(self, value):
        self._luminosities = value
Yannick Roehlly's avatar
Yannick Roehlly committed
117
118


119
120
121
    @property
    def fnu(self):
        """Total Fν flux density of the SED
Yannick Roehlly's avatar
Yannick Roehlly committed
122

123
124
        Return the total Fν density vector, i.e the total luminosity converted
        to Fν flux in mJy.
Yannick Roehlly's avatar
Yannick Roehlly committed
125
126
        """

Yannick Roehlly's avatar
Yannick Roehlly committed
127
        # Fλ flux density in W/m²/nm
128
129
        f_lambda = utils.luminosity_to_flux(self.luminosity,
                                            self.info['universe.luminosity_distance'])
Yannick Roehlly's avatar
Yannick Roehlly committed
130

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

134
        return f_nu
Yannick Roehlly's avatar
Yannick Roehlly committed
135

136
    def add_info(self, key, value, mass_proportional=False, force=False):
Yannick Roehlly's avatar
Yannick Roehlly committed
137
138
139
        """
        Add a key / value to the information dictionary

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

Yannick Roehlly's avatar
Yannick Roehlly committed
144
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
145
        ----------
146
        key: any immutable
Yannick Roehlly's avatar
Yannick Roehlly committed
147
           The key used to retrieve the information.
148
        value: anything
Yannick Roehlly's avatar
Yannick Roehlly committed
149
           The information.
150
        mass_proportional: boolean
151
152
           If True, the added variable is set as proportional to the
           mass.
153
        force: boolean
154
           If false (default), adding a key that already exists in the info
155
156
           dictionary will raise an error. If true, doing this will update
           the associated value.
Yannick Roehlly's avatar
Yannick Roehlly committed
157
158

        """
159
        if (key not in self.info) or force:
Yannick Roehlly's avatar
Yannick Roehlly committed
160
            self.info[key] = value
161
162
            if mass_proportional:
                self.mass_proportional_info.append(key)
Yannick Roehlly's avatar
Yannick Roehlly committed
163
        else:
164
            raise KeyError("The information %s is already present "
Yannick Roehlly's avatar
Yannick Roehlly committed
165
166
                           "in the SED. " % key)

167
168
169
    def add_module(self, module_name, module_conf):
        """Add a new module information to the SED.

Yannick Roehlly's avatar
Yannick Roehlly committed
170
        Parameters
171
        ----------
172
        module_name: string
173
174
            Name of the module. This name can be suffixed with anything
            using a dot.
175
        module_conf: dictionary
Yannick Roehlly's avatar
Yannick Roehlly committed
176
            Dictionary containing the module parameters.
177

Yannick Roehlly's avatar
Yannick Roehlly committed
178
        TODO: Complete the parameter dictionary with the default values from
179
180
              the module if they are not present.

Yannick Roehlly's avatar
Yannick Roehlly committed
181
        """
182
        self.modules.append((module_name, module_conf))
Yannick Roehlly's avatar
Yannick Roehlly committed
183

184
185
186
187
    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
188
189
190
191
192
193
194
195

        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
196
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
197
        ----------
198
        contribution_name: string
Yannick Roehlly's avatar
Yannick Roehlly committed
199
200
201
202
            Name of the contribution added. This name is used to retrieve the
            luminosity contribution and allows one module to add more than
            one contribution.

203
        results_wavelengths: array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
204
205
            The vector of the wavelengths of the module results (in nm).

206
        results_lumin: array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
207
208
209
210
211
212
213
            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.
214
        if self.luminosity is None:
Yannick Roehlly's avatar
Yannick Roehlly committed
215
            self.wavelength_grid = np.copy(results_wavelengths)
216
            self.luminosity = np.copy(results_lumin)
Yannick Roehlly's avatar
Yannick Roehlly committed
217
            self.luminosities = np.copy(results_lumin)
Yannick Roehlly's avatar
Yannick Roehlly committed
218
        else:
219
220
            # If the added luminosity contribution changes the SED wavelength
            # grid, we interpolate everything on a common wavelength grid.
221
222
223
224
225
226
227
228
229
230
231
            if (results_wavelengths.size != self.wavelength_grid.size or
                    not np.all(results_wavelengths == self.wavelength_grid)):
                # 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.
                new_luminosities = interp1d(self.wavelength_grid,
                                            self.luminosities,
                                            bounds_error=False,
232
                                            assume_sorted=True,
233
234
235
236
                                            fill_value=0.)(new_wavelength_grid)

                # Interpolate the added luminosity array to the new wavelength
                # grid
237
238
239
                interp_lumin = np.interp(new_wavelength_grid,
                                         results_wavelengths, results_lumin,
                                         left=0., right=0.)
Yannick Roehlly's avatar
Yannick Roehlly committed
240

241
242
                self.wavelength_grid = new_wavelength_grid
                self.luminosities = np.vstack((new_luminosities, interp_lumin))
243
                self.luminosity = self.luminosities.sum(0)
244
245
246
            else:
                self.luminosities = np.vstack((self.luminosities,
                                               results_lumin))
247
                self.luminosity += results_lumin
Yannick Roehlly's avatar
Yannick Roehlly committed
248
249
250
251
252
253
254

    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
255
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
256
        ----------
257
        name: string
Yannick Roehlly's avatar
Yannick Roehlly committed
258
259
260
261
            Name of the contribution

        Returns
        -------
262
        luminosities: array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
263
264
265
266
267
268
269
            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
270
        return self.luminosities[idx]
Yannick Roehlly's avatar
Yannick Roehlly committed
271

272
    def compute_fnu(self, transmission, lambda_eff):
Yannick Roehlly's avatar
Yannick Roehlly committed
273
274
275
276
277
278
279
280
281
282
283
284
        """
        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
285
        Fλ is in W/m²/nm. At redshift 0, the flux is computed at 10 pc. Then,
Yannick Roehlly's avatar
Yannick Roehlly committed
286
287
288
289
290
291
292
293
294
        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
295
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
296
        ----------
297
        transmission: 2D array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
298
299
300
            A numpy 2D array containing the filter response profile
            wavelength[nm] vs transmission).

301
        lambda_eff: float
Yannick Roehlly's avatar
Yannick Roehlly committed
302
303
304
305
            Effective wavelength of the filter in nm.

        Return
        ------
306
        fnu: float
Yannick Roehlly's avatar
Yannick Roehlly committed
307
308
309
            The integrated Fν density in mJy.
        """

310
        # Filter limits
311
312
        lambda_min = transmission[0][0]
        lambda_max = transmission[0][-1]
Yannick Roehlly's avatar
Yannick Roehlly committed
313

314
315
316
317
        wavelength = self.wavelength_grid
        l_lambda = self.luminosity

        # Test if the spectrum cover all the filter extend
318
319
        if ((wavelength[0] > lambda_min) or
                (wavelength[-1] < lambda_max)):
Yannick Roehlly's avatar
Yannick Roehlly committed
320
321
322
323
324
325
            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
326
327
328
329
            # filter one.

            w = np.where((wavelength >= lambda_min)&(wavelength <= lambda_max))
            wavelength_r = utils.best_grid(wavelength[w], transmission[0])
Yannick Roehlly's avatar
Yannick Roehlly committed
330
331
332
333
334

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

335
336
            if 'universe.luminosity_distance' in self.info.keys():
                dist = self.info['universe.luminosity_distance']
337
            else:
338
                dist = 10. * parsec
339

Yannick Roehlly's avatar
Yannick Roehlly committed
340
            f_lambda = utils.luminosity_to_flux(
341
342
                np.trapz(transmission_r * l_lambda_r, wavelength_r),
                dist)
Yannick Roehlly's avatar
Yannick Roehlly committed
343

344
345
            # Fν in W/m²/Hz. The 1e-9 factor is because λ is in nm.
            f_nu = lambda_eff * f_lambda * lambda_eff * 1e-9 / c
Yannick Roehlly's avatar
Yannick Roehlly committed
346
347

            # Conversion from W/m²/Hz to mJy
348
            f_nu *= 1e+29
Yannick Roehlly's avatar
Yannick Roehlly committed
349
350

        return f_nu
351
352
353
354
355
356
357

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

        Parameters
        ----------
358
        filename: string
359
            Name of the VO-table file
360
        mass: float
361
362
363
364
365
            Galaxy mass in solar mass. When need, the saved data will be
            multiplied by this mass.

        """
        save_sed_to_vo(self, filename, mass)