__init__.py 14.4 KB
Newer Older
Yannick Roehlly's avatar
Yannick Roehlly committed
1
# -*- coding: utf-8 -*-
Yannick Roehlly's avatar
Yannick Roehlly committed
2
"""Copyright (C) 2012, 2013 Centre de données Astrophysiques de Marseille
Yannick Roehlly's avatar
Yannick Roehlly committed
3
4
5
6
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
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

26
27
- lines: a dictionnary containing the emission lines associated with the SED.
  A dictionnary 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 dictionnary containing various information about the SED.
Yannick Roehlly's avatar
Yannick Roehlly committed
33

Yannick Roehlly's avatar
Yannick Roehlly committed
34
35
- mass_proportional_info: the list of keys in the info dictionnary whose value
  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
39


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

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
55
56
57
58
59
60
61
    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
62
63
64
        self.modules = []
        self.wavelength_grid = None
        self.contribution_names = []
Yannick Roehlly's avatar
Yannick Roehlly committed
65
        self.luminosities = None
66
        self.lines = {}
Yannick Roehlly's avatar
Yannick Roehlly committed
67
        self.info = {}
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

Yannick Roehlly's avatar
Yannick Roehlly committed
83
84
    @property
    def wavelength_grid(self):
Yannick Roehlly's avatar
Yannick Roehlly committed
85
        """ Return a copy of the wavelength grid
Yannick Roehlly's avatar
Yannick Roehlly committed
86
87
88
89
90
91
92
93
94
95
96
        """
        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
97
98
    def luminosities(self):
        """ Return a copy of the luminosity contributions
Yannick Roehlly's avatar
Yannick Roehlly committed
99
        """
Yannick Roehlly's avatar
Yannick Roehlly committed
100
        if self._luminosities is None:
Yannick Roehlly's avatar
Yannick Roehlly committed
101
102
            return None
        else:
Yannick Roehlly's avatar
Yannick Roehlly committed
103
            return np.copy(self._luminosities)
Yannick Roehlly's avatar
Yannick Roehlly committed
104

Yannick Roehlly's avatar
Yannick Roehlly committed
105
106
107
    @luminosities.setter
    def luminosities(self, value):
        self._luminosities = value
Yannick Roehlly's avatar
Yannick Roehlly committed
108
109
110

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

Yannick Roehlly's avatar
Yannick Roehlly committed
113
114
115
        Return the total luminosity density vector, i.e. the sum of all the
        contributions in W/nm.
        """
Yannick Roehlly's avatar
Yannick Roehlly committed
116
117
118
119
        if self._luminosities is None:
            return None
        else:
            return self._luminosities.sum(0)
Yannick Roehlly's avatar
Yannick Roehlly committed
120

Yannick Roehlly's avatar
Yannick Roehlly committed
121
    def lambda_fnu(self, redshift=0, apply_redshift=False):
Yannick Roehlly's avatar
Yannick Roehlly committed
122
123
124
125
        """
        Return the (redshifted if asked) total Fν flux density vs wavelength
        spectrum of the SED.

Yannick Roehlly's avatar
Yannick Roehlly committed
126
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
127
128
        ----------
        redshift : float, default = 0
Yannick Roehlly's avatar
Yannick Roehlly committed
129
130
            If 0 (the default), the flux at 10 pc is computed.
        apply_redshift : boolean, default = None
131
132
133
            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
134
135
136
137
138
139
140

        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
141
142
        wavelength = self.wavelength_grid
        luminosity = self.luminosity
Yannick Roehlly's avatar
Yannick Roehlly committed
143

Yannick Roehlly's avatar
Yannick Roehlly committed
144
145
146
        if apply_redshift:
            wavelength, luminosity = utils.redshift_spectrum(
                wavelength, luminosity, redshift)
Yannick Roehlly's avatar
Yannick Roehlly committed
147

Yannick Roehlly's avatar
Yannick Roehlly committed
148
149
150
151
152
        # 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
153
154
155

        return wavelength, f_nu

156
    def add_info(self, key, value, mass_proportional=False, force=False):
Yannick Roehlly's avatar
Yannick Roehlly committed
157
158
159
        """
        Add a key / value to the information dictionary

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

Yannick Roehlly's avatar
Yannick Roehlly committed
164
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
165
166
167
168
169
        ----------
        key : any immutable
           The key used to retrieve the information.
        value : anything
           The information.
170
171
172
        mass_proportional : boolean
           If True, the added variable is set as proportional to the
           mass.
173
        force : boolean
174
           If false (default), adding a key that already exists in the info
175
176
           dictionary will raise an error. If true, doing this will update
           the associated value.
Yannick Roehlly's avatar
Yannick Roehlly committed
177
178

        """
179
        if (key not in self.info) or force:
Yannick Roehlly's avatar
Yannick Roehlly committed
180
            self.info[key] = value
181
182
            if mass_proportional:
                self.mass_proportional_info.append(key)
Yannick Roehlly's avatar
Yannick Roehlly committed
183
        else:
184
            raise KeyError("The information %s is already present "
Yannick Roehlly's avatar
Yannick Roehlly committed
185
186
                           "in the SED. " % key)

187
188
189
    def add_module(self, module_name, module_conf):
        """Add a new module information to the SED.

Yannick Roehlly's avatar
Yannick Roehlly committed
190
        Parameters
191
192
193
194
195
        ----------
        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
196
            Dictionary containing the module parameters.
197

Yannick Roehlly's avatar
Yannick Roehlly committed
198
        TODO: Complete the parameter dictionary with the default values from
199
200
              the module if they are not present.

Yannick Roehlly's avatar
Yannick Roehlly committed
201
        """
202
        self.modules.append((module_name, module_conf))
Yannick Roehlly's avatar
Yannick Roehlly committed
203

204
205
206
207
    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
208
209
210
211
212
213
214
215

        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
216
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
        ----------
        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
234
        if self.luminosities is None:
Yannick Roehlly's avatar
Yannick Roehlly committed
235
            self.wavelength_grid = np.copy(results_wavelengths)
Yannick Roehlly's avatar
Yannick Roehlly committed
236
            self.luminosities = np.copy(results_lumin)
Yannick Roehlly's avatar
Yannick Roehlly committed
237
238
239
240
241
242
243
        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
244
245
246
247
248
249
250
            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
251
252
253
254
            interp_lumin = interp1d(results_wavelengths,
                                    results_lumin,
                                    bounds_error=False,
                                    fill_value=0)(new_wavelength_grid)
Yannick Roehlly's avatar
Yannick Roehlly committed
255
256

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

259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
    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
285
286
287
288
289
290
    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
291
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
292
293
294
295
296
297
298
299
300
301
302
303
304
305
        ----------
        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
306
        return self.luminosities[idx]
Yannick Roehlly's avatar
Yannick Roehlly committed
307

308
    def compute_fnu(self, transmission, lambda_eff,
309
310
                    redshift=0, apply_redshift=False,
                    add_line_fluxes=True):
Yannick Roehlly's avatar
Yannick Roehlly committed
311
312
313
314
315
316
317
318
319
320
321
322
        """
        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
323
        Fλ is in W/m²/nm. At redshift 0, the flux is computed at 10 pc. Then,
Yannick Roehlly's avatar
Yannick Roehlly committed
324
325
326
327
328
329
330
331
332
        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
333
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
334
335
336
337
338
339
340
341
342
        ----------
        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
343
            The redshift of the galaxy. If 0, the flux is computed at 10 pc.
Yannick Roehlly's avatar
Yannick Roehlly committed
344

Yannick Roehlly's avatar
Yannick Roehlly committed
345
        apply_redshift : boolean
346
347
348
349
            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.

350
351
352
353
        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
354
355
356
357
358
359
        Return
        ------
        fnu : float
            The integrated Fν density in mJy.
        """

360
        # Filter limits
Yannick Roehlly's avatar
Yannick Roehlly committed
361
362
363
        lambda_min = min(transmission[0])
        lambda_max = max(transmission[0])

364
365
366
367
368
369
370
        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
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
        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
395
            # already normalised?
Yannick Roehlly's avatar
Yannick Roehlly committed
396
397
398
399
400
401
            f_lambda = utils.luminosity_to_flux(
                (np.trapz(transmission_r * l_lambda_r, wavelength_r) /
                 np.trapz(transmission_r, wavelength_r)),
                redshift
            )

402
            # Add the Fλ fluxes from the spectral lines.
403
            # TODO write the code
404

Yannick Roehlly's avatar
Yannick Roehlly committed
405
406
407
408
409
410
411
            # 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