diff --git a/measureComb.py b/measureComb.py
new file mode 100755
index 0000000000000000000000000000000000000000..83c1bd8c4f19127f6ac2fe779092450245fee180
--- /dev/null
+++ b/measureComb.py
@@ -0,0 +1,239 @@
+#!/usr/bin/python
+import math, os, re
+import matplotlib as matplotlib
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+import matplotlib.cm as cm
+import numpy as np
+
+def combSpectrum(targetDirectory, flag):
+
+ # Grant David Meadors
+ # g m e a d o r s @ u m i c h . e d u
+ # 02012-11-28 (JD 2456260)
+ # combSpectrum
+ #
+ # Scan output log files for parts of the spectrum where feedforward is
+ # adversely affecting Hoft.
+ #
+ # Input argument: targetDirectory
+ # Choose a directory, e.g. "/home/pulsar/public_html/feedforward/diagnostics/LHO/H-H1_AMPS_C02_L2-9531",
+ # into which comb files have been written.
+ # Input argument: flag
+ # Flag whether the directory itself contains comb files ('one') o
+ # is a directory of directories containing comb files ('all')
+
+ # The results of the frequency comb ratio test are stored in files of the form
+ # EleutheriaComb-startGPS-stopGPS.txt
+ # and thus can be found be searching for (substitute LLO if appropriate)
+ # /home/pulsar/public_html/feedforward/diagnostics/LHO/*/*Comb*
+ # In each file, the comb frequencies are listed as the values of the columns in
+ # row 5 (using Python indexing, starting from 0), and the ratios as the values
+ # in row 6; the differences are in row 7.
+
+ if flag == 'all':
+ # The comb files are in the directories listed under the target directory.
+ highDirectoryList = os.listdir(targetDirectory)
+ highDirectoryListDirOnly = sorted([x for x in highDirectoryList if x.find('.') == -1])
+ # Because the GPS time 953100000 is used so often for testing,
+ # it is temporarily excluded:
+ highDirectoryListDirOnly = sorted([x for x in highDirectoryListDirOnly if x.find('9531') == -1])
+ combFiles = []
+ for x in highDirectoryListDirOnly:
+ targetDirectoryFiles = os.listdir(targetDirectory + x)
+ for file in targetDirectoryFiles:
+ if file.find('EleutheriaComb') > -1:
+ combFiles.append(x + '/' + file)
+ combFiles = sorted(combFiles)
+ if flag == 'one':
+ # The comb files are in the target directory.
+ # Pull a list of all the comb files.
+ targetDirectoryFiles = sorted(os.listdir(targetDirectory))
+ combFiles = []
+ for file in targetDirectoryFiles:
+ if file.find('EleutheriaComb') > -1:
+ combFiles.append(file)
+ combFiles = sorted(combFiles)
+
+ # Now write a function that can read an individual file into memory.
+ def combReader(targetDirectory, eachCombFile):
+ combFileName = targetDirectory + eachCombFile
+ try:
+ combObject = open(combFileName)
+ combLines = combObject.readlines()
+ combObject.close()
+ except IOError:
+ print 'File not found or accessible; skipping ' +\
+ combFileName
+ return combLines[5:8] + combLines[11:13]
+ # Apply the function to all the comb files in the directory.
+ for k, eachCombFile in enumerate(combFiles):
+ # Each file will have four properties:
+ # (using zero-indexed terminology)
+ # start time, read in the GPS time of the file name
+ # frequency bins, read from the fifth row of the file (0th of combRaw),
+ # ratios, read from the sixth (1st of combRaw),
+ # and difference, read from the seventh (2nd of combRaw).
+ # Then, after some comment lines,
+ # before feedforward in the eleventh line (3rd of combRaw)
+ # and after feedforward in the twelfth line (4th of combRaw)
+ if k == 0:
+ times = re.search('-(?P<GPS>\d+)-(\d+)\.', eachCombFile)
+ timeArray = np.asarray(times.group(1), dtype=np.int32)
+ combRaw = combReader(targetDirectory, eachCombFile)
+ frequencyArray = np.asarray(str(combRaw[0]).split(), dtype=np.float32)
+ frequencyArray = frequencyArray.astype(int)
+ ratioArray = np.asarray(str(combRaw[1]).split(), dtype=np.float32)
+ differenceArray = np.asarray(str(combRaw[2]).split(), dtype=np.float32)
+ beforeArray = np.asarray(str(combRaw[3]).split(), dtype=np.float32)
+ afterArray = np.asarray(str(combRaw[4]).split(), dtype=np.float32)
+ if k > 0:
+ times = re.search('-(?P<GPS>\d+)-(\d+)\.', eachCombFile)
+ timeArray = np.vstack([timeArray, np.asarray(times.group(1), dtype=np.int32)])
+ combRaw = combReader(targetDirectory, eachCombFile)
+ frequencyArray = np.vstack([frequencyArray, np.asarray(str(combRaw[0]).split(), dtype=np.float32)])
+ ratioArray = np.vstack([ratioArray, np.asarray(str(combRaw[1]).split(), dtype=np.float32)])
+ differenceArray = np.vstack([differenceArray, np.asarray(str(combRaw[2]).split(), dtype=np.float32)])
+ beforeArray = np.vstack([beforeArray, np.asarray(str(combRaw[3]).split(), dtype=np.float32)])
+ afterArray = np.vstack([afterArray, np.asarray(str(combRaw[4]).split(), dtype=np.float32)])
+
+ # Now we are going to plot these arrays.
+ # First, set a directory for the output. At least for the time being,
+ # that can be the same directory as the target directory.
+ # Note the inelegant way from extracting the start time as a label.
+ graphTitleRatio = targetDirectory + "EleutheriaPostPlotCombRatio" +\
+ '-' + str(timeArray[0]).strip("[").strip("]").strip("'")
+ graphTitleDiff = targetDirectory + "EleutheriaPostPlotCombDiff" +\
+ '-' + str(timeArray[0]).strip("[").strip("]").strip("'")
+ graphTitleBeforeAfter = targetDirectory + "EleutheriaPostPlotCombBeforeAfter" +\
+ '-' + str(timeArray[0]).strip("[").strip("]").strip("'")
+ graphWholeRun = "EleutheriaWholeRun"
+ x, y = np.meshgrid(timeArray, frequencyArray[0])
+ extensions = [x[0, 0], x[-1, -1], y[0, 0], y[-1, -1]]
+ # Define a helpfully broken color-scheme as found on the
+ # SciPy Matplotlib Cookbook,
+ # http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps
+ #cdict = {'red': ((0.0, 0.0, 0.0),\
+ #(0.5, 1.0, 0.7),\
+ #(1.0, 1.0, 1.0)),\
+ #'green': ((0.0, 0.0, 0.0),\
+ #(0.5, 1.0, 0.0),\
+ #(1.0, 1.0, 1.0)),\
+ #'blue': ((0.0, 0.0, 0.0),\
+ #(0.5, 1.0, 0.0),\
+ #(1.0, 0.5, 1.0))}
+ #my_cmap = matplotlib.colors.LinearSegmentedColormap('my_colormap', cdict, 256)
+ #plt.figure()
+ #CSratio = plt.contourf(x.T, y.T, ratioArray, \
+ #np.asarray([0.8, 0.85, 0.9, 0.95, \
+ #1, 1.05, 1.1, 1.15, 1.2]), cmap=my_cmap)
+ #plt.colorbar(CSratio, shrink=0.8, extend='both')
+ #plt.xlabel('GPS time (s)')
+ #plt.ylabel('Frequency (Hz)')
+ #plt.title('Post/pre-filtering Hoft ratio (lower is better)')
+ #plt.savefig(graphTitleRatio + 'Contour.png')
+ #plt.savefig(graphTitleRatio + 'Contour.pdf')
+ #plt.close()
+ #plt.clf()
+ #fig = plt.figure()
+ #ax = fig.add_subplot(111)
+ #CSratioIm = ax.imshow(ratioArray.T, origin='lower', \
+ #interpolation = 'nearest', extent=extensions, vmin=0.8, vmax=1.2,\
+ #cmap=my_cmap)
+ #CSratioIm = fig.colorbar(CSratioIm, shrink = 0.8, extend = 'both')
+ #ax.set_aspect('auto')
+ #ax.set_xlabel('GPS time (s)')
+ #ax.set_ylabel('Frequency (Hz)')
+ #ax.set_title('Post/pre-filtering Hoft ratio (lower is better)')
+ #fig.savefig(graphTitleRatio + '.png')
+ #fig.savefig(graphTitleRatio + '.pdf')
+ #plt.clf()
+ #fig = plt.figure()
+ #ax = fig.add_subplot(111)
+ #CSdiffIm = ax.imshow(-differenceArray.T, origin='lower', \
+ #interpolation = 'nearest' , extent=extensions, vmin=-2e-24, vmax=2e-24, \
+ #cmap=my_cmap)
+ #plt.show()
+ #CSdiffIm = fig.colorbar(CSdiffIm, shrink = 0.8, extend = 'both')
+ #ax.set_aspect('auto')
+ #ax.set_xlabel('GPS time (s)')
+ #ax.set_ylabel('Frequency (Hz)')
+ #ax.set_title('Post - pre Hoft difference (lower is better)')
+ #fig.savefig(graphTitleDiff + '.png')
+ #fig.savefig(graphTitleDiff + '.pdf')
+ #plt.clf()
+ #plt.figure()
+ #CSdiff = plt.contourf(x.T, y.T, -differenceArray, \
+ #np.asarray([-2e-24, -1.5e-24, -1e-24, -5e-25,\
+ #0, 5e-25, 1e-24, 1.5e-24, 2e-24]), cmap=my_cmap)
+ #plt.colorbar(CSdiff, shrink=0.8, extend='both')
+ #plt.xlabel('GPS time (s)')
+ #plt.ylabel('Frequency (Hz)')
+ #plt.title('Post - pre Hoft difference (lower is better)')
+ #plt.savefig(graphTitleDiff + 'Contour.png')
+ #plt.savefig(graphTitleDiff + 'Contour.pdf')
+ #plt.close()
+ #plt.figure()
+ #freqList = [25, 27, 29]
+ #p0 = plt.scatter(timeArray, -differenceArray[:, 25], color='b')
+ #mean0 = np.mean(-differenceArray[:, 25])
+ #p1 = plt.scatter(timeArray, -differenceArray[:, 27], color='r')
+ #mean1 = np.mean(-differenceArray[:, 27])
+ #p2 = plt.scatter(timeArray, -differenceArray[:, 29], color='k')
+ #mean2 = np.mean(-differenceArray[:, 29])
+ #plt.grid(True)
+ #plt.xlabel('GPS time (s)')
+ #plt.ylabel('Post - pre Hoft difference (lower is better)')
+ #plt.title('Difference vs time by frequency')
+ #plt.legend([p0, p1, p2], \
+ #[str(frequencyArray[0, freqList[0]]) + ' Hz, mean: ' + str(mean0),\
+ #str(frequencyArray[0, freqList[1]]) + ' Hz, mean: ' + str(mean1),\
+ #str(frequencyArray[0, freqList[2]]) + ' Hz, mean: ' + str(mean2)])
+ #plt.savefig(graphTitleDiff + 'Freq.png')
+ #plt.savefig(graphTitleDiff + 'Freq.pdf')
+ #plt.close()
+ #plt.figure()
+ #p0 = plt.scatter(timeArray, beforeArray[:, 27], color='b')
+ #mean0 = np.mean(beforeArray[:, 27])
+ #numberofthings = np.shape(beforeArray[:,27])[0]
+ #mean0array = mean0 * np.ones((np.shape(beforeArray[:, 27])[0], 1))
+ #p1 = plt.scatter(timeArray, afterArray[:, 27], color='r')
+ #mean1 = np.mean(afterArray[:, 27])
+ #mean1array = mean1 * np.ones((np.shape(afterArray[:, 27])[0], 1))
+ #p2 = plt.plot(timeArray, mean0array, 'b^')
+ #p3 = plt.plot(timeArray, mean1array, 'r+')
+ #plt.grid(True)
+ #plt.xlabel('GPS time (s)')
+ #plt.ylabel('Hoft')
+ #plt.title('Before and after feedforward')
+ #plt.legend([p0, p1, p2, p3], \
+ #[str(frequencyArray[0, freqList[1]]) + ' Hz before, scatterplot', \
+ #str(frequencyArray[0, freqList[1]]) + ' Hz before, mean: ' + str(mean0), \
+ #str(frequencyArray[0, freqList[1]]) + ' Hz after, scatterplot', \
+ #str(frequencyArray[0, freqList[1]]) + ' Hz after, mean: ' + str(mean1), \
+ #])
+ #plt.savefig(graphTitleBeforeAfter + '.png')
+ #plt.savefig(graphTitleBeforeAfter + '.pdf')
+ #plt.legend
+ #plt.close()
+ plt.figure()
+ p0 = plt.scatter(timeArray, beforeArray[:, 27], color='b')
+ plt.grid(True)
+ plt.savefig(graphWholeRun + '.png')
+ plt.savefig(graphWholeRun + '.pdf')
+ plt.close()
+
+# Uncomment below to test on one directory only:
+combSpectrum('/home/pulsar/public_html/feedforward/diagnostics/LHO/H-H1_AMPS_C02_L2-9326/', 'one')
+#combSpectrum('/home/pulsar/feedforward/2012/10/25/AMPS/files_for_comb/', 'one')
+# Uncomment below to test all directories and produce a whole-run overview.
+#grandTarget = '/home/pulsar/public_html/feedforward/diagnostics/LHO/'
+#combSpectrum(grandTarget, 'all')
+# Uncomment below to test each directory, making plots one-by-one.
+#grandTarget = '/home/pulsar/public_html/feedforward/diagnostics/LHO'
+#highDirectoryList = os.listdir(grandTarget)
+#highDirectoryListDirOnly = [x for x in highDirectoryList if x.find('.') == -1]
+#for x in highDirectoryListDirOnly:
+# combSpectrum(grandTarget + '/' + x + '/', 'one')
+
+