utils.py 5.85 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
# Authors: Yannick Roehlly, Médéric Boquien
Yannick Roehlly's avatar
Yannick Roehlly committed
5
6

import numpy as np
7
from scipy.constants import c, pi
Yannick Roehlly's avatar
Yannick Roehlly committed
8
9
10
11
12


def lambda_to_nu(wavelength):
    """Convert wavelength (nm) to frequency (Hz)

Yannick Roehlly's avatar
Yannick Roehlly committed
13
    Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
14
    ----------
15
    wavelength: float or array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
16
17
18
19
        The wavelength(s) in nm.

    Returns
    -------
20
    nu: float or array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
21
22
23
24
25
26
27
28
29
        The frequency(ies) in Hz.

    """
    return c / (wavelength * 1.e-9)


def nu_to_lambda(frequency):
    """Convert frequency (Hz) to wavelength (nm)

Yannick Roehlly's avatar
Yannick Roehlly committed
30
    Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
31
    ----------
32
    frequency: float or numpy.array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
33
34
35
36
        The frequency(ies) in Hz.

    Returns
    -------
37
    wavelength: float or numpy.array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
38
39
40
41
42
        The wavelength(s) in nm.

    """
    return 1.e-9 * c / frequency

Médéric Boquien's avatar
Médéric Boquien committed
43

44
45
46
47
48
49
50
51
52
53
54
def best_grid_memoise(f):
    """
    Memoisation helper for the best_grid() function. Given that best_grid takes
    numpy arrays that are not hashable as input, we cannot use the standard
    memoisation functions. Here we use the array sizes, minimum and maximum
    values as hashes.

    Parameters
    ----------
    f: function
        Function to memoise.
Yannick Roehlly's avatar
Yannick Roehlly committed
55

56
57
58
59
60
61
    Returns
    -------
    best_grid_helper: function
        Meomoised best_grid function
    """
    memo = {}
Médéric Boquien's avatar
Médéric Boquien committed
62
63

    def best_grid_helper(x, y):
64
65
66
67
68
69
70
71
72
73
74
        sx = x.size
        minx = x[0]
        maxx = x[-1]
        sy = y.size
        miny = y[0]
        maxy = y[-1]
        if (sx, minx, maxx, sy, miny, maxy) not in memo:
            memo[(sx, minx, maxx, sy, miny, maxy)] = f(x, y)
        return memo[(sx, minx, maxx, sy, miny, maxy)]
    return best_grid_helper

Médéric Boquien's avatar
Médéric Boquien committed
75

76
@best_grid_memoise
Yannick Roehlly's avatar
Yannick Roehlly committed
77
78
79
80
def best_grid(wavelengths1, wavelengths2):
    """
    Return the best wavelength grid to regrid to arrays

Yannick Roehlly's avatar
Yannick Roehlly committed
81
    Considering the two wavelength grids passed in parameters, this function
Yannick Roehlly's avatar
Yannick Roehlly committed
82
    compute the best new grid that will be used to regrid the two spectra
83
84
    before combining them. We do not use np.unique as it is much slowe than
    finding the unique elements by hand.
Yannick Roehlly's avatar
Yannick Roehlly committed
85

Yannick Roehlly's avatar
Yannick Roehlly committed
86
    Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
87
    ----------
88
    wavelengths1, wavelengths2: array of floats
89
        The wavelength grids to be 'regridded'.
Yannick Roehlly's avatar
Yannick Roehlly committed
90
91
92

    Returns
    -------
93
    new_grid: array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
94
95
96
        Array containing all the wavelengths found in the input arrays.

    """
97

98
    wl = np.concatenate((wavelengths1, wavelengths2))
99
    wl.sort(kind='mergesort')
100
101
102
    flag = np.ones(len(wl), dtype=bool)
    np.not_equal(wl[1:], wl[:-1], out=flag[1:])
    return wl[flag]
Yannick Roehlly's avatar
Yannick Roehlly committed
103

Médéric Boquien's avatar
Médéric Boquien committed
104

105
def luminosity_to_flux(luminosity, dist):
Yannick Roehlly's avatar
Yannick Roehlly committed
106
107
108
109
110
    """
    Convert a luminosity (or luminosity density) to a flux (or flux density).

    F = L / (4πDl2)

Yannick Roehlly's avatar
Yannick Roehlly committed
111
    Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
112
    ----------
113
    luminosity: float or array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
114
        Luminosity (typically in W) or luminosity density (W/nm or W/Hz).
115
116
    dist: float
        Luminosity distance of the object in metres
Yannick Roehlly's avatar
Yannick Roehlly committed
117
118
119

    Returns
    -------
120
    flux: float or array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
121
122
123
124
        The flux (typically in W/m²) of flux density (W/m²/nm or W/m²/Hz).

    """

125
    flux = luminosity / (4. * pi * dist * dist)
126

127
    return flux
Yannick Roehlly's avatar
Yannick Roehlly committed
128
129


Yannick Roehlly's avatar
Yannick Roehlly committed
130
def lambda_flambda_to_fnu(wavelength, flambda):
Yannick Roehlly's avatar
Yannick Roehlly committed
131
132
133
    """
    Convert a Fλ vs λ spectrum to Fν vs λ

Yannick Roehlly's avatar
Yannick Roehlly committed
134
    Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
135
    ----------
136
    wavelength: list-like of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
137
        The wavelengths in nm.
138
    flambda: list-like of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
139
        Fλ flux density in W/m²/nm (or Lλ luminosity density in W/nm).
Yannick Roehlly's avatar
Yannick Roehlly committed
140
141
142

    Returns
    -------
143
    fnu: array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
144
145
        The Fν flux density in mJy (or the Lν luminosity density in
        1.e-29 W/Hz).
Yannick Roehlly's avatar
Yannick Roehlly committed
146
147

    """
Yannick Roehlly's avatar
Yannick Roehlly committed
148
149
150
151
    wavelength = np.array(wavelength, dtype=float)
    flambda = np.array(flambda, dtype=float)

    # Factor 1e+29 is to switch from W/m²/Hz to mJy
Yannick Roehlly's avatar
Yannick Roehlly committed
152
    # Factor 1e-9 is to switch from nm to m (only one because the other nm
Yannick Roehlly's avatar
Yannick Roehlly committed
153
154
    # wavelength goes with the Fλ in W/m²/nm).
    fnu = 1e+29 * 1e-9 * flambda * wavelength * wavelength / c
Yannick Roehlly's avatar
Yannick Roehlly committed
155

Yannick Roehlly's avatar
Yannick Roehlly committed
156
    return fnu
Yannick Roehlly's avatar
Yannick Roehlly committed
157
158


Yannick Roehlly's avatar
Yannick Roehlly committed
159
def lambda_fnu_to_flambda(wavelength, fnu):
Yannick Roehlly's avatar
Yannick Roehlly committed
160
161
162
    """
    Convert a Fν vs λ spectrum to Fλ vs λ

Yannick Roehlly's avatar
Yannick Roehlly committed
163
    Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
164
    ----------
165
    wavelength: list-like of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
166
        The wavelengths in nm.
167
    fnu: list-like of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
168
169
        The Fν flux density in mJy (of the  Lν luminosity density in
        1.e-29 W/Hz).
Yannick Roehlly's avatar
Yannick Roehlly committed
170
171
172

    Returns
    -------
173
    flambda: array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
174
        Fλ flux density in W/m²/nm (or Lλ luminosity density in W/nm).
Yannick Roehlly's avatar
Yannick Roehlly committed
175
176

    """
Yannick Roehlly's avatar
Yannick Roehlly committed
177
178
    wavelength = np.array(wavelength, dtype=float)
    fnu = np.array(fnu, dtype=float)
Yannick Roehlly's avatar
Yannick Roehlly committed
179

Yannick Roehlly's avatar
Yannick Roehlly committed
180
181
182
    # Factor 1e-29 is to switch from Jy to W/m²/Hz
    # Factor 1e+9 is to switch from m to nm
    flambda = 1e-29 * 1e+9 * fnu / (wavelength * wavelength) * c
Yannick Roehlly's avatar
Yannick Roehlly committed
183

Yannick Roehlly's avatar
Yannick Roehlly committed
184
    return flambda
Yannick Roehlly's avatar
Yannick Roehlly committed
185
186


Yannick Roehlly's avatar
Yannick Roehlly committed
187
188
def redshift_spectrum(wavelength, flux, redshift, is_fnu=False):
    """Redshit a spectrum
Yannick Roehlly's avatar
Yannick Roehlly committed
189

Yannick Roehlly's avatar
Yannick Roehlly committed
190
    Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
191
    ----------
192
    wavelength: array like of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
193
        The wavelength in nm.
194
    flux: array like of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
195
        The flux or luminosity density.
196
    redshift: float
Yannick Roehlly's avatar
Yannick Roehlly committed
197
        The redshift.
198
    is_fnu: boolean
Yannick Roehlly's avatar
Yannick Roehlly committed
199
200
201
        If false (default) the flux is a Fλ density in W/m²/nm (or a Lλ
        luminosity density in W/nm). If true, the flux is a Fν density in mJy
        (or a Lν luminosity density in 1.e-29 W/Hz).
Yannick Roehlly's avatar
Yannick Roehlly committed
202
203
204

    Results
    -------
205
    wavelength, flux: tuple of numpy arrays of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
206
207
        The redshifted spectrum with the same kind of flux (or luminosity)
        density as the input.
Yannick Roehlly's avatar
Yannick Roehlly committed
208
209

    """
Yannick Roehlly's avatar
Yannick Roehlly committed
210
211
212
    wavelength = np.array(wavelength, dtype=float)
    flux = np.array(flux, dtype=float)
    redshift = float(redshift)
Yannick Roehlly's avatar
Yannick Roehlly committed
213

Yannick Roehlly's avatar
Yannick Roehlly committed
214
215
216
217
    if redshift < 0:
        redshift_factor = 1. / (1. - redshift)
    else:
        redshift_factor = 1. + redshift
Yannick Roehlly's avatar
Yannick Roehlly committed
218

Yannick Roehlly's avatar
Yannick Roehlly committed
219
220
221
    if is_fnu:
        # Switch to Fλ
        flux = lambda_fnu_to_flambda(wavelength, flux)
Yannick Roehlly's avatar
Yannick Roehlly committed
222

Yannick Roehlly's avatar
Yannick Roehlly committed
223
224
    wavelength *= redshift_factor
    flux /= redshift_factor
Yannick Roehlly's avatar
Yannick Roehlly committed
225

Yannick Roehlly's avatar
Yannick Roehlly committed
226
227
228
    if is_fnu:
        # Switch back to Fλ
        flux = lambda_flambda_to_fnu(wavelength, flux)
Yannick Roehlly's avatar
Yannick Roehlly committed
229

Yannick Roehlly's avatar
Yannick Roehlly committed
230
    return wavelength, flux