__init__.py 12.4 KB
Newer Older
Yannick Roehlly's avatar
Yannick Roehlly committed
1
2
# -*- coding: utf-8 -*-
"""
3
Copyright (C) 2012, 2013 Centre de données Astrophysiques de Marseille
Yannick Roehlly's avatar
Yannick Roehlly committed
4
5
6
7
8
9
10
11
12
13
14
15
16
Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt

@author: Yannick Roehlly <yannick.roehlly@oamp.fr>

"""


import numpy as np
from . import utils
from scipy.constants import c


class SED(object):
17
    """Spectral Energy Distribution with associated information
Yannick Roehlly's avatar
Yannick Roehlly committed
18
19
20
21

    This class represents a Spectral Energy Distribution (SED) as constructed
    by pCigale. Such a SED is characterised by:

Yannick Roehlly's avatar
Yannick Roehlly committed
22
    - A list of tuples (module name, parameter dictionary) describing all the
Yannick Roehlly's avatar
Yannick Roehlly committed
23
24
25
26
27
28
29
30
31
32
33
34
35
      pCigale modules the SED 'went through'.

    - The wavelengths grid (in nm).

    - An array of luminosity densities (in W/nm) containing the luminosity
      contribution (in positive or negative) of each module the SED 'went
      through'. The first axis corresponds to the index in the contribution
      list (see below) and the second axis corresponds to the wavelength grid.

    - The list of the contributions that are in the above array. This list is
      separated from the list of the modules so that one module can result in
      various contributions to the SED.

Yannick Roehlly's avatar
Yannick Roehlly committed
36
    - A dictionary of arbitrary information associated with the SED.
Yannick Roehlly's avatar
Yannick Roehlly committed
37

38
39
40
    - The list of the keys from the info dictionary whose value is
      proportional to the galaxy mass.

Yannick Roehlly's avatar
Yannick Roehlly committed
41
42
43
44
45
46
47
48
    """

    def __init__(self):
        self.modules = []
        self.wavelength_grid = None
        self.lumin_contributions = None
        self.contribution_names = []
        self.info = {}
49
        self.mass_proportional_info = []
Yannick Roehlly's avatar
Yannick Roehlly committed
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86

    @property
    def wavelength_grid(self):
        """ Return a copy of the wavelength grid to avoid involuntary
        modifications.
        """
        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
    def lumin_contributions(self):
        """ Return a copy of the luminosity contribution table to avoid
        involuntary modifications.
        """
        if self._lumin_contributions is None:
            return None
        else:
            return np.copy(self._lumin_contributions)

    @lumin_contributions.setter
    def lumin_contributions(self, value):
        self._lumin_contributions = value

    @property
    def luminosity(self):
        """
        Return the total luminosity density vector, i.e. the sum of all the
        contributions in W/nm.
        """
        return self.lumin_contributions.sum(0)

87
    def lambda_fnu(self, redshift=0, redshift_spectrum=False):
Yannick Roehlly's avatar
Yannick Roehlly committed
88
89
90
91
        """
        Return the (redshifted if asked) total Fν flux density vs wavelength
        spectrum of the SED.

Yannick Roehlly's avatar
Yannick Roehlly committed
92
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
93
94
95
        ----------
        redshift : float, default = 0
            If 0 (the default), the flux at 10 pc is computed.
96
97
98
99
        redshift_spectrum : boolean, default = None
            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
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114

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

        """
        # Fλ flux density in W/m²/nm
        f_lambda = utils.luminosity_to_flux(self.luminosity, redshift)

        # The 1.e29 factor is to convert from W/m²/Hz to mJy
        f_nu = (f_lambda * self.wavelength_grid
                * 1.e-9 * self.wavelength_grid / c
                * 1.e29)

115
116
117
118
119
        if redshift_spectrum:
            wavelength = utils.redshift_wavelength(self.wavelength_grid,
                                                   redshift)
        else:
            wavelength = np.copy(self.wavelength_grid)
Yannick Roehlly's avatar
Yannick Roehlly committed
120
121
122

        return wavelength, f_nu

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

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

Yannick Roehlly's avatar
Yannick Roehlly committed
131
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
132
133
134
135
136
        ----------
        key : any immutable
           The key used to retrieve the information.
        value : anything
           The information.
137
138
139
        mass_proportional : boolean
           If True, the added variable is set as proportional to the
           mass.
Yannick Roehlly's avatar
Yannick Roehlly committed
140
141
142
143

        """
        if key not in self.info:
            self.info[key] = value
144
145
            if mass_proportional:
                self.mass_proportional_info.append(key)
Yannick Roehlly's avatar
Yannick Roehlly committed
146
147
148
149
        else:
            raise KeyError("The information %s is yet present "
                           "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
156
157
158
        ----------
        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
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
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
        ----------
        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.
        if self.lumin_contributions is None:
            self.wavelength_grid = np.copy(results_wavelengths)
            self.lumin_contributions = np.array([np.copy(results_lumin)])
        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.
            new_lumin_table = None
            for old_lumin in self.lumin_contributions:
                interp_lumin = np.interp(
                    new_wavelength_grid,
                    self.wavelength_grid,
                    old_lumin,
                    right=0,
                    left=0)

                if new_lumin_table is None:
                    new_lumin_table = np.array([interp_lumin])
                else:
                    new_lumin_table = np.vstack((new_lumin_table,
                                                 interp_lumin))

            # Add the new module luminosities
            interp_lumin = np.interp(
                new_wavelength_grid,
                results_wavelengths,
                results_lumin,
                right=0,
                left=0)
            new_lumin_table = np.vstack((new_lumin_table, interp_lumin))

            self.wavelength_grid = new_wavelength_grid
            self.lumin_contributions = new_lumin_table

    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
240
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
        ----------
        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))
        return self.lumin_contributions[idx]

257
258
    def compute_fnu(self, transmission, lambda_eff,
                    redshift=0, redshift_spectrum=False):
Yannick Roehlly's avatar
Yannick Roehlly committed
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
        """
        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λ ) )

        Fλ is in W/m²/nm. At redshift 0, the flux is computed at 10 pc. Then,
        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
281
        Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
282
283
284
285
286
287
288
289
290
291
292
        ----------
        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
            The redshift of the galaxy. If 0, the flux is computed at 10 pc.

293
294
295
296
297
        redshift_spectrum : boolean
            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
298
299
300
301
302
303
304
305
306
307
308
309
310
311
        Return
        ------
        fnu : float
            The integrated Fν density in mJy.
        """

        lambda_min = min(transmission[0])
        lambda_max = max(transmission[0])

        if ((min(self.wavelength_grid) > lambda_min) or
                (max(self.wavelength_grid) < lambda_max)):
            f_nu = -99.

        else:
312
313
314
315
316
317
            if redshift_spectrum:
                wavelength = utils.redshift_wavelength(self.wavelength_grid,
                                                       redshift)
            else:
                wavelength = np.copy(self.wavelength_grid)

Yannick Roehlly's avatar
Yannick Roehlly committed
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
            l_lambda = self.luminosity

            # 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
            # yet normalised?
            f_lambda = utils.luminosity_to_flux(
                (np.trapz(transmission_r * l_lambda_r, wavelength_r) /
                 np.trapz(transmission_r, wavelength_r)),
                redshift
            )

            # 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