#--------------------------------------------------------------
#Author: Lisa Dang
#Created: 2016-11-09 1:21 AM EST
#Last Modified:
#Title: Photometry with Companion (PSF + Aperture)
#--------------------------------------------------------------
import numpy as np
from scipy import interpolate
from scipy import stats as st
from scipy import optimize as opt
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.patches
from astropy.io import fits
from astropy.stats import sigma_clip
from photutils import aperture_photometry
from photutils import CircularAperture, EllipticalAperture, RectangularAperture
from photutils.utils import calc_total_error
import glob
import csv
import time as tim
import os, sys
import warnings
[docs]def get_fnames(directory, AOR_snip, ch):
'''
Find paths to all the fits files.
Parameters
----------
directory : string object
Path to the directory containing all the Spitzer data.
AOR_snip : string object
Common first characters of data directory eg. 'r579'
ch : string objects
Channel used for the observation eg. 'ch1' for channel 1
Returns
-------
fname : list
List of paths to all bcd.fits files.
len(fnames): int
Number of fits file found.
'''
lst = os.listdir(directory)
AOR_list = [k for k in lst if AOR_snip in k]
#AOR_list = np.delete(AOR_list, [0, 4]) # used to ignore calibration data sets
fnames = []
for i in range(len(AOR_list)):
path = directory + '/' + AOR_list[i] + '/' + ch +'/bcd'
fnames.extend([filename for filename in glob.glob(os.path.join(path, '*bcd.fits'))])
fnames.sort()
return fnames, len(fnames)
[docs]def get_time(hdu_list, time):
'''
Gets the time stamp for each image.
Parameters
----------
hdu_list :
content of fits file.
time : 1D array
Time array to append new time stamps to.
Returns
-------
time : 1D array
Time array with new values appended.
'''
h, w, l = hdu_list[0].data.shape
t = np.linspace(hdu_list[0].header['AINTBEG'] + hdu_list[0].header['EXPTIME']/2, hdu_list[0].header['ATIMEEND'] - hdu_list[0].header['EXPTIME']/2, h)
sec2day = 1.0/(3600.0*24.0)
t = sec2day*t
time.extend(t)
return time
[docs]def sigma_clipping(image_data, filenb = 0 , fname = ['not provided'], tossed = 0, badframetable = [], bounds = (11, 19, 11, 19), **kwarg):
'''
Sigma clips bad pixels and mask entire frame if the sigma clipped
pixel is too close to the target.
Parameters
----------
image_data: 3D array
Data cube of images (2D arrays of pixel values).
filenb : int (optional)
Index of current file in the 'fname' list (list of names of files)
to keep track of the files that were tossed out. Default is 0.
fname : list (optional)
list (list of names of files) to keep track of the files that were
tossed out.
tossed : int (optional)
Total number of image tossed out. Default is 0 if none provided.
badframetable: list (optional)
List of file names and frame number of images tossed out from 'fname'.
bounds : tuple (optional)
Bounds of box around the target. Default is (11, 19 ,11, 19).
Returns
-------
sig_clipped_data: 3D array
Data cube of sigma clipped images (2D arrays of pixel values).
tossed : int
Updated total number of image tossed out.
badframetable: list
Updated list of file names and frame number of images tossed
out from 'fname'.
'''
lbx, ubx, lby, uby = bounds
h, w, l = image_data.shape
# mask invalids
image_data2 = np.ma.masked_invalid(image_data)
# make mask to mask entire bad frame
x = np.ones(shape = (w, l))
mask = np.ma.make_mask(x)
sig_clipped_data = sigma_clip(image_data2, sigma=4, iters=5, cenfunc=np.ma.median, axis = 0)
for i in range (h):
oldstar = image_data[i, lbx:ubx, lby:uby]
newstar = sig_clipped_data[i, lbx:ubx, lby:uby]
truth = newstar==oldstar
if(truth.sum() < truth.size):
sig_clipped_data[i,:,:] = np.ma.masked_array(sig_clipped_data[i,:,:], mask = mask)
badframetable.append([i,filenb,fname])
tossed += 1
return sig_clipped_data, tossed, badframetable
[docs]def bgsubtract(img_data, bg_flux = [], bg_err = [], bounds = (11, 19, 11, 19)):
'''
Measure the background level and subtracts the background from
each frame.
Parameters
----------
img_data : 3D array
Data cube of images (2D arrays of pixel values).
bg_err : 1D array (optional)
Array of uncertainties on background measurements for previous images.
Default if none given is an empty list
bounds : tuple (optional)
Bounds of box around the target to exclude from the background level
measurements. Default is (11, 19 ,11, 19).
Returns
-------
bgsub_data: 3D array
Data cube of sigma clipped images (2D arrays of pixel values).
bg_flux : 1D array
Updated array of background flux measurements for previous
images.
bg_err : 1D array
Updated array of uncertainties on background measurements for previous
images.
'''
lbx, ubx, lby, uby = bounds
image_data = np.copy(img_data)
h, w, l = image_data.shape
x = np.zeros(shape = image_data.shape)
x[:, lbx:ubx,lby:uby] = 1
mask = np.ma.make_mask(x)
masked = np.ma.masked_array(image_data, mask = mask)
masked = np.reshape(masked, (h, w*l))
bg_med = np.reshape(np.ma.median(masked, axis=1), (h, 1, 1))
bgsub_data = image_data - bg_med
bgsub_data = np.ma.masked_invalid(bgsub_data)
bg_flux.extend(bg_med.ravel())
bg_err.extend(np.ma.std(masked, axis=1))
return bgsub_data, bg_flux, bg_err
[docs]def oversampling(image_data, a = 2):
'''
First, substitutes all invalid/sigmaclipped pixel by interpolating the value.
Then oversamples the image.
Parameters
----------
image_data: 3D array
Data cube of images (2D arrays of pixel values).
a : int (optional)
Sampling factor, e.g. if a = 2, there will be twice as much data points in
the x and y axis. Default is 2. (Do not recommend larger than 2)
Returns
-------
image_over: 3D array
Data cube of oversampled images (2D arrays of pixel values).
'''
l, h, w = image_data.shape
gridx, gridy = np.mgrid[0:h:1/a, 0:w:1/a]
image_over = np.empty((l, h*a, w*a))
for i in range(l):
image_masked = np.ma.masked_invalid(image_data[i,:,:])
points = np.ma.nonzero(image_masked)
image_compre = np.ma.compressed(image_masked)
image_over[i,:,:] = interpolate.griddata(points, image_compre, (gridx, gridy), method = 'linear')
return image_over
[docs]def Gaussian2D(position, amp, xo, yo, sigx, sigy):
'''
Parameters
----------
position : 3D array
Meshgrids of x and y indices of pixels. position[:,:,0] = x and
position[:,:,1] = y.
amp : float
Amplitude of the 2D Gaussian.
xo : float
x value of the peak of the 2D Gaussian.
yo : float
y value of the peak of the 2D Gaussian.
sigx : float
Width of the 2D Gaussian along the x axis.
sigy : float
Width of the 2D Gaussian along the y axis.
Returns
-------
PSF.ravel(): 1D array
z values of the 2D Gaussian raveled.
'''
centroid = [yo, xo]
cov = [[sigx**2, 0],[0, sigy**2]]
rv = st.multivariate_normal(mean = centroid, cov = cov)
PSF = amp*(rv.pdf(position))
return PSF
[docs]def image_template(position, amp, xo, yo, sigx, sigy, amp2, xo2, yo2, sigx2, sigy2):
'''
Parameters
----------
position : 3D array
Meshgrids of x and y indices of pixels. position[:,:,0] = x and
position[:,:,1] = y.
amp : float
Amplitude of the 2D Gaussian.
xo : float
x value of the peak of the target 2D Gaussian.
yo : float
y value of the peak of the target 2D Gaussian.
sigx : float
Width of the 2D Gaussian target along the x axis.
sigy : float
Width of the 2D Gaussian target along the y axis.
amp2 : float
Amplitude of the 2D Gaussian.
xo2 : float
x value of the peak of the companion 2D Gaussian.
yo2 : float
y value of the peak of the companion 2D Gaussian.
sigx2 : float
Width of the 2D Gaussian companion along the x axis.
sigy2 : float
Width of the 2D Gaussian companion along the y axis.
Returns
-------
PSF.ravel(): 1D array
z values of the 2D Gaussian raveled.
'''
target = Gaussian2D(position, amp, xo, yo, sigx, sigy)
companion = Gaussian2D(position, amp2, xo2, yo2, sigx2, sigy2)
image = target + companion
return image.ravel()
[docs]def imagefit(image_data, tossed, popt, pcov, scale = 1,
tbounds = (9, 20, 9, 20),
pinit = (18000, 15, 15, 2, 2, 2000, 12, 13, 1, 1),
ub = (25000, 15.3, 15.3, 5, 5, 5000, 12.4, 13.6, 2, 2),
lb = (3000, 14.6, 14.6, 0.5, 0.5, 200, 11.7, 12.7, 0, 0)):
lbx, ubx, lby, uby = tbounds
lbx, ubx, lby, uby = lbx*scale, ubx*scale, lby*scale, uby*scale
l, h, w = image_data.shape
x, y = np.mgrid[lbx/scale:ubx/scale:1/scale, lby/scale:uby/scale:1/scale]
position = np.empty(x.shape + (2,))
position[:,:,0] = x
position[:,:,1] = y
popt_tmp = np.empty((l,len(pinit)))
pcov_tmp = np.empty((l,len(pinit), len(pinit)))
for i in range(l):
dataravel = image_data[i,lbx:ubx,lby:uby].ravel()
try:
popt_tmp[i,:], pcov_tmp[i,:,:] = opt.curve_fit(image_template, position, dataravel, p0=pinit, bounds = (lb,ub))
except RuntimeError:
print("Error - curve_fit failed")
popt_tmp[i,:] = np.nan
tossed += 1
popt = np.append(popt, popt_tmp, axis = 0)
pcov = np.append(pcov, pcov_tmp, axis = 0)
return popt, pcov, tossed
[docs]def companion_subtraction(image_data, c_popt):
l, h, w = image_data.shape
x, y = np.mgrid[:h:,:w:]
position = np.empty(x.shape + (2,))
position[:,:,0] = x
position[:,:,1] = y
companion = np.empty(image_data.shape)
mask = np.ones((h,w))
for i in range(l):
if (np.isnan(c_popt[i,0])==False):
companion[i,:,:] = Gaussian2D(position, c_popt[i,0], c_popt[i,1], c_popt[i,2], c_popt[i,3], c_popt[i,4])
else:
image_data[i,:,:] = np.ma.masked_array(image_data[i,:,:], mask=mask)
companion[i,:,:] = np.ma.masked_array(companion[i,:,:], mask=mask)
return image_data - companion
[docs]def centroid_FWM(image_data, xo = [], yo = [], wx = [], wy = [], scale = 1, bounds = (14, 18, 14, 18)):
'''
Gets the centroid of the target by flux weighted mean and the PSF width
of the target.
Parameters
----------
img_data : 3D array
Data cube of images (2D arrays of pixel values).
xo : list (optional)
List of x-centroid obtained previously. Default if none given is an
empty list.
yo : list (optional)
List of y-centroids obtained previously. Default if none given is an
empty list.
wx : list (optional)
List of PSF width (x-axis) obtained previously. Default if none given
is an empty list.
wy : list (optional)
List of PSF width (x-axis) obtained previously. Default if none given
is an empty list.
scale : int (optional)
If the image is oversampled, scaling factor for centroid and bounds,
i.e, give centroid in terms of the pixel value of the initial image.
bounds : tuple (optional)
Bounds of box around the target to exclude background . Default is (11, 19 ,11, 19).
Returns
-------
xo : list
Updated list of x-centroid obtained previously.
yo : list
Updated list of y-centroids obtained previously.
wx : list
Updated list of PSF width (x-axis) obtained previously.
wy : list
Updated list of PSF width (x-axis) obtained previously.
'''
lbx, ubx, lby, uby = bounds
lbx, ubx, lby, uby = lbx*scale, ubx*scale, lby*scale, uby*scale
starbox = image_data[:, lbx:ubx, lby:uby]
h, w, l = starbox.shape
# get centroid
X, Y = np.mgrid[:w,:l]
cx = (np.sum(np.sum(X*starbox, axis=1), axis=1)/(np.sum(np.sum(starbox, axis=1), axis=1))) + lbx
cy = (np.sum(np.sum(Y*starbox, axis=1), axis=1)/(np.sum(np.sum(starbox, axis=1), axis=1))) + lby
xo.extend(cx/scale)
yo.extend(cy/scale)
# get PSF widths
X, Y = np.repeat(X[np.newaxis,:,:], h, axis=0), np.repeat(Y[np.newaxis,:,:], h, axis=0)
cx, cy = np.reshape(cx, (h, 1, 1)), np.reshape(cy, (h, 1, 1))
X2, Y2 = (X + lbx - cx)**2, (Y + lby - cy)**2
widx = np.sqrt(np.sum(np.sum(X2*starbox, axis=1), axis=1)/(np.sum(np.sum(starbox, axis=1), axis=1)))
widy = np.sqrt(np.sum(np.sum(Y2*starbox, axis=1), axis=1)/(np.sum(np.sum(starbox, axis=1), axis=1)))
wx.extend(widx/scale)
wy.extend(widy/scale)
return xo, yo, wx, wy
[docs]def A_photometry(image_data, bg_err, ape_sum = [], ape_sum_err = [],
cx = 15, cy = 15, r = 2.5, a = 5, b = 5, w_r = 5, h_r = 5,
theta = 0, shape = 'Circular', method='center'):
'''
Performs aperture photometry, first by creating the aperture (Circular,
Rectangular or Elliptical), then it sums up the flux that falls into the
aperture.
Parameters
----------
image_data: 3D array
Data cube of images (2D arrays of pixel values).
bg_err : 1D array
Array of uncertainties on pixel value.
ape_sum : 1D array (optional)
Array of flux to append new flux values to. If 'None', the new values
will be appended to an empty array
ape_sum_err: 1D array (optional)
Array of flux uncertainty to append new flux uncertainty values to. If
'None', the new values will be appended to an empty array.
cx : int (optional)
x-coordinate of the center of the aperture. Default is 15.
cy : int (optional)
y-coordinate of the center of the aperture. Default is 15.
r : int (optional)
If phot_meth is 'Aperture' and ap_shape is 'Circular', c_radius is
the radius for the circular aperture. Default is 2.5.
a : int (optional)
If phot_meth is 'Aperture' and ap_shape is 'Elliptical', e_semix is
the semi-major axis for elliptical aperture (x-axis). Default is 5.
b : int (optional)
If phot_meth is 'Aperture' and ap_shape is 'Elliptical', e_semiy is
the semi-major axis for elliptical aperture (y-axis). Default is 5.
w_r : int (optional)
If phot_meth is 'Aperture' and ap_shape is 'Rectangular', r_widthx is
the full width for rectangular aperture (x-axis). Default is 5.
h_r : int (optional)
If phot_meth is 'Aperture' and ap_shape is 'Rectangular', r_widthy is
the full height for rectangular aperture (y-axis). Default is 5.
theta : int (optional)
If phot_meth is 'Aperture' and ap_shape is 'Elliptical' or
'Rectangular', theta is the angle of the rotation angle in radians
of the semimajor axis from the positive x axis. The rotation angle
increases counterclockwise. Default is 0.
shape : string object (optional)
If phot_meth is 'Aperture', ap_shape is the shape of the aperture.
Possible aperture shapes are 'Circular', 'Elliptical', 'Rectangular'.
Default is 'Circular'.
method : string object (optional)
If phot_meth is 'Aperture', apemethod is the method used to
determine the overlap of the aperture on the pixel grid. Possible
methods are 'exact', 'subpixel', 'center'. Default is 'center'.
Returns
-------
ape_sum : 1D array
Array of flux with new flux appended.
ape_sum_err: 1D array
Array of flux uncertainties with new flux uncertainties appended.
'''
l, h, w = image_data.shape
position=[cx, cy]
if (shape == 'Circular'):
aperture = CircularAperture(position, r=r)
elif (shape == 'Elliptical'):
aperture = EllipticalAperture(position, a=a, b=b, theta=theta)
elif (shape == 'Rectangular'):
aperture = RectangularAperture(position, w=w_r, h=h_r, theta=theta)
for i in range(l):
data_error = calc_total_error(image_data[i,:,:], bg_err[i], effective_gain=1)
phot_table = aperture_photometry(image_data[i,:,:],aperture, error=data_error, pixelwise_error=False)
if (phot_table['aperture_sum_err'] > 0.000001):
ape_sum.extend(phot_table['aperture_sum'])
ape_sum_err.extend(phot_table['aperture_sum_err'])
else:
ape_sum.extend([np.nan])
ape_sum_err.extend([np.nan])
return ape_sum, ape_sum_err
[docs]def binning_data(data, size):
'''
Median bin an array.
Parameters
----------
data : 1D array
Array of data to be binned.
size : int
Size of bins.
Returns
-------
binned_data: 1D array
Array of binned data.
binned_data: 1D array
Array of standard deviation for each entry in binned_data.
'''
data = np.ma.masked_invalid(data)
reshaped_data = data.reshape((len(data)/size, size))
binned_data = np.ma.median(reshaped_data, axis=1)
binned_data_std = np.std(reshaped_data, axis=1)
return binned_data, binned_data_std
[docs]def get_lightcurve(datapath, savepath, AOR_snip, channel, subarray,
save = True, save_full = '/ch2_datacube_full_AORs579.dat', bin_data = True,
bin_size = 128, save_bin = '/ch2_datacube_binned128_AORs579.dat', plot = True,
plot_name= 'CoRoT-2b.pdf', oversamp = False, savecomp = True, **kwargs):
'''
Given a directory, looks for data (bcd.fits files), opens them and performs photometry.
Parameters
----------
datapath : string object
Directory where the spitzer data is stored.
savepath : string object
Directory the outputs will be saved.
AORsnip : string objects
Common first characters of data directory eg. 'r579'
channel : string objects
Channel used for the observation eg. 'ch1' for channel 1
subarray : bool
True if observation were taken in subarray mode. False if
observation were taken in full-array mode.
save : bool (optional)
True if you want to save the outputs. Default is True.
save_full: string object (optional)
Filename of the full unbinned output data. Default is
'/ch2_datacube_full_AORs579.dat'.
bin_data : bool (optional)
True you want to get binned data. Default is True.
bin_size : int (optional)
If bin_data is True, the size of the bins. Default is 64.
save_bin : string object (optional)
Filename of the full binned output data. Default is
'/ch2_datacube_binned_AORs579.dat'.
plot : bool (optional)
True if you want to plot the time resolved lightcurve.
Default is True.
plot_name: string object (optional)
If plot and save is True, the filename of the plot to be
saved as. Default is True.
oversamp : bool (optional)
True if you want to oversample you image. Default is True.
savecomp : boold (optional)
True is you want to save companion subtracted image. Default is True.
**kwargs : dictionary
Argument passed onto other functions.
Raises
------
Error :
If Photometry method is not supported/recognized by this pipeline.
'''
print(kwargs)
# Ignore warning and starts timing
warnings.filterwarnings('ignore')
tic = tim.clock()
# get list of filenames and nb of files
fnames, nfiles = get_fnames(datapath, AOR_snip, channel)
# variables declaration
percent = 0 # to show progress while running the code
tossed = 0 # Keep tracks of number of frame discarded
badframetable = [] # list of filenames of the discarded frames
time = [] # time array
bg_flux = [] # background flux
bg_err = [] # background flux error
popt = np.empty(shape = (0,10))
pcov = np.empty(shape = (0,10,10))
xo = [] # centroid value along the x-axis
yo = [] # centroid value along the y-axis
xw = [] # PSF width along the x-axis
yw = [] # PSF width along the y-axis
aperture_sum = [] # flux obtained from aperture photometry
aperture_sum_err = [] # error on flux obtained from aperture photometry
#image_data_full=np.zeros((64*nfiles, 32, 32))
#factor = np.zeros(64*nfiles)
# data reduction & aperture photometry part
if (subarray == True):
for i in range(nfiles):
# open fits file
hdu_list = fits.open(fnames[i])
image_data0 = hdu_list[0].data
h, w, l = image_data0.shape
# get time
time = get_time(hdu_list, time)
# convert MJy/str to electron count
convfact = hdu_list[0].header['GAIN']*hdu_list[0].header['EXPTIME']/hdu_list[0].header['FLUXCONV']
image_data1 = convfact*image_data0
# sigma clip
fname = fnames[i]
image_data2, tossed, badframetable = sigma_clipping(image_data1, i ,fname[fname.find('ch2/bcd/')+8:], tossed, badframetable, **kwargs)
# bg subtract
image_data3, bg_flux, bg_err = bgsubtract(image_data2, bg_flux, bg_err)
# oversampling
if (oversamp == True):
image_data3 = np.ma.masked_invalid(oversampling(image_data3))
popt, pcov, tossed = imagefit(image_data3, tossed, popt, pcov, scale = 2)
else:
popt, pcov, tossed = imagefit(image_data3, tossed, popt, pcov)
# companion subtraction
image_data4 = companion_subtraction(image_data3, popt[-h:,5:10])
# save companion subtracted image
if (savecomp == True):
savename = fnames[i]
new_name = savepath + '/Companion_subtracted/' + savename[57:86]
image_data4.dump(new_name)
# Aperture Photometry
if (oversamp == True):
# get centroids & PSF width
xo, yo, xw, yw = centroid_FWM(image_data3, xo, yo, xw, yw, scale = 2)
# aperture photometry
aperture_sum, aperture_sum_err = A_photometry(image_data3, bg_err[-h:], aperture_sum, aperture_sum_err, cx = 2*15, cy = 2*15, r = 2*2.5, a = 2*5, b = 2*5, w_r = 2*5, h_r = 2*5, **kwargs)
else :
# get centroids & PSF width
xo, yo, xw, yw = centroid_FWM(image_data4, xo, yo, xw, yw)
# aperture photometry
aperture_sum, aperture_sum_err = A_photometry(image_data4, bg_err[-h:], aperture_sum, aperture_sum_err, **kwargs)
print('Status:', i, 'out of', nfiles)
print(np.shape(time))
print(np.shape(bg_flux))
print(np.shape(bg_err))
elif (subarray == False):
print('Sorry this part is undercontruction!')
if (bin_data == True):
binned_flux, binned_flux_std = binning_data(np.asarray(aperture_sum), bin_size)
binned_time, binned_time_std = binning_data(np.asarray(time), bin_size)
binned_xo, binned_xo_std = binning_data(np.asarray(xo), bin_size)
binned_yo, binned_yo_std = binning_data(np.asarray(yo), bin_size)
binned_xw, binned_xw_std = binning_data(np.asarray(xw), bin_size)
binned_yw, binned_yw_std = binning_data(np.asarray(yw), bin_size)
binned_bg, binned_bg_std = binning_data(np.asarray(bg_flux), bin_size)
binned_bg_err, binned_bg_err_std = binning_data(np.asarray(bg_err), bin_size)
if (plot == True):
fig, axes = plt.subplots(nrows=3, ncols=1, sharex=True, figsize=(15,5))
fig.suptitle("CoRoT-2b", fontsize="x-large")
axes[0].plot(binned_time, binned_flux,'k+')
axes[0].set_ylabel("Stellar Flux (MJy/pixel)")
axes[1].plot(binned_time, binned_xo, '+')
axes[1].set_ylabel("$x_0$")
axes[2].plot(binned_time, binned_yo, 'r+')
axes[2].set_xlabel("Time since IRAC turn-on (days)")
axes[2].set_ylabel("$y_0$")
fig.subplots_adjust(wspace=0)
if (save == True):
pathplot = savepath + '/' + plot_name
fig.savefig(pathplot)
else :
plt.show()
if (save == True):
POPT_data = popt.T
POPT_head = 'Flux, Time, x-centroid, y-centroid, x-PSF width, y-PSF width, Flux2, Time2, x-centroid2, y-centroid2, x-PSF width2, y-PSF width2'
FULL_data = np.c_[aperture_sum, aperture_sum_err, time, xo, yo, xw, yw, bg_flux, bg_err]
FULL_head = 'Flux, Flux Uncertainty, Time, x-centroid, y-centroid, x-PSF width, y-PSF width, background flux, background flux err'
BINN_data = np.c_[binned_flux, binned_flux_std, binned_time, binned_time_std, binned_xo, binned_xo_std, binned_yo, binned_yo_std, binned_xw, binned_xw_std, binned_yw, binned_yw_std, binned_bg, binned_bg_std, binned_bg_err, binned_bg_err_std]
BINN_head = 'Flux, Flux std, Time, Time std, x-centroid, x-centroid std, y-centroid, y-centroid std, x-PSF width, x-PSF width std, y-PSF width, y-PSF width std, bg, bg std, bg_err bg_err std'
pathPOPT = savepath + '/popt.dat'
pathFULL = savepath + save_full
pathBINN = savepath + save_bin
np.savetxt(pathPOPT, POPT_data, header = POPT_head)
np.savetxt(pathFULL, FULL_data, header = FULL_head)
np.savetxt(pathBINN, BINN_data, header = BINN_head)
toc = tim.clock()
print('Number of discarded frames:', tossed)
print('Time:', toc-tic, 'seconds')
if __name__=='__main__': main()