utils.py 5.84 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

43
44
45
46
47
48
49
50
51
52
53
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
54

55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
    Returns
    -------
    best_grid_helper: function
        Meomoised best_grid function
    """

    memo = {}
    def best_grid_helper(x,y):
        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

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

Yannick Roehlly's avatar
Yannick Roehlly committed
79
    Considering the two wavelength grids passed in parameters, this function
Yannick Roehlly's avatar
Yannick Roehlly committed
80
    compute the best new grid that will be used to regrid the two spectra
81
82
    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
83

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

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

    """
95

96
    wl = np.concatenate((wavelengths1, wavelengths2))
97
    wl.sort(kind='mergesort')
98
99
100
    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
101

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

    F = L / (4πDl2)

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

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

    """

122
    flux = luminosity / (4. * pi * dist * dist)
123

124
    return flux
Yannick Roehlly's avatar
Yannick Roehlly committed
125
126


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

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

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

    """
Yannick Roehlly's avatar
Yannick Roehlly committed
145
146
147
148
    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
149
    # Factor 1e-9 is to switch from nm to m (only one because the other nm
Yannick Roehlly's avatar
Yannick Roehlly committed
150
151
    # 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
152

Yannick Roehlly's avatar
Yannick Roehlly committed
153
    return fnu
Yannick Roehlly's avatar
Yannick Roehlly committed
154
155


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

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

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

    """
Yannick Roehlly's avatar
Yannick Roehlly committed
174
175
    wavelength = np.array(wavelength, dtype=float)
    fnu = np.array(fnu, dtype=float)
Yannick Roehlly's avatar
Yannick Roehlly committed
176

Yannick Roehlly's avatar
Yannick Roehlly committed
177
178
179
    # 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
180

Yannick Roehlly's avatar
Yannick Roehlly committed
181
    return flambda
Yannick Roehlly's avatar
Yannick Roehlly committed
182
183


Yannick Roehlly's avatar
Yannick Roehlly committed
184
185
def redshift_spectrum(wavelength, flux, redshift, is_fnu=False):
    """Redshit a spectrum
Yannick Roehlly's avatar
Yannick Roehlly committed
186

Yannick Roehlly's avatar
Yannick Roehlly committed
187
    Parameters
Yannick Roehlly's avatar
Yannick Roehlly committed
188
    ----------
189
    wavelength: array like of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
190
        The wavelength in nm.
191
    flux: array like of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
192
        The flux or luminosity density.
193
    redshift: float
Yannick Roehlly's avatar
Yannick Roehlly committed
194
        The redshift.
195
    is_fnu: boolean
Yannick Roehlly's avatar
Yannick Roehlly committed
196
197
198
        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
199
200
201

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

    """
Yannick Roehlly's avatar
Yannick Roehlly committed
207
208
209
    wavelength = np.array(wavelength, dtype=float)
    flux = np.array(flux, dtype=float)
    redshift = float(redshift)
Yannick Roehlly's avatar
Yannick Roehlly committed
210

Yannick Roehlly's avatar
Yannick Roehlly committed
211
212
213
214
    if redshift < 0:
        redshift_factor = 1. / (1. - redshift)
    else:
        redshift_factor = 1. + redshift
Yannick Roehlly's avatar
Yannick Roehlly committed
215

Yannick Roehlly's avatar
Yannick Roehlly committed
216
217
218
    if is_fnu:
        # Switch to Fλ
        flux = lambda_fnu_to_flambda(wavelength, flux)
Yannick Roehlly's avatar
Yannick Roehlly committed
219

Yannick Roehlly's avatar
Yannick Roehlly committed
220
221
    wavelength *= redshift_factor
    flux /= redshift_factor
Yannick Roehlly's avatar
Yannick Roehlly committed
222

Yannick Roehlly's avatar
Yannick Roehlly committed
223
224
225
    if is_fnu:
        # Switch back to Fλ
        flux = lambda_flambda_to_fnu(wavelength, flux)
Yannick Roehlly's avatar
Yannick Roehlly committed
226

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