firefly instrument

This class handles instrumental properties that influence the FIREFLY fitting process.

Provides a set of functions to handle instrumental effects.

log_rebin() has been pulled from mangadap.contrib.ppxf_util.py and modified.

Source location:

$MANGADAP_DIR/python/mangadap/util/instrument.py

Python2/3 compliance:

from __future__ import division
from __future__ import print_function
from __future__ import absolute_import

import sys
if sys.version > '3':
    long = int

Imports:

import warnings
import numpy
from scipy.interpolate import InterpolatedUnivariateSpline
import astropy.constants

Revision history:

**24 May 2017 - Implemented the new downgrader version (D. Goddard)

class firefly_instrument.VariableGaussianKernel(sigma, minsig=0.01, nsig=3.0, integral=False)

Support class for variable sigma convolution. See convolution_variable_sigma().

Stolen from M. Cappellari’s gaussian_filter1d function.

Args:

y (np.ndarray): Vector to convolve sigma (np.ndarray): Coordinate-dependent standard deviation of the

Gaussian kernel

ye (np.ndarray): (Optional) Error in the vector to

convolve

Raises:

ValueError: Raised if y is not a 1D vector, or if the shape of

y and sigma (and ye if provided) are different.

Attributes:

x (np.ndarray): Pixel coordinate vector y (np.ndarray): Vector to convolve ye (np.ndarray): Error in the vector to convolve sigma (np.ndarray): Coordinate-dependent standard deviation of the

Gaussian kernel

norm (np.ndarray): Gaussian normalization; calculated once for

efficiency

class firefly_instrument.convolution_integral_element(y, sigma, ye=None)

Support class for variable sigma convolution. See convolution_variable_sigma().

OUT OF DATE; DO NOT USE

Args:

y (np.ndarray): Vector to convolve sigma (np.ndarray): Coordinate-dependent standard deviation of the

Gaussian kernel

ye (np.ndarray): (Optional) Error in the vector to

convolve

Raises:

ValueError: Raised if y is not a 1D vector, or if the shape of

y and sigma (and ye if provided) are different.

Attributes:

x (np.ndarray): Pixel coordinate vector y (np.ndarray): Vector to convolve ye (np.ndarray): Error in the vector to convolve sigma (np.ndarray): Coordinate-dependent standard deviation of the

Gaussian kernel

norm (np.ndarray): Gaussian normalization; calculated once for

efficiency

error(xc)

Calculates the error in the weighted mean of y using nominal error propagation. The weights are defined by a Gaussian with standard deviation sigma and centered at xc.

Args:

xc (float): Center for the Gaussian weighting function

Returns:

float: The error in the weighted mean of y

firefly_instrument.convolution_variable_sigma(y, sigma, ye=None, integral=False)

Convolve a discretely sampled function \(y(x)\) with a Gaussian kernel, \(g\), where the standard deviation of the kernel is a function of \(x\), \(\sigma(x)\). Nominal calculations can be performed to propagate the error in the result; these calculations do not include the covariance between the pixels, which will mean that the calculations likely have significant error!

The convolution is defined as:

\[\begin{split}(y\ast g)(x) &= \int_{-\infty}^{\infty} y(X)\ g(\sigma,x-X)\ dX \\ &= \int_{-\infty}^{\infty} \frac{y(X)}{\sqrt{2\pi}\ \sigma(X)}\ \exp\left(-\frac{(x-X)^2}{2\ \sigma(X)^2}\right) dX .\end{split}\]

To minimize edge effects and account for the censoring of the data (finite range in \(x\)), the convolution is actually calculated as a definite integral and normalized as follows:

\[(y\ast g)(x) \sim\frac{ \int_{x-n_\sigma*\sigma(x)}^{x+n_\sigma*\sigma(x)} y(X)\ g(\sigma,x-X)\ dX}{ \int_{x-n_\sigma*\sigma(x)}^{x+n_\sigma*\sigma(x)} g(\sigma,x-X)\ dX} .\]

The above is identical to getting the weighted mean of \(y\) at each position \(x\), where the weights are defined by a Gaussian kernel centered at \(x\) with a variable dispersion.

Use of this function requires:

• y and sigma must be 1D vectors
• y and sigma must be uniformly sampled on the same grid
• sigma must be in pixel units.

Args:

y (np.ndarray): A uniformly sampled function to convolve. sigma (np.ndarray): The standard deviation of the Gaussian

kernel sampled at the same positions as y. The units of sigma must be in pixels.

ye (np.ndarray): (Optional) Errors in the function

\(y(x)\).

Returns:

np.ndarray: Arrays with the convolved function \((y\ast g)(x)\) sampled at the same positions as the input \(x\) vector and its error. The second array will be returned as None if the error vector is not provided.

firefly_instrument.downgrade(wave, flux, deltal_in, sigma_galaxy, wave_instrument, r_instrument)

Adapted from the manga DAP downgrader from Kyle Westfall.

Downgrades an input spectrum to a given galaxy velocity dispersion using the input SEDs resolution and the resolution of the observation.

Returns flux of downgraded SED.

firefly_instrument.match_spectral_resolution(wave, flux, sres, new_sres_wave, new_sres, ivar=None, mask=None, min_sig_pix=0.0, no_offset=True, variable_offset=False, log10=False, new_log10=False)

Adjust the existing spectral resolution of a spectrum to a lower resolution as best as possible. The primary functionality is in spectral_resolution, which determines the Gaussian kernel parameters needed to match the resolution, and convolve_variable_sigma(), which actually performs the convolution to match the resolution.

In particular, see spectral_resolution.GaussianKernelDifference() for a description of how the kernel parameters are determined and how regions where the target resolution is higher than the existing resolution. In this case, one of the options is to adopt an offset of the resolution (in km/s) that could be corrected for in subsequent analysis. In this case, setting variable_offset to True allows the offset to be different for all input spectra. If one expects to combine the spectra, the default behavior should be used, forcing all the spectra to have a constant offset.

Args:

wave (np.ndarray): A 1D or 2D (:math:`N_{rm spec}times

N_{rm pix}`) array with the wavelength in angstroms for a set of spectra. The sampling may be either in linear steps of wavelength or \(\log_{10}\) steps, as set using log10.

flux (np.ndarray): A 1D or 2D (:math:`N_{rm spec}times

N_{rm pix}`) array with the flux sampled at the provided wavelengths.

sres (np.ndarray): A 1D or 2D (:math:`N_{rm spec}times

N_{rm pix}`) array with the spectral resolution, \(R\), at the provided wavelengths.

new_sres_wave (np.ndarray): A 1D vector with the wavelength

in angstroms at which the new resolution of the input spectra has been sampled. The sampling may be either in linear steps of wavelength or \(\log_{10}\) steps, as set using new_log10.

new_sres (np.ndarray): A 1D vector with the new resolution

for the input spectra.

ivar (np.ndarray): (Optional) A 1D or 2D (:math:`N_{rm

spec}times N_{rm pix}`) array with the inverse variance of the flux sampled at the provided wavelengths. This vector is used to estimate the noise in the resolution-matched spectra.

Warning

The accuracy of the errors still remain untested!

mask (np.ndarray): (Optional) A 1D or 2D (:math:`N_{rm

spec}times N_{rm pix}`) array with a uint mask for the flux sampled at the provided wavelengths.

no_offset (bool): (Optional) Force :math:`sigma^2_{v,o} =

0` by masking regions with \(\sigma_{p,d} < -\epsilon_\sigma\); i.e., the value of this arguments selects Option 1 (True) or Option 2 (False). See spectral_resolution.GaussianKernelDifference().

min_sig_pix (float): (Optional) Minimum value of the

standard deviation in pixels allowed before assuming the kernel is a Delta function.

variable_offset (bool): (Optional) Flag to allow the offset

applied to each spectrum (when the input contains more than one spectraum) to be tailored to each spectrum. Otherwise (variable_offset=False) the offset is forced to be the same for all spectra.

log10 (bool): (Optional) Flag that the spectrum has been

binned logarithmically (base 10) in wavelength

new_log10 (bool): (Optional) Flag that the coordinates of

the new spectral resolution are spectrum as been binned logarithmically (base 10) in wavelength.

Returns:

np.ndarray: Five objects are returned:

• A 1D or 2D (\(N_{\rm spec}\times N_{\rm pix}\)) array with the resolution-matched flux sampled at the input wavelengths.
• A 1D or 2D (\(N_{\rm spec}\times N_{\rm pix}\)) array with the spectral resolution, \(R\), of the resolution-matched spectra at the provided wavelengths.
• A 1D vector with any constant offset in resolution in km/s between the targetted value and the result. See spectral_resolution.GaussianKernelDifference().
• A 1D or 2D (\(N_{\rm spec}\times N_{\rm pix}\)) array with a uint mask for the resolution-matched flux sampled at the input wavelengths. This is returned regardless of whether an input mask was provided. Any pixel that had a resolution that was lower than the target resolution (up to some tolerance defined by min_sig_pix) is returned as masked.
• A 1D or 2D (\(N_{\rm spec}\times N_{\rm pix}\)) array with the inverse variance of the resolution-matched flux sampled at the input wavelengths. If ivar is not provided, a ‘None’ returned as the last element

Raises:

ValueError: Raised if:

• the input wave array is 2D and the sres array is not; a 1D wavelength array is allowed for a 2D sres array but not vice versa
• the number of spectral pixels in wave, flux, and sres is not the same
• the shape of the flux, mask (if provided), and ivar (if provided) are not the same
• the shape of the new_sres_wave and new_sres arrays are not the same and/or not 1D

firefly_instrument.resample_vector(y, xRange=None, inLog=False, newRange=None, newpix=None, newLog=True, dx=None, base=10.0, ext_value=0.0, conserve=False, flat=True)

Resample the provided vector to a new grid using integration.

This is a generalization of the routine log_rebin() provided by Michele Cappellari in the pPXF package.

Args:

y (np.ndarray): Vector to resample. Must be 1-D. xRange (array): (Optional) A two-element array with the

starting and ending value for the coordinates of the centers of the first and last pixels in y. If not provided, the pixel coordinates are used; i.e., xRange = [0,y.size-1].

inLog (bool): (Optional) Flag that the input vector is

logarithmically spaced within xRange. Cannot be used if xRange is not provided!

newRange (array): (Optional) Coordinates for the centers of

the first and last pixel in the output vector. If not provided, assumed to be the same as the input range.

newpix (int): (Optional) Number of pixels for the output

vector. If not provided, assumed to be the same as the input vector.

newLog (bool): (Optional) The output vector should be

logarithmically binned in the x-coordinates.

dx (float): (Optional) The sampling step for the output

vector. If newLog is True, this has to be the change in the logarithm of x for the output vector! If not provided, the sampling is set by the output range (see newRange above) and number of pixels (see newpix above).

base (float): (Optional) When logarithmically binning the

output vector, use this as the base. The default is 10.0; use np.exp(1) for natural logarithm.

ext_value (float): (Optional) Set extrapolated values to the

provided float.

conserve (bool): (Optional) Conserve the integral of the

input vector. For example, if the input vector is a spectrum in flux units, you should conserve the flux in the resampling; if the spectrum is in units of flux density, you do not want to conserve the integral.

flat (bool): (Optional) Assume the ‘true’ y function is flat

across a pixel, as is done in M. Cappellari’s log_rebin routine; this is the default behavior. If set to False, the integration follows a basic linear interpolation across the pixel.

Returns:

np.ndarray: Two numpy arrays with the new x coordinates and new y values for the resampled vector.

Raises:

ValueError: Raised if y is not of type np.ndarray, if y

is not one-dimensional, or if xRange is not provided and the input vector is logarithmically binned (see inLog above).

firefly_instrument.resample_vector_npix(outRange=None, dx=None, log=False, base=10.0, default=None)

Determine the number of pixels needed to resample the vector.

Args:

outRange (list or np.ndarray) : Two-element array with the

starting and ending x coordinate of the pixel centers to divide into pixels of a given width. If log is True, this must still be the linear value of the x coordinate, not log(x)!.

dx (float) : Linear or logarithmic pixel width. log (bool) : Flag that the range should be logarithmically

binned.

base (float) : Base for the logarithm default (int) : Default number of pixels to use. The default is

returned if either outRange or dx are not provided.

Returns:

int, np.ndarray: Returns two objects: The number of pixels to cover outRange with pixels of width dx and the adjusted range such that number of pixels of size dx is the exact integer.

Raises:

ValueError: Raised if the range is not a two-element vector

firefly_instrument.spectral_coordinate_step(wave, log=False, base=10.0)

Return the sampling step for the input wavelength vector.

If the sampling is logarithmic, return the change in the logarithm of the wavelength; otherwise, return the linear step in angstroms.

Args:

wave (np.ndarray): Wavelength coordinates of each spectral

channel in angstroms.

log (bool): (Optional) Input spectrum has been sampled

geometrically.

base (float): (Optional) If sampled geometrically, the

sampling is done using a logarithm with this base. For natural logarithm, use np.exp(1).

Returns:

float: Spectral sampling step in either angstroms (log=False) or the step in log(angstroms).

class firefly_instrument.spectral_resolution(wave, sres, log10=False, interp_ext='extrapolate')

Container class for the resolution, \(R = \lambda/\Delta\lambda\), of a spectrum. The primary functionality is to determine the parameters necessary to match the resolution of one spectrum to another. It can also be used as a function to interpolate the spectral resolution at a given wavelenth.

Args:

wave (np.ndarray): A 1D vector with the wavelength in

angstroms. The sampling may be either in linear steps of wavelength or \(\log_{10}\) steps.

sres (np.ndarray): A 1D vector with the spectral resolution,

\(R\), sampled at the positions of the provided wavelength vector.

log10 (bool): (Optional) Flag that the spectrum has been

binned logarithmically (base 10) in wavelength

interp_ext (int or str): (Optional) The value to pass as

ext to the interpolator, which defines its behavior when attempting to sample the spectral resolution beyond where it is defined. See scipy.interpolate.InterpolatedUnivariateSpline. Default is to extrapolate.

Raises:

ValueError: Raised if wave is not a 1D vector or if wave and

sres do not have the same shape.

Attributes:

interpolator

(scipy.interpolate.InterpolatedUnivariateSpline): An object used to interpolate the spectral resolution at any given wavelength. The interpolation is hard-wired to be linear and its extrapolation behavior is defined by interp_ext. The wavelength and resolution vectors are held by this object for later reference if needed.

log10 (bool): Flag that the spectrum has been binned

logarithmically (base 10) in wavelength

cnst (mangadap.util.constants): Used to define the

conversion factor between a Gaussian sigma and FWHM

c (float): The speed of light; provided by astropy.constants. dv (float): The velocity step per pixel in km/s. Defined using

spectrum_velocity_scale() if log10 is True; otherwise set to None.

dw (float): The wavelength step per pixel in angstroms. Defined

as the wavelength step between the first two pixels if log10 is False; otherwise set to None.

min_sig (float): Minimum standard deviation allowed for the

kernel used to match two spectral resolutions. See GaussianKernelDifference().

sig_pd (np.ndarray): The standard deviation in pixels

required to match the spectral resolution of this object to the resolution defined by a different spectral_resolution object. See GaussianKernelDifference().

sig_mask (np.ndarray): A uint vector used to identify

measurements of sig_pd that should not be used to match the spectral resolution of this object to the resolution defined by a different spectral_resolution object. See GaussianKernelDifference().

sig_vo (float): A constant offset of the kernal standard

deviation in km/s that has been applied to sig_pd. See GaussianKernelDifference().

Warning

The default behavior of the interpolator is to extrapolate the input spectral resolution vector when trying to sample from regions beyond where it is sampled. Use interp_ext change this; see scipy.interpolate.InterpolatedUnivariateSpline.

GaussianKernelDifference(new_sres, no_offset=True, min_sig_pix=0.0)

Determine the parameters of a Gaussian kernel required to convert the resolution of this object to the resolution of a different the spectral_resolution object, new_sres.

The spectral resolution is defined as \(R = \lambda/\Delta\lambda\), where \(\Delta\lambda\) is the FWHM of the spectral resolution element. The standard deviation of the resolution element in angstroms is then

\[\sigma_\lambda = \frac{\lambda}{f R}, \ \ {\rm where} \ \ f = \frac{{\rm FWHM_\lambda}}{\sigma_\lambda}.\]

Assuming a Gaussian (in angstroms) line-spread function:

\[\sigma^2_{\lambda,2} = \sigma^2_{\lambda,1} + \sigma^2_{\lambda,d}\]

such that

\[\sigma^2_{\lambda,d} = \left(\frac{\lambda}{f}\right)^2 (R^{-2}_2 - R^{-2}_1)\]

is the defining parameter of the Gaussian kernel needed to take a spectrum of resolution \(R_1\) to one with a resolution of \(R_2\).

For input to convolution_variable_sigma(), the wavelength-dependent parameter of the Gaussian kernel is converted to pixels. This function allows for the input spectra to be linearly sampled in angstroms or log10(angstroms). For the former (log10=False),

\[\sigma^2_{p,d} = \left(\frac{\lambda}{f\ \delta\lambda}\right)^2 (R^{-2}_2 - R^{-2}_1)\]

where \(\delta\lambda\) is the size of the pixel in angstroms. If the units are log10(angstrom) (log10=True), we approximate the velocity width of each pixel to be \(\delta v = c \ln(10.0) (\log\lambda[1]-\log\lambda[0])\), such that

\[\begin{split}\sigma^2_{p,d} &= \left(\frac{c}{ \delta v \lambda}\right)^2 \sigma^2_{\lambda,d} \\ &= \left(\frac{c}{ f\ \delta v}\right)^2 (R^{-2}_2 - R^{-2}_1)\ ;\end{split}\]

\(c\) is the speed of light in km/s.

The nominal use of this algorithm assumes \(R_1 \geq R_2\). However, in practice, convolution_variable_sigma() only uses a Gaussian kernel up to some minimum value of \(\epsilon_\sigma\); below this, the kernel is assumed to be a Delta function. Therefore, as long as

\[\sigma_{p,d} \equiv \sigma^2_{p,d}/\sqrt{|\sigma^2_{p,d}|} \geq -\epsilon_\sigma\ ,\]

the behavior of convolution_variable_sigma() should not be affected.

Even so, there may be spectral regions that do not have \(\sigma_{p,d} \geq -\epsilon_\sigma\); for such spectral regions there are three choices:

(Option 1) trim the spectral range to only those spectral regions where the existing resolution is better than the target resolution,

(Option 2) match the existing resolution to the target resolution up to some constant offset that must be accounted for in subsequent analyses, or

(Option 3) allow for a wavelength dependent difference in the spectral resolution that must be accounted for in subsequent analyses.

The choice of either Option 1 or 2 is selected by setting no_offset to, respectively, True or False; Option 1 is the default behavior. Currently, Option 3 is not allowed.

For Option 1, pixels with \(\sigma_{p,d} < -\epsilon_\sigma\) are masked (sigma_mask = 1); however, the returned values of \(\sigma_{p,d}\) are left unchanged.

For Option 2, we define

\[\sigma^2_{v,o} = -{\rm min}(\sigma^2_{v,d}) - {\rm max}(\epsilon_\sigma \delta v)^2\]

where \(\delta v\) is constant for the logarithmically binned spectrum and is wavelength dependent for the linearly binned spectra; in the latter case, the velocity step is determined for each pixel:

_wave = self.wave()
dv = self.c * (2.0*(_wave[1:] - _wave[0:-1]) / (_wave[1:] + _wave[0:-1]))

If \(\sigma^2_{v,o} > 0.0\), it must be that \({\rm min}(\sigma^2_{v,d}) < -{\rm max}(\epsilon_\sigma \delta v)^2\), such that an offset should be applied. In that case, the returned kernel parameters are

\[\sigma^\prime_{v,d} = \sqrt{\sigma^2_{v,d} + \sigma^2_{v,o}}\ .\]

with the units converted to pixels using the equations above, no pixels are masked, and \(\sqrt{\sigma^2_{v,o}}\) is returned for the offset. Otherwise, the offset is set to 0.

Args:

new_sres (spectral_resolution): Spectral resolution

to match to.

no_offset (bool): (Optional) Force :math:`sigma^2_{v,o}

= 0` by masking regions with \(\sigma_{p,d} < -\epsilon_\sigma\); i.e., the value of this arguments selects Option 1 (True) or Option 2 (False).

min_sig_pix (float): (Optional) Minimum value of the

standard deviation allowed before assuming the kernel is a Delta function.

adjusted_resolution(indx=None)

Return the resolution that should result from applying convolution_variable_sigma() to the spectrum associated with this spectral resolution object using sigma_pd. I.e., calculate:

\[R_2 = \left[ \left(\frac{f}{c}\right)^2 \sigma^2_{v,d} + R^{-2}_1\right]^{-1/2}\ . \]

Args:

indx (tuple): (Optional) Selection tuple used to return

a subset of the full resolution vector.

Returns:

np.ndarray: The (full or selected) vector with the adjusted resolution.

match(new_sres, no_offset=True, min_sig_pix=0.0)

Currently only an alias for GaussianKernelDifference().

offset_GaussianKernelDifference(offset)

If the properties required to match the resolution of one spectrum to another has already been calculated (see GaussianKernelDifference()), this allows for one to apply an additional offset. The additional offset must be in km/s (not pixels).

The offset is applied in quadrature; however, the offset can be negative such that one can reduce sig_pd. Once converted to km/s, the offset is applied by calculating:

\[\sigma^{\prime\ 2}_{v,d} = \sigma^{2}_{v,d} + \sigma_{off}|\sigma_{off}|\ .\]

sig_vo is adjusted in the same way, and the change in \(\sigma^{\prime\ 2}_{v,d}\) is then propagated to sig_pd and sig_mask.

Args:

offset (float): Value of the standard deviation in km/s to

add in quadrature to a previously calculated sig_pd.

Raises:

ValueError: Raised if the kernel properties have not yet

been defined.

sres()

Return the resolution vector; held by interpolator.

wave()

Return the wavelength vector; held by interpolator.

firefly_instrument.spectrum_velocity_scale(wave)

Determine the velocity sampling of an input wavelength coordinate vector.

Note

The wavelength vector is assumed to be geometrically sampled! However, the input units expected to be in angstroms, not, e.g., log(angstrom).

Args:

wave (np.ndarray): Wavelength coordinates of each spectral

channel in angstroms. It is expected that the spectrum has been sampled geometrically

Returns:

float: Velocity scale of the spectrum in km/s.

firefly_instrument.where_not(indx, size)

Return a tuple with the indices of a vector that were not selected by a call to `np.where`_. The function currently only works for 1D vectors.

Args:

indx (tuple): Tuple returned by a call to `np.where`_ for a

1D vector.

size (int): Length of the original vector in the call to

`np.where`_.

Warning

Performs no checks of the input.