licks.py 文件源码

def reduce_resolution(wi, fi, sigma0=0.55, sigma_floor=0.2):
    """ Adapt the resolution of the spectra to match the lick definitions

        Lick definitions have different resolution elements as function
        of wavelength. These definition are hard-coded in this function

    Parameters
    ---------
    wi: ndarray (n, )
        wavelength definition
    fi: ndarray (nspec, n) or (n, )
        spectra to convert
    sigma0: float
        initial broadening in the spectra `fi`
    sigma_floor: float
        minimal dispersion to consider

    Returns
    -------
    flux_red: ndarray (nspec, n) or (n, )
        reduced spectra
    """

    # all in AA
    w_lick_res = (4000., 4400., 4900., 5400., 6000.)
    lick_res = (11.5, 9.2, 8.4, 8.4, 9.8)   # FWHM in AA

    w = np.asarray(wi)
    flux = np.atleast_2d(fi)

    # Linear interpolation of lick_res over w
    # numpy interp does constant instead of extrapolation
    # res = np.interp(w, w_lick_res, lick_res)

    # spline order: 1 linear, 2 quadratic, 3 cubic ...
    from scipy.interpolate import InterpolatedUnivariateSpline
    res = InterpolatedUnivariateSpline(w_lick_res, lick_res, k=1)(w)

    # Compute width from fwhm
    const = 2. * np.sqrt(2. * np.log(2))  # conversion fwhm --> sigma
    lick_sigma = np.sqrt((res ** 2 - sigma0 ** 2)) / const

    # Convolution by g=1/sqrt(2*pi*sigma^2) * exp(-r^2/(2*sigma^2))
    flux_red = np.zeros(flux.shape, dtype=flux.dtype)

    for i, sigma in enumerate(lick_sigma):
        maxsigma = 3. * sigma
        # sampling floor: min (0.2, sigma * 0.1)
        delta = min(sigma_floor, sigma * 0.1)
        delta_wj = np.arange(-maxsigma, + maxsigma, delta)
        wj = delta_wj + w[i]
        for k, fk in enumerate(flux):
            fluxj = np.interp(wj, w, fk, left=0., right=0.)
            flux_red[k, i] = np.sum(fluxj * delta * np.exp(-0.5 * (delta_wj / sigma) ** 2))

    flux_red /= lick_sigma * const

    return flux_red.reshape(np.shape(fi))