From 6855e4803ff8873e08ccc6903294dd422e333a00 Mon Sep 17 00:00:00 2001
From: Aaron Viets <aaron.viets@ligo.org>
Date: Tue, 16 Feb 2021 06:33:13 -0800
Subject: [PATCH] gstlal-calibration:  plotting script updates

---
 .../tests/check_calibration/Makefile          |   8 +-
 .../check_calibration/compute_fs_over_Q.py    |  44 +++
 .../tests/check_calibration/compute_tau.py    | 110 ++++++++
 .../check_calibration/plot_kappas_from_txt.py | 265 ++++++++++++++++++
 4 files changed, 426 insertions(+), 1 deletion(-)
 create mode 100644 gstlal-calibration/tests/check_calibration/compute_fs_over_Q.py
 create mode 100644 gstlal-calibration/tests/check_calibration/compute_tau.py
 create mode 100644 gstlal-calibration/tests/check_calibration/plot_kappas_from_txt.py

diff --git a/gstlal-calibration/tests/check_calibration/Makefile b/gstlal-calibration/tests/check_calibration/Makefile
index 6018213f60..818c4e4c54 100644
--- a/gstlal-calibration/tests/check_calibration/Makefile
+++ b/gstlal-calibration/tests/check_calibration/Makefile
@@ -253,10 +253,16 @@ kappastimeseries_GDS: $(IFO)1_hoft_GDS_frames.cache $(IFO)1_easy_raw_frames.cach
 
 exactkappastimeseries: $(IFO)1_hoft_DCS_APPROXKAPPAS_frames.cache $(IFO)1_hoft_DCS_EXACTKAPPAS_frames.cache
 	python3 frame_manipulator.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --frame-cache $(IFO)1_hoft_DCS_EXACTKAPPAS_frames.cache --output-path txt --channel-list 'DCS-CALIB_KAPPA_TST_REALEXACTKAPPAS,DCS-CALIB_KAPPA_TST_IMAGINARYEXACTKAPPAS,DCS-CALIB_KAPPA_PUM_REALEXACTKAPPAS,DCS-CALIB_KAPPA_PUM_IMAGINARYEXACTKAPPAS,DCS-CALIB_KAPPA_UIM_REALEXACTKAPPAS,DCS-CALIB_KAPPA_UIM_IMAGINARYEXACTKAPPAS,DCS-CALIB_KAPPA_CEXACTKAPPAS,DCS-CALIB_F_CCEXACTKAPPAS,DCS-CALIB_F_S_SQUAREDEXACTKAPPAS,DCS-CALIB_SRC_Q_INVERSEEXACTKAPPAS'
+	python3 compute_tau.py --kappar-txt $(IFO)1-DCS-CALIB_KAPPA_TST_REALEXACTKAPPAS.txt --kappai-txt $(IFO)1-DCS-CALIB_KAPPA_TST_IMAGINARYEXACTKAPPAS.txt --config-file $(DCSEXACTKAPPASCONFIGS) --kappa-filename $(IFO)1-DCS-CALIB_KAPPA_TST_EXACTKAPPAS.txt --tau-filename $(IFO)1-DCS-CALIB_TAU_TST_EXACTKAPPAS.txt
+	python3 compute_tau.py --kappar-txt $(IFO)1-DCS-CALIB_KAPPA_PUM_REALEXACTKAPPAS.txt --kappai-txt $(IFO)1-DCS-CALIB_KAPPA_PUM_IMAGINARYEXACTKAPPAS.txt --config-file $(DCSEXACTKAPPASCONFIGS) --kappa-filename $(IFO)1-DCS-CALIB_KAPPA_PUM_EXACTKAPPAS.txt --tau-filename $(IFO)1-DCS-CALIB_TAU_PUM_EXACTKAPPAS.txt
+	python3 compute_tau.py --kappar-txt $(IFO)1-DCS-CALIB_KAPPA_UIM_REALEXACTKAPPAS.txt --kappai-txt $(IFO)1-DCS-CALIB_KAPPA_UIM_IMAGINARYEXACTKAPPAS.txt --config-file $(DCSEXACTKAPPASCONFIGS) --kappa-filename $(IFO)1-DCS-CALIB_KAPPA_UIM_EXACTKAPPAS.txt --tau-filename $(IFO)1-DCS-CALIB_TAU_UIM_EXACTKAPPAS.txt
 	python3 compute_fs_over_Q.py --fs-squared-txt $(IFO)1-DCS-CALIB_F_S_SQUAREDEXACTKAPPAS.txt --Qinv-txt $(IFO)1-DCS-CALIB_SRC_Q_INVERSEEXACTKAPPAS.txt --filename $(IFO)1-DCS-CALIB_F_S_OVER_QEXACTKAPPAS.txt
 	python3 frame_manipulator.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --frame-cache $(IFO)1_hoft_DCS_APPROXKAPPAS_frames.cache --output-path txt --channel-list 'DCS-CALIB_KAPPA_TST_REALAPPROXKAPPAS,DCS-CALIB_KAPPA_TST_IMAGINARYAPPROXKAPPAS,DCS-CALIB_KAPPA_PUM_REALAPPROXKAPPAS,DCS-CALIB_KAPPA_PUM_IMAGINARYAPPROXKAPPAS,DCS-CALIB_KAPPA_UIM_REALAPPROXKAPPAS,DCS-CALIB_KAPPA_UIM_IMAGINARYAPPROXKAPPAS,DCS-CALIB_KAPPA_CAPPROXKAPPAS,DCS-CALIB_F_CCAPPROXKAPPAS,DCS-CALIB_F_S_SQUAREDAPPROXKAPPAS,DCS-CALIB_SRC_Q_INVERSEAPPROXKAPPAS'
+	python3 compute_tau.py --kappar-txt $(IFO)1-DCS-CALIB_KAPPA_TST_REALAPPROXKAPPAS.txt --kappai-txt $(IFO)1-DCS-CALIB_KAPPA_TST_IMAGINARYAPPROXKAPPAS.txt --config-file $(DCSAPPROXKAPPASCONFIGS) --kappa-filename $(IFO)1-DCS-CALIB_KAPPA_TST_APPROXKAPPAS.txt --tau-filename $(IFO)1-DCS-CALIB_TAU_TST_APPROXKAPPAS.txt
+	python3 compute_tau.py --kappar-txt $(IFO)1-DCS-CALIB_KAPPA_PUM_REALAPPROXKAPPAS.txt --kappai-txt $(IFO)1-DCS-CALIB_KAPPA_PUM_IMAGINARYAPPROXKAPPAS.txt --config-file $(DCSAPPROXKAPPASCONFIGS) --kappa-filename $(IFO)1-DCS-CALIB_KAPPA_PUM_APPROXKAPPAS.txt --tau-filename $(IFO)1-DCS-CALIB_TAU_PUM_APPROXKAPPAS.txt
+	python3 compute_tau.py --kappar-txt $(IFO)1-DCS-CALIB_KAPPA_UIM_REALAPPROXKAPPAS.txt --kappai-txt $(IFO)1-DCS-CALIB_KAPPA_UIM_IMAGINARYAPPROXKAPPAS.txt --config-file $(DCSAPPROXKAPPASCONFIGS) --kappa-filename $(IFO)1-DCS-CALIB_KAPPA_UIM_APPROXKAPPAS.txt --tau-filename $(IFO)1-DCS-CALIB_TAU_UIM_APPROXKAPPAS.txt
 	python3 compute_fs_over_Q.py --fs-squared-txt $(IFO)1-DCS-CALIB_F_S_SQUAREDAPPROXKAPPAS.txt --Qinv-txt $(IFO)1-DCS-CALIB_SRC_Q_INVERSEAPPROXKAPPAS.txt --filename $(IFO)1-DCS-CALIB_F_S_OVER_QAPPROXKAPPAS.txt
-	python3 plot_kappas_from_txt.py --txt-list '$(IFO)1-DCS-CALIB_KAPPA_TST_REALAPPROXKAPPAS.txt,$(IFO)1-DCS-CALIB_KAPPA_PUM_REALAPPROXKAPPAS.txt,$(IFO)1-DCS-CALIB_KAPPA_UIM_REALAPPROXKAPPAS.txt;$(IFO)1-DCS-CALIB_KAPPA_TST_REALEXACTKAPPAS.txt,$(IFO)1-DCS-CALIB_KAPPA_PUM_REALEXACTKAPPAS.txt,$(IFO)1-DCS-CALIB_KAPPA_UIM_REALEXACTKAPPAS.txt:$(IFO)1-DCS-CALIB_KAPPA_TST_IMAGINARYAPPROXKAPPAS.txt,$(IFO)1-DCS-CALIB_KAPPA_PUM_IMAGINARYAPPROXKAPPAS.txt,$(IFO)1-DCS-CALIB_KAPPA_UIM_IMAGINARYAPPROXKAPPAS.txt;$(IFO)1-DCS-CALIB_KAPPA_TST_IMAGINARYEXACTKAPPAS.txt,$(IFO)1-DCS-CALIB_KAPPA_PUM_IMAGINARYEXACTKAPPAS.txt,$(IFO)1-DCS-CALIB_KAPPA_UIM_IMAGINARYEXACTKAPPAS.txt' --labels 'Approx;Exact' --filename actuation_TDCFs.png
+	python3 plot_kappas_from_txt.py --txt-list '$(IFO)1-DCS-CALIB_KAPPA_TST_APPROXKAPPAS.txt,$(IFO)1-DCS-CALIB_KAPPA_PUM_APPROXKAPPAS.txt,$(IFO)1-DCS-CALIB_KAPPA_UIM_APPROXKAPPAS.txt;$(IFO)1-DCS-CALIB_KAPPA_TST_EXACTKAPPAS.txt,$(IFO)1-DCS-CALIB_KAPPA_PUM_EXACTKAPPAS.txt,$(IFO)1-DCS-CALIB_KAPPA_UIM_EXACTKAPPAS.txt:$(IFO)1-DCS-CALIB_TAU_TST_APPROXKAPPAS.txt,$(IFO)1-DCS-CALIB_TAU_PUM_APPROXKAPPAS.txt,$(IFO)1-DCS-CALIB_TAU_UIM_APPROXKAPPAS.txt;$(IFO)1-DCS-CALIB_TAU_TST_EXACTKAPPAS.txt,$(IFO)1-DCS-CALIB_TAU_PUM_EXACTKAPPAS.txt,$(IFO)1-DCS-CALIB_TAU_UIM_EXACTKAPPAS.txt' --labels 'Approx;Exact' --filename actuation_TDCFs.png
 	python3 plot_kappas_from_txt.py --txt-list '$(IFO)1-DCS-CALIB_KAPPA_CAPPROXKAPPAS.txt;$(IFO)1-DCS-CALIB_KAPPA_CEXACTKAPPAS.txt:$(IFO)1-DCS-CALIB_F_CCAPPROXKAPPAS.txt;$(IFO)1-DCS-CALIB_F_CCEXACTKAPPAS.txt' --labels 'Approx;Exact' --filename sensing_TDCFs.png
 	python3 plot_kappas_from_txt.py --txt-list '$(IFO)1-DCS-CALIB_F_S_SQUAREDAPPROXKAPPAS.txt;$(IFO)1-DCS-CALIB_F_S_SQUAREDEXACTKAPPAS.txt:$(IFO)1-DCS-CALIB_F_S_OVER_QAPPROXKAPPAS.txt;$(IFO)1-DCS-CALIB_F_S_OVER_QEXACTKAPPAS.txt' --labels 'Approx;Exact' --filename SRC_TDCFs.png
 
diff --git a/gstlal-calibration/tests/check_calibration/compute_fs_over_Q.py b/gstlal-calibration/tests/check_calibration/compute_fs_over_Q.py
new file mode 100644
index 0000000000..d3dd4ab99c
--- /dev/null
+++ b/gstlal-calibration/tests/check_calibration/compute_fs_over_Q.py
@@ -0,0 +1,44 @@
+#!/usr/bin/env python3
+# Copyright (C) 2021  Aaron Viets
+#
+# This program is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 2 of the License, or (at your
+# option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
+# Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with this program; if not, write to the Free Software Foundation, Inc.,
+# 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+
+
+import numpy as np
+from optparse import OptionParser, Option
+
+
+parser = OptionParser()
+parser.add_option("--fs-squared-txt", metavar = "name", help = "Name of txt file with a time series of fs^2 values")
+parser.add_option("--Qinv-txt", metavar = "name", help = "Name of txt file with a time series of 1/Q values")
+parser.add_option("--filename", metavar = "name", type = str, default = "fs_over_Q.txt", help = "Name of output file")
+
+options, filenames = parser.parse_args()
+
+
+fs2 = np.loadtxt(options.fs_squared_txt)
+Qinv = np.loadtxt(options.Qinv_txt)
+t = np.transpose(fs2)[0]
+fs2 = np.transpose(fs2)[1]
+Qinv = np.transpose(Qinv)[1]
+
+fs_over_Q = np.zeros(len(fs2))
+
+for i in range(len(fs2)):
+	fs_over_Q[i] = np.sign(fs2[i]) * np.sqrt(abs(fs2[i])) * Qinv[i]
+
+np.savetxt(options.filename, np.transpose([t, fs_over_Q]), fmt = ['%f', '%8e'], delimiter = "   ")
+
+
diff --git a/gstlal-calibration/tests/check_calibration/compute_tau.py b/gstlal-calibration/tests/check_calibration/compute_tau.py
new file mode 100644
index 0000000000..8ef3a56c2a
--- /dev/null
+++ b/gstlal-calibration/tests/check_calibration/compute_tau.py
@@ -0,0 +1,110 @@
+#!/usr/bin/env python3
+# Copyright (C) 2021  Aaron Viets
+#
+# This program is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 2 of the License, or (at your
+# option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
+# Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with this program; if not, write to the Free Software Foundation, Inc.,
+# 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+
+
+import numpy as np
+import os
+from optparse import OptionParser, Option
+import configparser
+
+
+parser = OptionParser()
+parser.add_option("--kappar-txt", metavar = "name", help = "Name of txt file with a time series of real parts of kappa")
+parser.add_option("--kappai-txt", metavar = "name", help = "Name of txt file with a time series of imaginary parts of kappa")
+parser.add_option("--config-file", metavar = "name", help = "Name of config file used to produce the TDCF data")
+parser.add_option("--kappa-filename", metavar = "name", type = str, help = "Name of output file with magnitude of TDCF")
+parser.add_option("--tau-filename", metavar = "name", type = str, help = "Name of output file with tau, the time advance associated with the TDCF")
+
+options, filenames = parser.parse_args()
+
+kr = np.loadtxt(options.kappar_txt)
+ki = np.loadtxt(options.kappai_txt)
+t = np.transpose(kr)[0]
+kr = np.transpose(kr)[1]
+ki = np.transpose(ki)[1]
+
+kappa = np.zeros(len(t))
+tau = np.zeros(len(t))
+
+
+# Find the frequency of the actuator line used to measure this TDCF.
+# First, we need the config file.
+def ConfigSectionMap(section):
+	dict1 = {}
+	options = Config.options(section)
+	for option in options:
+		try:
+			dict1[option] = Config.get(section, option)
+			if dict1[option] == -1:
+				DebugPrint("skip: %s" % option)
+		except:
+			print("exception on %s!" % option)
+			dict[option] = None
+	return dict1
+
+Config = configparser.ConfigParser()
+Config.read(options.config_file)
+InputConfigs = ConfigSectionMap("InputConfigurations")
+
+# Next, we need the filters file.
+# Search the directory tree for files with names matching the one we want.
+filters_name = InputConfigs["filtersfilename"]
+if filters_name.count('/') > 0:
+	# Then the path to the filters file was given
+	filters = np.load(filters_name)
+else:
+	# We need to search for the filters file
+	filters_paths = []
+	print("\nSearching for %s ..." % filters_name)
+	# Check the user's home directory
+	for dirpath, dirs, files in os.walk(os.environ['HOME']):
+		if filters_name in files:
+			# We prefer filters that came directly from a GDSFilters directory of the calibration SVN
+			if dirpath.count("GDSFilters") > 0:
+				filters_paths.insert(0, os.path.join(dirpath, filters_name))
+			else:
+				filters_paths.append(os.path.join(dirpath, filters_name))
+	# Check if there is a checkout of the entire calibration SVN
+	for dirpath, dirs, files in os.walk('/ligo/svncommon/CalSVN/aligocalibration/trunk/Runs/'):
+		if filters_name in files:
+			# We prefer filters that came directly from a GDSFilters directory of the calibration SVN
+			if dirpath.count("GDSFilters") > 0:
+				filters_paths.insert(0, os.path.join(dirpath, filters_name))
+			else:
+				filters_paths.append(os.path.join(dirpath, filters_name))
+	if not len(filters_paths):
+		raise ValueError("Cannot find filters file %s in home directory %s or in /ligo/svncommon/CalSVN/aligocalibration/trunk/Runs/*/GDSFilters", (filters_name, os.environ['HOME']))
+	print("Loading calibration filters from %s\n" % filters_paths[0])
+	filters = np.load(filters_paths[0])
+
+# Next, we need to infer what TDCF we are dealing with to choose the line.
+if 'TST' in options.kappar_txt or 'tst' in options.kappar_txt:
+	freq = filters['ktst_esd_line_freq']
+elif 'PUM' in options.kappar_txt or 'pum' in options.kappar_txt:
+	freq = filters['pum_act_line_freq']
+else:
+	freq = filters['uim_act_line_freq']
+
+
+for i in range(len(t)):
+	kappa[i] = abs(kr[i] + 1j * ki[i])
+	tau[i] = np.angle(kr[i] + 1j * ki[i]) / 2 / np.pi / freq * 1000000.0
+
+np.savetxt(options.kappa_filename, np.transpose([t, kappa]), fmt = ['%f', '%8e'], delimiter = "   ")
+np.savetxt(options.tau_filename, np.transpose([t, tau]), fmt = ['%f', '%8e'], delimiter = "   ")
+
+
diff --git a/gstlal-calibration/tests/check_calibration/plot_kappas_from_txt.py b/gstlal-calibration/tests/check_calibration/plot_kappas_from_txt.py
new file mode 100644
index 0000000000..d6555a1b61
--- /dev/null
+++ b/gstlal-calibration/tests/check_calibration/plot_kappas_from_txt.py
@@ -0,0 +1,265 @@
+#!/usr/bin/env python3
+# Copyright (C) 2020  Aaron Viets
+#
+# This program is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 2 of the License, or (at your
+# option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
+# Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with this program; if not, write to the Free Software Foundation, Inc.,
+# 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+
+#
+# =============================================================================
+#
+#				   Preamble
+#
+# =============================================================================
+#
+
+
+import matplotlib; matplotlib.use('Agg')
+import sys
+import os
+import numpy as np
+import time
+from math import pi
+import resource
+import datetime
+from matplotlib import rc
+rc('text', usetex = True)
+matplotlib.rcParams['font.family'] = 'Times New Roman'
+matplotlib.rcParams['font.size'] = 32
+matplotlib.rcParams['legend.fontsize'] = 20
+matplotlib.rcParams['mathtext.default'] = 'regular'
+import glob
+import matplotlib.pyplot as plt
+from ticks_and_grid import ticks_and_grid
+
+from optparse import OptionParser, Option
+
+from lal import LIGOTimeGPS
+
+parser = OptionParser()
+parser.add_option("--txt-list", metavar = "list", help = "colon-, semicolon-, and comma-separated list of txt files to read in.  Commas separate different datasets from the same calibration type to be placed on the same plot.  Semicolons separate different calibration types to be placed on the same plot.  Colons separate datasets to be placed on different plots arranged vertically on the same figure.")
+parser.add_option("--labels", metavar = "list", help = "Semicolon-separated list of names to use in plot legends for different calibrated data sets")
+parser.add_option("--plot-title", metavar = "name", default = None, help = "Title to appear at the top of the figure")
+parser.add_option("--filename", metavar = "name", type = str, default = "", help = "Name of output file.  GPS times will be added at the end of the filename.")
+
+options, filenames = parser.parse_args()
+
+# Set up list of labels to be used in plot legends and filenames
+short_labels = options.labels.split(';')
+
+# Arrange list of data sets
+t_start = np.inf
+t_end = 0
+data = []
+labels = []
+ylabels = []
+data_min = []
+data_max = []
+txt_list = options.txt_list.split(':')
+for i in range(len(txt_list)):
+	data.append([])
+	labels.append([])
+	ylabels.append([])
+	txt_list[i] = txt_list[i].split(';')
+	for j in range(len(txt_list[i])):
+		data[i].append([])
+		labels[i].append([])
+		ylabels[i].append(None)
+		txt_list[i][j] = txt_list[i][j].split(',')
+		for k in range(len(txt_list[i][j])):
+			data[i][j].append(np.loadtxt(txt_list[i][j][k]).transpose())
+			labels[i][j].append(short_labels[j])
+			t_start = min(t_start, data[i][j][k][0][0])
+			t_end = max(t_end, data[i][j][k][0][-1])
+			# Find out which TDCF this is and add that to the label
+			if 'tst' in txt_list[i][j][k] or 'TST' in txt_list[i][j][k]:
+				if 'imag' in txt_list[i][j][k] or 'IMAG' in txt_list[i][j][k]:
+					labels[i][j][k] += r': $\Im(\kappa_{\rm T})$'
+					ylabels[i][j] = "Correction"
+					if j == 0 and k == 0:
+						data_min.append(-0.1)
+						data_max.append(0.1)
+				elif 'real' in txt_list[i][j][k] or 'REAL' in txt_list[i][j][k]:
+					labels[i][j][k] += r': $\Re(\kappa_{\rm T})$'
+					ylabels[i][j] = "Correction"
+					if j == 0 and k == 0:
+						data_min.append(0.9)
+						data_max.append(1.1)
+				elif 'tau' in txt_list[i][j][k] or 'TAU' in txt_list[i][j][k]:
+					labels[i][j][k] += r': $\tau_{\rm T}$'
+					ylabels[i][j] = "Time ($\mu$s)"
+					if j == 0 and k == 0:
+						data_min.append(-1000)
+						data_max.append(1000)
+				else:
+					# Assume it's the magnitude
+					labels[i][j][k] += r': $\kappa_{\rm T}$'
+					ylabels[i][j] = "Correction"
+					if j == 0 and k == 0:
+						data_min.append(0.9)
+						data_max.append(1.1)
+			elif 'pum' in txt_list[i][j][k] or 'PUM' in txt_list[i][j][k]:
+				if 'imag' in txt_list[i][j][k] or 'IMAG' in txt_list[i][j][k]:
+					labels[i][j][k] += r': $\Im(\kappa_{\rm P})$'
+					ylabels[i][j] = "Correction"
+					if j == 0 and k == 0:
+						data_min.append(-0.1)
+						data_max.append(0.1)
+				elif 'real' in txt_list[i][j][k] or 'REAL' in txt_list[i][j][k]:
+					labels[i][j][k] += r': $\Re(\kappa_{\rm P})$'
+					ylabels[i][j] = "Correction"
+					if j == 0 and k == 0:
+						data_min.append(0.9)
+						data_max.append(1.1)
+				elif 'tau' in txt_list[i][j][k] or 'TAU' in txt_list[i][j][k]:
+					labels[i][j][k] += r': $\tau_{\rm P}$'
+					ylabels[i][j] = "Time ($\mu$s)"
+					if j == 0 and k == 0:
+						data_min.append(-1000)
+						data_max.append(1000)
+				else:
+					# Assume it's the magnitude
+					labels[i][j][k] += r': $\kappa_{\rm P}$'
+					ylabels[i][j] = "Correction"
+					if j == 0 and k == 0:
+						data_min.append(0.9)
+						data_max.append(1.1)
+			elif 'uim' in txt_list[i][j][k] or 'UIM' in txt_list[i][j][k]:
+				if 'imag' in txt_list[i][j][k] or 'IMAG' in txt_list[i][j][k]:
+					labels[i][j][k] += r': $\Im(\kappa_{\rm U})$'
+					ylabels[i][j] = "Correction"
+					if j == 0 and k == 0:
+						data_min.append(-0.1)
+						data_max.append(0.1)
+				elif 'real' in txt_list[i][j][k] or 'REAL' in txt_list[i][j][k]:
+					labels[i][j][k] += r': $\Re(\kappa_{\rm U})$'
+					ylabels[i][j] = "Correction"
+					if j == 0 and k == 0:
+						data_min.append(0.9)
+						data_max.append(1.1)
+				elif 'tau' in txt_list[i][j][k] or 'TAU' in txt_list[i][j][k]:
+					labels[i][j][k] += r': $\tau_{\rm U}$'
+					ylabels[i][j] = "Time ($\mu$s)"
+					if j == 0 and k == 0:
+						data_min.append(-1000)
+						data_max.append(1000)
+				else:
+					# Assume it's the magnitude
+					labels[i][j][k] += r': $\kappa_{\rm U}$'
+					ylabels[i][j] = "Correction"
+					if j == 0 and k == 0:
+						data_min.append(0.9)
+						data_max.append(1.1)
+			elif 'pu' in txt_list[i][j][k] or 'PU' in txt_list[i][j][k]:
+				if 'imag' in txt_list[i][j][k] or 'IMAG' in txt_list[i][j][k]:
+					labels[i][j][k] += r': $\Im(\kappa_{\rm PU})$'
+					ylabels[i][j] = "Correction"
+					if j == 0 and k == 0:
+						data_min.append(-0.1)
+						data_max.append(0.1)
+				elif 'real' in txt_list[i][j][k] or 'REAL' in txt_list[i][j][k]:
+					labels[i][j][k] += r': $\Re(\kappa_{\rm PU})$'
+					ylabels[i][j] = "Correction"
+					if j == 0 and k == 0:
+						data_min.append(0.9)
+						data_max.append(1.1)
+				elif 'tau' in txt_list[i][j][k] or 'TAU' in txt_list[i][j][k]:
+					labels[i][j][k] += r': $\tau_{\rm PU}$'
+					ylabels[i][j] = "Time ($\mu$s)"
+					if j == 0 and k == 0:
+						data_min.append(-1000)
+						data_max.append(1000)
+				else:
+					# Assume it's the magnitude
+					labels[i][j][k] += r': $\kappa_{\rm PU}$'
+					ylabels[i][j] = "Correction"
+					if j == 0 and k == 0:
+						data_min.append(0.9)
+						data_max.append(1.1)
+			elif 'kc' in txt_list[i][j][k] or 'KC' in txt_list[i][j][k] or 'kappac' in txt_list[i][j][k] or 'KAPPAC' in txt_list[i][j][k] or 'kappa_c' in txt_list[i][j][k] or 'KAPPA_C' in txt_list[i][j][k] or 'kappa_C' in txt_list[i][j][k]:
+				labels[i][j][k] += r': $\kappa_{\rm C}$'
+				ylabels[i][j] = "Correction"
+				if j == 0 and k == 0:
+					data_min.append(0.95)
+					data_max.append(1.05)
+			elif 'fc' in txt_list[i][j][k] or 'FC' in txt_list[i][j][k] or 'f_c' in txt_list[i][j][k] or 'F_C' in txt_list[i][j][k]:
+				labels[i][j][k] += r': $f_{\rm cc}$'
+				ylabels[i][j] = "Frequency (Hz)"
+				if j == 0 and k == 0:
+					data_min.append(390)
+					data_max.append(450)
+			elif ('fs_over_Q' in txt_list[i][j][k] or 'FS_OVER_Q' in txt_list[i][j][k] or 'f_s_over_Q' in txt_list[i][j][k] or 'F_S_OVER_Q' in txt_list[i][j][k]):
+				labels[i][j][k] += r': $f_{\rm s} / Q$'
+				ylabels[i][j] = "Frequency (Hz)"
+				if j == 0 and k == 0:
+					data_min.append(-2)
+					data_max.append(0)
+			elif 'fs' in txt_list[i][j][k] or 'FS' in txt_list[i][j][k] or 'f_s' in txt_list[i][j][k] or 'F_S' in txt_list[i][j][k]:
+				labels[i][j][k] += r': $f_{\rm s}^2$'
+				ylabels[i][j] = "Square Frequency (Hz$^2$)"
+				if j == 0 and k == 0:
+					data_min.append(-10)
+					data_max.append(200)
+			elif 'Q' in txt_list[i][j][k] or 'q' in txt_list[i][j][k]:
+				labels[i][j][k] += r': $Q^{-1}$'
+				ylabels[i][j] = "Inverse Quality Factor"
+				if j == 0 and k == 0:
+					data_min.append(-1)
+					data_max.append(1)
+			else:
+				labels[i][j][k] += 'TDCF'
+				ylabels[i][j] = "Correction"
+				if j == 0 and k == 0:
+					data_min.append(min(data[i][j][k][1]))
+					data_max.append(max(data[i][j][k][1]))
+				print('unknown legend label for %s' % txt_list[i][j][k])
+
+# Read data from files and plot it
+colors = [['red', 'darkred', 'salmon'], ['limegreen', 'darkgreen', 'lightgreen'], ['blue', 'darkblue', 'lightblue']] # Hopefully the user will not compare more than 3 versions of calibration or plot more than 3 TDCFs on the same plot.
+
+# Figure out the (horizontal) time axis
+dur = t_end - t_start
+t_unit = 'seconds'
+sec_per_t_unit = 1.0
+if dur > 60 * 60 * 100:
+	t_unit = 'days'
+	sec_per_t_unit = 60.0 * 60.0 * 24.0
+elif dur > 60 * 100:
+	t_unit = 'hours'
+	sec_per_t_unit = 60.0 * 60.0
+elif dur > 100:
+	t_unit = 'minutes'
+	sec_per_t_unit = 60.0
+
+# Make plots
+plt.figure(figsize = (18, len(data) * 6))
+for i in range(len(data)):
+	ax = plt.subplot(len(data), 1, i + 1)
+	for j in range(len(data[i])):
+		if ylabels[i][j] is not None:
+			plt.ylabel(ylabels[i][j])
+		for k in range(len(data[i][j])):
+			plt.plot((data[i][j][k][0] - t_start) / sec_per_t_unit, data[i][j][k][1], colors[j % 3][k % 3], linewidth = 2.0, label = labels[i][j][k])
+	leg = plt.legend(fancybox = True)
+	leg.get_frame().set_alpha(0.5)
+
+	if i == 0 and options.plot_title is not None:
+		plt.title(options.plot_title)
+
+	if i == len(data) - 1:
+		plt.xlabel(r'${\rm Time \  in \  %s \  since \  %s \  UTC}$' % (t_unit, time.strftime("%b %d %Y %H:%M:%S".replace(':', '{:}').replace('-', '\mbox{-}').replace(' ', '\ '), time.gmtime(t_start + 315964782))))
+	ticks_and_grid(ax, xscale = 'linear', yscale = 'linear', ymin = data_min[i], ymax = data_max[i])
+
+plt.savefig(options.filename.split('.')[0] + '_%d-%d.' % (t_start, dur) + options.filename.split('.')[1])
+
+
-- 
GitLab