added eooe calculation script

e7baef8d · Peter Jansweijer · 03f0809c · e7baef8d
Commit e7baef8d authored Nov 07, 2018 by Peter Jansweijer
Hide whitespace changes
Inline Side-by-side

Showing with 304 additions and 0 deletions

calc.py sw/eooe/calc.py +304 -0

No files found.
--- a/sw/eooe/calc.py
+++ b/sw/eooe/calc.py
+#!/usr/bin/python
+
+"""
+calc.py: Analyses the stat_result files that were output by AnalyzeScopeData.py
+and takes the linear correlation results to calculate the final optical to
+electrical delay of the photo detector.
+
+-------------------------------------------------------------------------------
+Copyright (C) 2017 Peter Jansweijer, Henk Peek
+
+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 3 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, see <http://www.gnu.org/licenses/>
+-------------------------------------------------------------------------------
+
+Usage:
+  calc.py     <root-name> (for example "180928_12_09_18_3_lambda_insitu_alpha_scan")
+  calc.py     -h | --help
+
+  <root-name> identifies the files to analyse. Four files have the root-name plus:
+              <root-name>_l1_a    lambda 1 (-> lambda 2)
+              <root-name>_l2_a    lambda 2
+              <root-name>_l2_b    lambda 2 (-> lambda 1)
+              <root-name>_l1_b    lambda 1
+              so two sets of files are recorded: 'a' setting first lambda 1 then lambda 2
+              and 'b' setting first lambda 2 then lambda 1. This is to evaluate whether
+              tuning in one or the other direction is of influence on the CRTT measurements.
+              So far this has not been an issue.
+
+Options:
+  -h --help    Show this screen.
+  -t           tolerance [ps], skip points outside fitted line, default: "200"
+  -i           use ITU channels instead of wavelength
+
+"""
+
+import os
+import sys
+import scipy
+import numpy
+import matplotlib.pyplot as plt
+
+import pdb
+
+# Add parent directory (containing 'lib') to the module search path
+lib_path = (os.path.dirname(os.path.abspath(__file__)))
+lib_path = os.path.join(lib_path,"..")
+sys.path.insert(0,lib_path)
+
+############################################################################
+
+
+def read_stat_result(dirname,sign="neg"):
+  """
+  read_stat_result
+  
+  dirname          -- pointer to the directory that contains the stat_result.txt file
+  sign             -- default "neg" i.e. stat_result lineair-mean values are negative due
+                   -- to the fact that channel 1 and 3 of the oscilloscope are swapped.
+
+  returns:
+  fail, lin_meas 
+                   -- <Bool> fail = True when an error occurred while reading the file
+                   -- <Tuple> lin_meas [0] = lineair-mean,
+                   --                  [1] = se,
+                   --                  [2] = sem,
+  """
+
+  if not(os.path.isdir(dirname)):
+    print(dirname+": is not a directory.")
+    return (True,nan)
+
+  stat_result_filename = os.path.join(dirname,"stat_result.txt")
+  if not(os.path.isfile(stat_result_filename)):
+    print(stat_result_filename + ": No stat_result.txt file found.")
+    return (True,nan)
+  
+  stat_result_file = open(stat_result_filename,"r")
+  while 1:
+    # Lines look like this:
+    # data/data_181106_1/181106_2 lineair-mean: -6.79100959999e-10 se: 2.16707944778e-14 sem: 4.33415889555e-15 topt: 1     if line == "":
+    line = stat_result_file.readline()
+    if line == "":
+      break #end of file
+    line_lst = line.strip().split(" ")
+
+    if "lineair" in line_lst[1]:
+      if sign == "neg":
+        lin_meas = -float(line_lst[2]), float(line_lst[4]), float(line_lst[6])
+      else:
+        lin_meas = +float(line_lst[2]), float(line_lst[4]), float(line_lst[6])
+      return (False,lin_meas)
+
+  return (True,nan)
+
+###########################################################################
+  
+  
+def write_meas(file_ptr,meas_txt,meas):
+  """
+  Test whether a measurement is present or not
+  If not then exit with a failure message
+  
+  file_ptr         -- <file> Pointer to output file
+  meas_txt         -- <Str> preceding string
+  meas             -- <Tuple> lin_meas [0] = lineair-mean,
+                   --                  [1] = se,
+                   --                  [2] = sem,
+
+  returns:
+  """
+
+  file_ptr.write(meas_txt + str(meas[0]) + ", " + str(meas[1]) + ", " + str(meas[2]) + "\n")
+  return
+
+###########################################################################
+#
+# If run from commandline, we can test the library
+#a_meas[
+"""
+Usage:
+  calc.py     <root-name> (for example "180928_12_09_18_3_lambda_insitu_alpha_scan")
+  calc.py     -h | --help
+
+  <root-name> identifies the files to analyse. Four files have the root-name plus:
+              <root-name>_l1_a    lambda 1 (-> lambda 2)
+              <root-name>_l2_a    lambda 2
+              <root-name>_l2_b    lambda 2 (-> lambda 1)
+              <root-name>_l1_b    lambda 1
+              so two sets of files are recorded: 'a' setting first lambda 1 then lambda 2
+              and 'b' setting first lambda 2 then lambda 1. This is to evaluate whether
+              tuning in one or the other direction is of influence on the CRTT measurements.
+              So far this has not been an issue.
+
+Options:
+  -h --help    Show this screen.
+  -t           tolerance [ps], skip points outside fitted line, default: "200"
+  -i           use ITU channels instead of wavelength
+
+"""
+if __name__ == "__main__":
+
+  import argparse
+  parser = argparse.ArgumentParser()
+  parser.add_argument("name", help="directory that contains the measurement files _a, _b1, _b2, _c, _d1, _d2, _e, _f")
+  args = parser.parse_args()
+
+  name = args.name
+
+  files = []
+
+  # Figure out if the input files exist. Extensions should be:
+  files_exist = True
+  
+  if not(os.path.isdir(name)):
+    print(name+": is not a directory.")
+    print("Provide the root that contains measurement subdirectories _a, _b1, _b2, _c, _d1, _d2, _e, _f")
+    sys.exit()
+  
+  meas_dirs = os.listdir(name)  
+
+  fail = True
+
+  for meas_dir in meas_dirs:
+    #print(os.path.join(name,meas_dir))
+    if meas_dir[-2:] == "_aux":
+      fail,aux_meas = read_stat_result(os.path.join(name,meas_dir))
+    elif meas_dir[-2:] == "_a":
+      fail,a_meas = read_stat_result(os.path.join(name,meas_dir))
+    elif meas_dir[-2:] == "b1":
+      fail,b1_meas = read_stat_result(os.path.join(name,meas_dir))
+    elif meas_dir[-2:] == "b2":
+      fail,b2_meas = read_stat_result(os.path.join(name,meas_dir))
+    elif meas_dir[-2:] == "_c":
+      fail,c_meas = read_stat_result(os.path.join(name,meas_dir))
+    elif meas_dir[-2:] == "d1":
+      fail,d1_meas = read_stat_result(os.path.join(name,meas_dir))
+    elif meas_dir[-2:] == "d2":
+      fail,d2_meas = read_stat_result(os.path.join(name,meas_dir))
+    elif meas_dir[-2:] == "_e":
+      fail,e_meas = read_stat_result(os.path.join(name,meas_dir))
+    elif meas_dir[-2:] == "_f":
+      fail,f_meas = read_stat_result(os.path.join(name,meas_dir))
+    
+    if fail:
+      print(" failed at file: " + os.path.join(name,meas_dir))
+      sys.exit()
+
+  # Check if all mandatory measurements are present:
+  try:
+    a_meas[0]
+  except:
+    print("FAIL: no file found for measurements _a")
+    sys.exit()
+
+  try:
+    b1_meas[0]
+  except:
+    print("FAIL: no file found for measurements _b1") 
+    sys.exit()
+
+  try:
+    b2_meas[0]
+  except:
+    print("FAIL: no file found for measurements _b2")
+    sys.exit()
+
+  try:
+    c_meas[0]
+  except:
+    print("FAIL: no file found for measurements _c")
+    sys.exit()
+
+  try:
+    d1_meas[0]
+  except:
+    print("FAIL: no file found for measurements _d1")
+    sys.exit()
+
+  try:
+    d2_meas[0]
+  except:
+    print("FAIL: no file found for measurements _d2")
+    sys.exit()
+
+  try:
+    e_meas[0]
+  except:
+    print("FAIL: no file found for measurements _e")
+    sys.exit()
+
+  try:
+    f_meas[0]
+  except:
+    print("FAIL: no file found for measurements _f")
+    sys.exit()
+   
+  # Now do the calculations:
+
+  # lineair-mean values are tuple [0] items
+  # lineair-se are tuple [1] items, lineair-sem are tuple[2] items
+  D_f = (85.36e-12,1.1e-13,1.1e-14)
+  Type_B_err = (1.1e-12,0,1.67929e-12)
+  # note that b and d are two measurements; take the mean value and
+  # calculate their standard deviation and standard error.
+  b_meas = ((b1_meas[0] + b2_meas[0])/2, numpy.sqrt(b1_meas[1]**2 + b2_meas[1]**2)/2, numpy.sqrt(b1_meas[2]**2 + b2_meas[2]**2)/2)
+  d_meas = ((d1_meas[0] + d2_meas[0])/2, numpy.sqrt(d1_meas[1]**2 + d2_meas[1]**2)/2, numpy.sqrt(d1_meas[2]**2 + d2_meas[2]**2)/2)
+
+  # https://doi.org/10.1364/OE.26.014650
+  # fill in Eq. (10) 
+  
+  D_OE = ((-a_meas[0] - b_meas[0] + c_meas[0] - d_meas[0] + e_meas[0] + f_meas[0] + D_f[0]) / 2) - Type_B_err[0]
+
+  # Eq. (10) fraction errors:
+  frac_se  = numpy.sqrt(a_meas[1]**2 + b_meas[1]**2 + c_meas[1]**2 + d_meas[1]**2 + e_meas[1]**2 + f_meas[1]**2 + D_f[1]**2) / 2
+  frac_sem = numpy.sqrt(a_meas[2]**2 + b_meas[2]**2 + c_meas[2]**2 + d_meas[2]**2 + e_meas[2]**2 + f_meas[2]**2 + D_f[2]**2) / 2
+  
+  # Eq. (10) total errors:
+  D_OE_se  = numpy.sqrt(frac_se**2 + Type_B_err[1]**2)
+  D_OE_sem = numpy.sqrt(frac_sem**2 + Type_B_err[2]**2)
+
+  result_file = open(os.path.join(name,"calc_result.txt"),"w")
+  result_file.write("# Calculated result for eooe measurement\n")
+  result_file.write("# file: " + os.path.join(name,"calc_result.txt") + "\n")
+  write_meas(result_file, "a,        ", a_meas)
+  write_meas(result_file, "b,        ", b_meas)
+  write_meas(result_file, "c,        ", c_meas)
+  write_meas(result_file, "d,        ", d_meas)
+  write_meas(result_file, "e,        ", e_meas)
+  write_meas(result_file, "f,        ", f_meas)
+  result_file.write("\n")
+  result_file.write("D_OE,     " + str(D_OE) + "\n")
+  result_file.write("D_OE_se,  " + str(D_OE_se) + "\n")
+  result_file.write("D_OE_sem, " + str(D_OE_sem) + "\n")
+  result_file.close()
+
+  print("a:  ", a_meas)
+  print("b: ", b_meas)
+  #print("b1: ", b1_meas)
+  #print("b2: ",b2_meas)
+  print("c:  ",c_meas)
+  print("d: ",d_meas)
+  #print("d1: ",d1_meas)
+  #print("d2: ",d2_meas)
+  print("e:  ",e_meas)
+  print("f:  ",f_meas)
+
+
+  print("D_OE", D_OE)
+  #print("frac_se", frac_se)
+  #print("frac_sem", frac_sem)
+  print("D_OE_se", D_OE_se)
+  print("D_OE_sem", D_OE_sem)
+  
+  sys.exit()