Resample a spectrum to a new wavelengths map while conserving total flux.

    :param old_wavelengths:
        The original wavelengths array.

    :type old_wavelengths:
        :class:`numpy.array`

    :param new_wavelengths:
        The new wavelengths array to resample onto.

    :type new_wavelengths:
        :class:`numpy.array`

def resample(old_wavelengths, new_wavelengths):
    """
    Resample a spectrum to a new wavelengths map while conserving total flux.

    :param old_wavelengths:
        The original wavelengths array.

    :type old_wavelengths:
        :class:`numpy.array`

    :param new_wavelengths:
        The new wavelengths array to resample onto.

    :type new_wavelengths:
        :class:`numpy.array`
    """

    data = []
    old_px_indices = []
    new_px_indices = []
    for i, new_wl_i in enumerate(new_wavelengths):

        # These indices should span just over the new wavelength pixel.
        indices = np.unique(np.clip(
            old_wavelengths.searchsorted(new_wavelengths[i:i + 2], side="left")\
                + [-1, +1], 0, old_wavelengths.size - 1))
        N = np.ptp(indices)

        if N == 0:
            # 'Fake' pixel.
            data.append(np.nan)
            new_px_indices.append(i)
            old_px_indices.extend(indices)
            continue

        # Sanity checks.
        assert (old_wavelengths[indices[0]] <= new_wl_i \
            or indices[0] == 0)
        assert (new_wl_i <= old_wavelengths[indices[1]] \
            or indices[1] == old_wavelengths.size - 1)

        fractions = np.ones(N)

        # Edges are handled as fractions between rebinned pixels.
        _ = np.clip(i + 1, 0, new_wavelengths.size - 1)
        lhs = old_wavelengths[indices[0]:indices[0] + 2]
        rhs = old_wavelengths[indices[-1] - 1:indices[-1] + 1]
        fractions[0]  = (lhs[1] - new_wavelengths[i])/np.ptp(lhs)
        fractions[-1] = (new_wavelengths[_] - rhs[0])/np.ptp(rhs)

        # Being binned to a single pixel. Prevent overflow from fringe cases.
        fractions = np.clip(fractions, 0, 1)
        fractions /= fractions.sum()

        data.extend(fractions) 
        new_px_indices.extend([i] * N) # Mark the new pixel indices affected.
        old_px_indices.extend(np.arange(*indices)) # And the old pixel indices.

    return sparse.csc_matrix((data, (old_px_indices, new_px_indices)),
        shape=(old_wavelengths.size, new_wavelengths.size))