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 splinePSF to a spectral crosssection
INPUTS: 


NAME: 
EXAMPLE:


pixeltopixel nonuniformities (i.e., traditional flatfielding)
detector nonlinearities
tilted spectral lines
nonuniform slit widths (which cause nonsmooth backgrounds)
Note that the dispersion direction should be ‘horizontal’ (i.e., parallel to rows) in this frame.
INPUTS: 


Rowindex in each subregion of the location of the spectraltraceofinterest. 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_bicoptimized.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 LombScargle periodograms) SciPy ‘signal’ module (for medianfiltering) 
Returns the bestfit sky frame.
Extract spectrum
INPUTS: 


OPTIONS: 

RETURNS: 

EXAMPLE: 

SEE_ALSO:  
NOTES:  Note that this is nonoptimal 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 splineclass spatial profile to a spectrum crosssection
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 MillimeterWave Propagation Model”). Liebe claims that this relation has an error of <0.2% from 40 C to +40 C.
None
20111008 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
20031021  Written by W D Vacca 20111007 20:19 IJMC: Converted to Python, adapted for single objects 20111014 14:01 IJMC: Added check for Prism mode (has
different slit dimension keywords) and different pyfits header read mode.
20111107 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 inputdriven. For the SpeXToolversion 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)
20031021  Written by W D Vacca 20111007 20:19 IJMC: Converted to Python, adapted for single objects 20111014 14:01 IJMC: Added check for Prism mode (has
different slit dimension keywords) and different pyfits header read mode.
20111107 15:53 IJMC: Added ‘retall’ keyword 20111107 21:17 IJMC: Cannibalized from SpeXTool version 20111125 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 multiobject spectrograph data. Specifically: corrections for irregular slit widths, and pixelbypixel detector sensitivity variations.
Master spectral sky frame, e.g. from medianstacking many sky frames or maskingandstacking 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 nonuniform 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 medianstacking 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 pixeltopixel 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 splinefunction 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 tophatlike background, with an additional constant (or polynomial) outofechelleorder 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 wellsuited 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. 
TODO:  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:
20000405  written by M. Cushing, Institute for Astronomy,UH 20020815  cleaned up a bit. 20031021  changed to pro; outputs zenith angle as well  WDV 20111007 17:58 IJMC: Converted to Python
INPUTS:  ax : (axes instance) – axes in which to pick a location


Compute the chisquared 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 solidbody rotation and linear limb darkening.
This uses Eq. 18.14 of Gray’s Photospheres, 2005, 3rd Edition.
INPUTS:  delta_epsilon : 2sequence


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 3tuple: 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 crosscorrelation to subtract tilted sky backgrounds.
RETURNS:  a model of the sky background, of the same shape as ‘subframe.’ 

Use differenceimaging 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 bestfit 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 spectraltrace regions in ‘mask’ but set the ‘variance’ of those regions extremely high (to deweight them in the leastsquares fit).
To use for multiobject data, consider running multiple times (once per order)
Returns the bestfit 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  Fullwidth at halfmaximum 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. 20111007 15:45 IJMC: Converted to Python 20111114 16:29 IJMC: Rewrote to use erf() rather than
Fix scattered light in SpeX/SXD Kband 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 Kband 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 noninteger. 
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='slowlinear'))
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((xp[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 arraybased, rather than loopbased, calculations. For large spectra computing the Cmatrix takes the most time; this should be optimized somehow. 
SEE_ALSO: 
Greypixel 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, size1) 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 skysubtracted, 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 PSFfitting photometrymem Numpy 