Contents:
Suite to reduce spectroscopic data.
clean – clean and replace bad pixels extract
trace – trace spectral orders makeprofile – compute mean spectral PSF (a spline) for an order fitprofile – fit given spline-PSF to a spectral cross-section
| INPUTS: |
|
|---|
| NAME: |
EXAMPLE:
|
|---|
pixel-to-pixel nonuniformities (i.e., traditional flat-fielding)
detector nonlinearities
tilted spectral lines
non-uniform slit widths (which cause non-smooth backgrounds)
Note that the dispersion direction should be ‘horizontal’ (i.e., parallel to rows) in this frame.
| INPUTS: |
|
|---|
Row-index in each subregion of the location of the spectral-trace-of-interest. Only used for rectifying of sky frame; leaving this at None will have at most a mild effect. If not None, should be the same length as subreg_corners. If subreg_corners[0] = [800, 950] then locs[0] might be set to, e.g., 75 if the trace lies in the middle of the subregion.
- gain: scalar
- Detector gain, in electrons per ADU
- readnoise: scalar
- Detector readnoise, in electrons
| EXAMPLE: | try:
from astropy.io import fits as pyfits
except:
import pyfits
import spec
import ir
obs = ir.initobs('20121010')
allinds = [[7, 294], [310, 518], [532, 694], [710, 960], [976, 1360], [1378, 1673], [1689, 2022]]
locs = [221, 77, 53, 88, 201, 96, 194]
unnorm_flat_fn = obs._raw + 'mudflat.fits'
if False:
sky = pyfits.getdata(obs._raw + 'masktersky.fits')
unnorm_domeflat = pyfits.getdata(unnorm_flat_fn)
skycorrect, pixcorrect = spec.make_spectral_flats(sky, unnorm_domeflat, allinds, badpixelmask=obs.badpixelmask, locs=locs)
else:
skycorrect = obs._raw + 'skycorrect.fits'
pixcorrect = obs._raw + 'pixcorrect.fits'
linearcoef='/Users/ianc/proj/transit/data/mosfire_linearity/linearity/mosfire_linearity_cnl_coefficients_bic-optimized.fits'
rawfn = obs.rawsci
outfn = obs.procsci
spec.calibrate_stared_mosfire_spectra(rawfn, outfn, skycorrect, pixcorrect, allinds, linearize=True, badpixelmask=obs.badpixelmask, locs=locs, polycoef=linearcoef, verbose=True, clobber=True, unnormalized_flat=unnorm_flat_fn)
|
|---|---|
| NOTES: | This routine is slow, mostly because of the call to defringe_sinusoid(). Consider running multiple processes in parallel, to speed things up! |
| SEE_ALSO: |
Use simple fitting to subtract fringes and sky background.
Which spatial rows (if dispaxis=0) to use when fitting the tilt of sky lines across the spectrum. If you want to use all, set to None. If you want to ignore some (e.g., because there’s a bright object’s spectrum there) then set those rows’ elements of spatial_index to ‘False’.
| NOTES: | Note that (in my experience!) this approach works better when you set all weights to unity, rather than using the suggested (photon + read noise) variances. |
|---|---|
| REQUIREMENTS: | Numerical analysis routines (for analysis.fmin()) Planetary phase curve routines (for phasecurves.errfunc()) lomb (for Lomb-Scargle periodograms) SciPy ‘signal’ module (for median-filtering) |
Returns the best-fit sky frame.
Extract spectrum
| INPUTS: |
|
|---|---|
| OPTIONS: |
|
| RETURNS: |
|
| EXAMPLE: |
|
| SEE_ALSO: | |
| NOTES: | Note that this is non-optimal with highly tilted or curved spectra, for the reasons described by Marsh (1989) and Mukai (1990). |
Fit two Gaussians simultaneously to an input data vector.
| INPUTS: |
|
|---|
SEE ALSO: analysis.gaussian(), fitGaussian()
Fit a Gaussian function to an input data vector.
Return the fit, and uncertainty estimates on that fit.
SEE ALSO: analysis.gaussian()
Fit a Gaussian function to an input data vector.
Return the fit, and uncertainty estimates on that fit.
SEE ALSO: analysis.gaussian()
Helper function to fit 1D PSF near a given region. Assumes spectrum runs horizontally across the frame!
Fit a given spatial profile to a spectrum
Helper function for extractSpectralProfiles()
Fit a 1D tophat function to an input data vector.
Return the fit, and uncertainty estimates on that fit.
SEE ALSO: analysis.gaussian()
Fit a spline-class spatial profile to a spectrum cross-section
| PURPOSE: | Compute the integral from -inf to x of the normalized Gaussian |
|---|---|
| INPUTS: |
|
| NOTES: | Designed to copy the IDL function of the same name. |
humidpressure
To convert relative humidity into a H2O vapor partial pressure
Spectroscopy
humidpressure(RH, 273.15)
RH - relative humidity, in percent T - temperature, in Kelvin
h2o_pp : water vapor partial pressure, in Pascals
As outlined in Butler (1998): “Precipitable Water at KP”, MMA Memo 238 (which refers in turn to Liebe 1989, “MPM - An Atmospheric Millimeter-Wave Propagation Model”). Liebe claims that this relation has an error of <0.2% from -40 C to +40 C.
None
2011-10-08 17:08 IJMC: Created.
lightloss
To determine the slit losses from a spectrum.
Spectroscopy
### TBD lightloss, obj, std, wguide, seeing, out, CANCEL=cancel
obj - FITS file of the object spectrum wguide - wavelength at which guiding was done seeing - seeing FWHM at the guiding wavelength
press - mm Hg typical value (615 for IRTF, unless set) water - mm Hg typical value (2 for IRTF, unless set) temp - deg C typical value (0 for IRTF, unless set) fco2 - relative concentration of CO2 (0.004, unless set) wobj - wavelength scale for data dx - horizontal offset of star from slit center dy - vertical offset of star from slit center retall- whether to return much diagnostic info, or just lightloss.
‘seeing’, ‘dx’, and ‘dy’ should all be in the same units, and also the same units used to define the slit dimensions in the obj FITS file header
CANCEL - Set on return if there is a problem
array : fractional slit loss at each wavelength value
tuple of arrays: (slitloss, disp_obj, diff, fwhm, dx_obj, dy_obj)
Reads a Spextool FITS file.
None
2003-10-21 - Written by W D Vacca 2011-10-07 20:19 IJMC: Converted to Python, adapted for single objects 2011-10-14 14:01 IJMC: Added check for Prism mode (has
different slit dimension keywords) and different pyfits header read mode.
2011-11-07 15:53 IJMC: Added ‘retall’ keyword
lightloss2
To determine the slit losses from an observation (no FITS file involved)
Spectroscopy
### TBD lightloss, obj, std, wguide, seeing, out, CANCEL=cancel
wobj - wavelength scale for data slitwd - width of slit, in arcsec slitht - height of slit, in arcsec slitPA - slit Position Angle, in radians targetPA - Parallactic Angle at target, in radians zenith_angle - Zenith Angle, in radians wguide - wavelength at which guiding was done seeing - seeing FWHM at the guiding wavelength
press - mm Hg typical value (615, unless set) water - mm Hg typical value (2 , unless set) temp - deg C typical value (0 , unless set) fco2 - relative concentration of CO2 (0.004, unless set) obsalt- observatory altitude, in km teldiam- observatory limiting aperture diameter, in m dx - horizontal offset of star from slit center dy - vertical offset of star from slit center retall- whether to return much diagnostic info, or just lightloss.
values of wobj. This should be an array of the same size as wobj, with zero corresponding to the vertical middle of the slit and positive values tending toward zenith. In this case xdisp will be computed as XXXX rather than from the calculated atmospheric dispersion; dx and dy will also be ignored.
all values of wobj. This should be an array of the same size as wobj, measured in arc seconds.
‘slitwidth’, ‘slitheight’, ‘seeing’, ‘dx’, ‘dy’, and ‘fwhm’ (if used) should all be in the same units: arc seconds.
array : fractional slit loss at each wavelength value
tuple of arrays: (slitloss, disp_obj, diff, fwhm, dx_obj, dy_obj)
All input-driven. For the SpeXTool-version analogue, see lightloss()
import numpy as np import spec w = np.linspace(.5, 2.5, 100) # Wavelength, in microns d2r = np.deg2rad(1.) #targetPA, za = spec.parangle(1.827, 29.67*d2r, lat=20.*d2r) targetPA, za = 105.3, 27.4 slitPA = 90. * d2r
spec.lightloss2(w, 3., 15., slitPA, targetPA*d2r, za*d2r, 2.2, 1.0)
2003-10-21 - Written by W D Vacca 2011-10-07 20:19 IJMC: Converted to Python, adapted for single objects 2011-10-14 14:01 IJMC: Added check for Prism mode (has
different slit dimension keywords) and different pyfits header read mode.
2011-11-07 15:53 IJMC: Added ‘retall’ keyword 2011-11-07 21:17 IJMC: Cannibalized from SpeXTool version 2011-11-25 15:06 IJMC: Added ydisp and fwhm options.
Generate a normalized Sky frame from SpeX slitless spectroscopy data.
| INPUTS: |
|
|---|---|
| OUTPUTS: | (sky, skyHeader) |
Construct appropriate corrective frames for multi-object spectrograph data. Specifically: corrections for irregular slit widths, and pixel-by-pixel detector sensitivity variations.
Master spectral sky frame, e.g. from median-stacking many sky frames or masking-and-stacking dithered science spectra frames. This frame is used to construct a map to correct science frames (taken with the identical slit mask!) for irregular sky backgrounds resulting from non-uniform slit widths (e.g., Keck/MOSFIRE).
Note that the dispersion direction should be ‘horizontal’ (i.e., parallel to rows) in this frames.
Master dome spectral flat, e.g. from median-stacking many dome spectra. This need not be normalized in the dispersion or spatial directions. This frame is used to construct a flat map of the pixel-to-pixel variations in detector sensitivity.
Note that the dispersion direction should be ‘horizontal’ (i.e., parallel to rows) in this frames.
| RETURNS: | wideslit_skyflat, narrowslit_domeflat |
|---|---|
| EXAMPLE: | try:
from astropy.io import fits as pyfits
except:
import pyfits
import spec
import ir
obs = ir.initobs('20121010')
sky = pyfits.getdata(obs._raw + 'masktersky.fits')
domeflat = pyfits.getdata(obs._raw + 'mudflat.fits')
allinds = [[7, 294], [310, 518], [532, 694], [710, 960], [976, 1360], [1378, 1673], [1689, 2022]]
locs = [221, 77, 53, 88, 201, 96, 194]
skycorrect, pixcorrect = spec.make_spectral_flats(sky, domeflat, allinds, obs.badpixelmask, locs=locs)
|
Helper function for wavelengthMatch(): generate a scaled, interpolative model of the template.
Make a spatial profile from a spectrum, given its traced location. We interpolate the PSF at each pixel to a common reference frame, and then average them.
a spline-function that interpolates pixel locations onto the mean profile
a stack of data slices
estimates of the uncertainties
good pixel flag
list of splines
Helper function for XXXX.
| INPUTS: |
|
|---|---|
| NOTES: | Helper function for functions XXXX |
Model a raw spectral trace!
| INPUTS: |
|---|
NOTE that most inputs should be in the _rectified_ frame.
Trace background pedestal level : 1D array
Width of background pedestarl level : scalar (for now)
Center of trace : 1D array
Offset of object spectrum, relative to center : scalar
width of 1D PSF : scalar
Area of 1D psf : 1D array
Distortion (x and y, somehow???)
Scattered light background : scalar (for now)
Model a spectral resolution element.
| INPUTS: | param : sequence
|
|---|---|
| OUTPUTS: |
|
| DESCRIPTION: | Model a spectral resolution element along the spatial direction. This consists of a (presumably Gaussian) line profile superimposed on the spectral trace’s top-hat-like background, with an additional constant (or polynomial) out-of-echelle-order background component. |
Generate a rotational profile, convolve it with a second input profile, normalize it (simply), and return.
| INPUTS: |
dv : velocity grid prof2 : second input profile |
|---|
Trace and normalize a spectroscopic flat field frame.
| INPUTS: |
|
|---|
Extract spectrum, following Horne 1986.
| INPUTS: |
|
|---|---|
| OPTIONS: |
|
| RETURNS: |
|
| EXAMPLE: | |
| SEE_ALSO: | |
| NOTES: | Horne’s classic optimal extraction algorithm is optimal only so long as the spectral traces are very nearly aligned with detector rows or columns. It is not well-suited for extracting substantially tilted or curved traces, for the reasons described by Marsh 1989, Mukai 1990. For extracting such spectra, see superExtract(). |
Run optimal spectral extraction in IDL; pass results to Python.
| INPUTS: |
|
|---|---|
| OPTIONS: |
|
| OUTPUTS: |
|
| NOTES: | Note that this more closely follows Horne et al. than does optimalExtract(), and is faster than both that function and (especially!) extractSpectralProfiles(). The only downside (if it is one) is that this function requires IDL. |
| TO-DO: | Add options for user input of a variance frame, or of sky variance. Allow more flexibility (tracing, input/output options, etc.) |
| REQUIREMENTS: | IDL |
- NAME:
parangle
- PURPOSE:
To compute the parallactic angle at a given position on the sky.
- CATEGORY:
Spectroscopy
- CALLING SEQUENCE:
eta, za = parangle(HA, DEC, lat)
- INPUTS:
HA - Hour angle of the object, in decimal hours (0,24) DEC - Declination of the object, in degrees lat - The latitude of the observer, in degrees
- KEYWORD PARAMETERS:
CANCEL - Set on return if there is a problem
- OUTPUTS:
eta - The parallactic angle za - The zenith angle
- PROCEDURE:
Given an objects HA and DEC and the observers latitude, the zenith angle and azimuth are computed. The law of cosines then gives the parallactic angle.
- EXAMPLE:
NA
- MODIFICATION HISTORY:
2000-04-05 - written by M. Cushing, Institute for Astronomy,UH 2002-08-15 - cleaned up a bit. 2003-10-21 - changed to pro; outputs zenith angle as well - WDV 2011-10-07 17:58 IJMC: Converted to Python
| INPUTS: | ax : (axes instance) – axes in which to pick a location
|
|---|
Compute the chi-squared error on a spectrum vs. a profile
Convert GMOS frames taken with custom ROIs into standard FITS frames.
| INPUTS: |
|
|---|---|
| OPTIONS: |
|
Resample a spectrum while conserving flux density.
| INPUTS: |
|
|---|---|
| NOTE: | Format is the same as numpy.interp() |
| REQUIREMENTS: |
Compute the rotational profile of a star, assuming solid-body rotation and linear limb darkening.
This uses Eq. 18.14 of Gray’s Photospheres, 2005, 3rd Edition.
| INPUTS: | delta_epsilon : 2-sequence
|
|---|---|
| EXAMPLE: | import pylab as py
import spec
dlam = py.np.linspace(-2, 2, 200) # Create wavelength grid
profile = spec.rotationalProfile([1, 0.6], dlam)
py.figure()
py.plot(dlam, profile)
|
Run LBLRTM to compute atmospheric transmittance and/or radiance.
| INPUTS: |
|
|---|---|
| OPTIONS: |
|
| OUTPUTS: | A 2- or 3-tuple: First element is wavelength in microns, second element is transmittance (if requested). Radiance will (if requested) always be the last element, and in f_nu units: W/cm2/sr/(cm^-1) |
| REQUIREMENTS: | SciPy OCTAVE or MATLAB D. Feldman’s set of MATLAB wrapper scripts |
Use cross-correlation to subtract tilted sky backgrounds.
| RETURNS: | a model of the sky background, of the same shape as ‘subframe.’ |
|---|
Use difference-imaging techniques to subtract moderately tilted sky background. Doesn’t work so well!
| NOTES: | Note that (in my experience!) this approach works better when you set all weights to unity, rather than using the suggested (photon + read noise) variances. |
|---|
Returns the best-fit sky frame.
Use PCA and blank sky frames to subtract
f0 = pyfits.getdata(odome.procsci[0]) mask = pyfits.getdata(odome._proc + ‘skyframes_samplemask.fits’).astype(bool) badpixelmask = pyfits.getdata( odome.badpixelmask).astype(bool)
The simplest way to fit sky to a ‘frame’ containing bright spectral is to include the spectral-trace regions in ‘mask’ but set the ‘variance’ of those regions extremely high (to de-weight them in the least-squares fit).
To use for multi-object data, consider running multiple times (once per order)
Returns the best-fit sky frame as determined from the first ‘npca’ PCA components.
| NAME: | slittrans |
|---|---|
| PURPOSE: | Compute flux passing through a slit assuming a gaussian PSF. |
| CATEGORY: | Spectroscopy |
| CALLING SEQUENCE: | |
result = slittrans(width,height,fwhm,xoffset,yoffset,CANCEL=cancel) |
|
| INPUTS: | width - Width of slit. height - Height of slit. fwhm - Full-width at half-maximum of the gaussian image. xoffset - Offset in x of the image from the center of the slit. yoffset - Offset in y of the image from the center of the slit.
|
result = slittrans(0.3,15,0.6,0,0)
Computes the fraction of the flux transmitted through a slit 0.3x15 arcseconds with a PSF of 0.6 arcseconds FWHM. The PSF is centered on the slit.
Based on M Buie program, 1991 Mar., Marc W. Buie, Lowell Observatory Modified 2000 Apr., M. Cushing to include y offsets. 2011-10-07 15:45 IJMC: Converted to Python 2011-11-14 16:29 IJMC: Rewrote to use erf() rather than
Fix scattered light in SpeX/SXD K-band and write a new file.
| INPUTS: |
|
|---|---|
| OPTIONS: |
Other options will be passed to spexsxd_scatter_model() |
| OUTPUTS: |
|
Model the scattered light seen in SpeX/SXD K-band frames.
| INPUTS: |
|
|---|---|
| OPTIONS: |
|
| OUTPUT: |
OR:
|
| REQUIREMENTS: | pyfits, numpy, fit_atmo, Numerical analysis routines, Planetary phase curve routines |
| TO_DO_LIST: | I could stand to be more clever in modeling the scattered light components – perhaps fitting for the width, or at least allowing the width to be non-integer. |
Extract a spectrum from a frame using one of several methods.
| INPUTS: |
|
|---|---|
| RETURNS: | spectrum, error or variance of spectrum, sky background, ... |
| NOTES: | When ‘optimalextract’ is used: if len(bkg_radii)==2 then the background apertures will be reset based on the median location of the trace. If extract_radius is a singleton, it will be similarly redefined. |
| EXAMPLE: | frame = pyfits.getdata('spectral_filename.fits')
gain, readnoise = 3.3, 30
output = spec.spextractor(frame, gain, readnoise, mode='superExtract', options=dict(bord=2, bkg_radii=[20, 30], extract_radius=15, polyspacing=1./3, pord=5, verbose=True, trace=trace, qmode='slow-linear'))
output2 = spec.spextractor(frame, gain, readnoise, mode='optimalExtract', options=dict(bkg_radii=[20,30], extract_radius=15, bord=2, bsigma=3, pord=3, psigma=8, csigma=5, verbose=1))
|
Optimally extract curved spectra, following Marsh 1989.
| INPUTS: |
|
|---|---|
| OPTIONS: |
|
| RETURNS: |
|
| EXAMPLE: | import spec
import numpy as np
import pylab as py
def gaussian(p, x):
if len(p)==3:
p = concatenate((p, [0]))
return (p[3] + p[0]/(p[1]*sqrt(2*pi)) * exp(-(x-p[2])**2 / (2*p[1]**2)))
# Model some strongly tilted spectral data:
nx, nlam = 80, 500
x0 = np.arange(nx)
gain, readnoise = 3.0, 30.
background = np.ones(nlam)*10
flux = np.ones(nlam)*1e4
center = nx/2. + np.linspace(0,10,nlam)
FWHM = 3.
model = np.array([gaussian([flux[ii]/gain, FWHM/2.35, center[ii], background[ii]], x0) for ii in range(nlam)])
varmodel = np.abs(model) / gain + (readnoise/gain)**2
observation = np.random.normal(model, np.sqrt(varmodel))
fitwidth = 60
xr = 15
output_spec = spec.superExtract(observation, varmodel, gain, readnoise, polyspacing=0.5, pord=1, bkg_radii=[10,30], extract_radius=5, dispaxis=1)
py.figure()
py.plot(output_spec.spectrum.squeeze() / flux)
py.ylabel('(Measured flux) / (True flux)')
py.xlabel('Photoelectrons')
|
| TO_DO: | Iterate background fitting and reject outliers; maybe first time would be unweighted for robustness. Introduce even more array-based, rather than loop-based, calculations. For large spectra computing the C-matrix takes the most time; this should be optimized somehow. |
| SEE_ALSO: |
Grey-pixel tophat function with set width param: [cen_pix, amplitude, background] newparam: [amplitude, full width, cen_pix, background] x : must be array of ints, arange(0, size-1) returns the model.
Standard tophat function (alternative version).
| INPUTS: |
|
|---|---|
| OUTPUTS: |
Trace spectral orders for a specified filename.
OPTIONAL INPUTS: pord : int
polynomial order of spectral order fit
full path and filename to a 2D uncertainties FITS file, _OR_ a 2D numpy array representing such a file.
If this is set, ‘g’ and ‘rn’ below are ignored. This is useful if, e.g., you are analyzing data which have already been sky-subtracted, nodded on slit, or otherwise altered. But note that this must be the same size as the input data!
| RETURNS: |
|
|---|---|
| NOTES: | If tracing fails, a common reason can be that fitwidth is too small. Try increasing it! |
Determine dispersion solution for a spectrum, from a template.
| INPUTS: |
|
|---|---|
| OPTIONS: |
|
| RETURNS: | (wavelength, wavelength_polynomial_coefficients, full_parameter_set) |
| NOTES: | This implementation uses a rather crude MCMC sampling approach to sample parameters space and ‘home in’ on better solutions. There is probably a way to do this that is both faster and more optimal... Note that if ‘spectrum’ and ‘template’ are of different lengths, the longer one will be trimmed at the end to make the lengths match. |
| REQUIREMENTS: |
Compute zenith angle (in degrees) for an observation.
| INPUTS: |
|
|---|---|
| OUTPUTS: | Zenith angle, in degrees, for the specified observation |
| REQUIREMENTS: | Aperture and simple PSF-fitting photometrymem Numpy |