Created
May 9, 2016 22:24
-
-
Save alexrudy/8d3ae55cc27f0dccf1ab671e4e4bd7bc to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment