Phase curve spectroscopic analysis, and associated subroutines.

Contents:

pcsa.detdata(block, t, mode='poly', ord=1, meansub=False, eblock=None, n=200, verbose=False, polyind=None)

Remove temporal trends from blocks of data.

INPUTS: (L observations of M echelle orders with N wavelength bins)

block: array of data from pcsa.prepdata, either (N x L) or (L x M x N) t: 1D array of times (e.g., HJD), of length (L)

OPTIONAL INPUTS:

mode=’poly’ – method to use to detrend. Either ‘poly’

or ‘pca’( for principal component analysis.)

ord=1 – if mode==’poly’: polynomial order to subtract

– if mode==’pca’: number of Principal Components to subtract

polyind – 1D numpy array of length L. Use only ‘True’ indices

in computing polynomial fits, then remove the fit from the entire data (e.g., to avoid fitting out a transit).

meansub=False – whether to zero-out the mean, or keep the means fixed

eblock=None – 1-sigma errors on the values in ‘block’ for PCA n=200 – number of wavelength bins to use for each PCA calculation

verbose=False – bool flag to spit out lots of verbose messages

OUTPUT:

newblock – detrended dataset

EXAMPLE:

TBW

pcsa.fscs(data)

Perform Swain et al.’s “Fourier Self-Complementary Spectrum” analysis.

Input: Data, numpy array – M x N: M observations of N variables

pcsa.pca(block, eblock, ord=1, verbose=False, trans=False)

Calculate principal component analysis for an N x M input array of N observations of M variables.

(Helped function for detdata, but could eventually be stand-alone)

ord – number of principal components to remove, sorted by eigenvalue size trans – take the transpose of block and eblock (i.e., they are

M x N instead of the desired N x M)

WARNING: do not trust any PC’s with coefficient magnitudes near-zero!

pcsa.pca2(block, eblock=None, ord=1, svd=True)

Compute Principal Component Analysis on a 2D array.

NOTE: I don’t trust the MDP module (or, I don’t understand what it’s doing), since it does not agree with my hand-coded (somewhat slower) PCA routine.

block – N x M numpy array – N observations of M variables.

eblock – dummy option kept for backwards compatibility; does nothing. ord – number of Principal Components to extract svd – bool flag; whether to use SVD to compute PCA

pcsa.pca_eigv(block, ord=None, svd=True, retnode=False)

Compute the eigenvalues of covariance matrix (i.e., the PC coefficients) for input data block.

ord – number of eigenvalues to return. If None, returns all. svd – bool flag; whether to use SVD to compute PCA retnode – also return the mdp.PCANode object

pcsa.pca_project(pc, block)

Project a given principal component eigenvector onto a block of data

pc – length-N 1D numpy array block – N x M input array of N observations of M variables.

pcsa.prepdata(date, trim=None, telcor=True, getcal=False, divop='median', interp=True)

Collect data and prepare for LSD analysis.

INPUTS:

date – string, observation date: ‘2008mar21’ – 55cnce

‘2008jun15’ – tauboob ‘2008jul12’ – 189733, L-band

OPTIONAL INPUTS:

trim – list or 1D array, indices of loaded data to keep, e.g,

[0,1,2,3,...,4000]

telcor=True – load telluric-corrected data getcal=False – load calibration spectra instead of sci. spectra divop=’median’ – by default, divide each spectrum by its

median to normalize it. Can also set to ‘none’

interp=True – whether to load wavelength-interpolated

(upsampled) spectra (if True) or the raw 1024-pixel spectra (if False)

OUTPUT: (assume L observations of M echelle orders with N wavelength bins)

w – N-length wavelength grid for data proc – processed data cubes, L x M x N eproc – computed _instrumental_ uncertainties, L x M x N sky – sky background intensity, L x M x N hjd – L-length time array of heliocentric Julian Date airmass – L-length airmass array p – “planet” object from analysis.getobj() rv – planet’s maximal Radial Velocity at specified HJD keyvals – various specified header keywords at specified HJD spec_location – location of star along slit at given HJD, L x M

pcsa.princomp(A, numpc=0)

computing eigenvalues and eigenvectors of covariance matrix

Taken from http://glowingpython.blogspot.it/2011/07/pca-and-image-compression-with-numpy.html

pcsa.salign(s1, s2=None, retshift=False, offset=None)

Shift and align spectra.

INPUT:

s1 – array. If 1D, a second spectrum (to be aligned) must also

be input. If 2D (of shape [nobs x npix]), the first spectrum is used to align each of the succeeding (nobs-1) spectra.

OPTIONAL INPUT:

retshift – return pixel shifts as second tuple element (e.g.,

for use in aligning related vector sets.)

offset – if None, calculate the aligning offsets. Otherwise,

use these offsets (via SOFFSET) to align the vectors.

OUTPUT:

salign – array. A 1D or 2D array, corresponding to the shape of s1.

NOTES:

• If 1D, s1 and s2 must be the same length.
• Uses cross-correlation and the nearest pixel. The disadvantage can be poor alignment; the advantage is speed.

pcsa.salign2(s1, s2=None, retshift=False, offset=None)

Shift and align spectra.

INPUT:

s1 – array. If 1D, a second spectrum (to be aligned) must also

be input. If 2D (of shape [nobs x npix]), the first spectrum is used to align each of the succeeding (nobs-1) spectra.

OPTIONAL INPUT:

retshift – return pixel shifts as second tuple element (e.g.,

for use in aligning related vector sets.)

offset – if None, calculate the aligning offsets. Otherwise,

use these offsets (via SOFFSET) to align the vectors.

OUTPUT:

salign – array. A 1D or 2D array, corresponding to the shape of s1.

NOTES:

• If 1D, s1 and s2 must be the same length.
• Uses cross-correlation and the nearest pixel. The disadvantage can be poor alignment; the advantage is speed.

pcsa.shiftvec(vec, offset)

Shift a 1D vector by an offset. Helper function for salign.

non-overlapping elements are replaced by zeros.

pcsa.soffset(s1, s2=None, plotalot=False, oversamp=16, searchrange=[-2, 2], retnew=False, maxshift=10, method='interp')

Get pixelwise offsets between two similar 1D vectors.

INPUT:

s1 – array. If 1D, a second vector (to be aligned) must also

be input. If 2D (of shape [nvec x npix]), the first vector is used to align each of the succeeding (nvec-1) spectra.

oversamp – scalar. Factor by which to oversample the inputs, to

permit computation of more precise offsets.

searchrange – list. Minimum and maximum number of pixels on

each side of the maximum cross-correlation value to use for parabolic fitting.

maxshift – int. Min/max of pixel correlation lags to compute.

method – str: ‘interp’ to use spline interpolation, or ‘lsq’

for least squares shift-and-fit

OUTPUT:

offset – array. A 1D array with elements corresponding to the

vector offsets necessary to align s1 and s2 via salign.

NOTES:

• If 1D, s1 and s2 must be the same length.
• Uses cross-correlation and the nearest pixel. The disadvantage can be poorer alignment; the advantage is speed.
• Helper function for salign.