resample_spec.py

def resample_spectrum(source_wavelength, source_spectrum, destination_wavelengths, destination_resolution):
    """Resample a spectrum to a new wavelength grid, possibly degrading the spectrum's resoltuion.
    
    :param source_wavelength: The source wavelength array.
    :param source_spectrum: The source spectrum array.
    :param destination_wavelengths: The destination wavelength array.
    :param destination_resolution: The destination resolution.
    :returns: The resampled spectrum.
    
    This function is definitely valid for making a spectrum worse. Doing anything else is questionable.
    
    """
    # The main resampling function.
    # A two dimensional grid is used (instead of two for-loops). The grid stretches across wavelength (data), wavelength (requested). The 
    # standard deviation of the blurring gaussian corresponds point-for-point to the requested wavelengths. As such, we will have columns which
    # we can sum to make the new wavelengths, and rows which show the distribution of flux from each old wavelength into each new wavelength.
    sigma = destination_wavelengths / destination_resolution / 2.35
    curve = lambda wl,center,sig : (1.0/np.sqrt(np.pi * sig ** 2.0 )) * np.exp( - 0.5 * (wl - center) ** 2.0 / (sig ** 2.0) ).astype(np.float)

    # This is the two dimensional grid for the resampling function.
    MWL,MCENT = np.meshgrid(source_wavelength,destination_wavelengths)
    MWL,MSIGM = np.meshgrid(source_wavelength,sigma)

    # DO THE MATH!
    curves = curve(MWL,MCENT,MSIGM)

    # We then must normalize the light spread across each aperture by the gaussian. This makes sure the blurring gaussian only distributes
    # the amont of flux under each wavelength.
    ones = np.ones(source_wavelength.shape)
    base = np.sum(curves * ones, axis=1)
    top  = np.sum(curves * source_spectrum,axis=1)

    # If we try to normalize by dividing by zero, we are doing something wrong.
    # Removing these data points should be okay, because they are data points which we calculated to
    # have zero total flux anyways, and so we can ignore them.
    zeros = base == np.zeros(base.shape)
    topzo = top == np.zeros(base.shape)    
    # We don't actually clip those zero data points, we just make them into dumb numbers so that we don't get divide-by-zero errors
    base[zeros] = np.ones(np.sum(zeros))
    top[zeros] = np.zeros(np.sum(zeros))

    # Do the actual normalization
    flux = top  / base
    return flux