Skip to content

W
White Rabbit Calibration
Commit 7ce67bc7 authored by Peter Jansweijer's avatar Peter Jansweijer
Browse files
Options

calculate std-dev for alpha and print in _err.out file

parent a726052b
Hide whitespace changes
Inline Side-by-side
Showing with 92 additions and 38 deletions
sw/insitu_alpha/analyze_sellmeier.py
View file @ 7ce67bc7
......@@ -248,7 +248,7 @@ def calc_alpha_tangent(l, crtt_l, tangent, fixed_lambda):
############################################################################
def calc_alpha_sellmeier(l, fixed_lambda, popt):
def calc_alpha_sellmeier(l, fixed_lambda, A, B, C, D, E):
"""
calc_alpha_sellmeier calculates the alpha factor for the wavelength lambda l
while using fixed_lambda for the return path.
......@@ -257,7 +257,7 @@ def calc_alpha_sellmeier(l, fixed_lambda, popt):
l -- <float> lambda [m]
fixed_lambda -- <float> fixed return channel lambda [m]
a..e -- <float> Sellmeier terms
A..E -- <float> Sellmeier terms
returns:
alpha -- <float> alpha value for lambda l when using fixed_lambda
......@@ -265,12 +265,6 @@ def calc_alpha_sellmeier(l, fixed_lambda, popt):
"""
A = popt[0]
B = popt[1]
C = popt[2]
D = popt[3]
E = popt[4]
bb = l**2 - fixed_lambda**2
cc = l**-2 - fixed_lambda**-2
dd = l**4 - fixed_lambda**4
......@@ -281,6 +275,54 @@ def calc_alpha_sellmeier(l, fixed_lambda, popt):
return(alpha)
############################################################################
def calc_alpha_sellmeier_stddev(l, l_stddev, fixed_lambda, fixed_lambda_stddev, popt, perr):
"""
calc_alpha_sellmeier calculates the alpha factor for the wavelength lambda l
while using fixed_lambda for the return path.
It takes the Sellmeier cooefcients that came out of a cable rount trip time
wavelength scan.
l -- <float> lambda [m]
l_stddev -- <float> std-dev for lambda [m]
fixed_lambda -- <float> fixed return channel lambda [m]
fixed_lambda_stddev -- <float> std-dev for fixed return channel lambda [m]
popt -- <float> Sellmeier terms
perr -- <float> std-dev for Sellmeier terms
returns:
alpha -- <float> alpha value for lambda l when using fixed_lambda
as return channel
alpha -- <float> std-dev of alpha value
"""
A = popt[0]
B = popt[1]
C = popt[2]
D = popt[3]
E = popt[4]
# Calculate the (mean) value of alpha
alpha = calc_alpha_sellmeier(l, fixed_lambda, A, B, C, D, E)
# Calculate the stddev contribution to alpha of each of the parameters l, fixed_lambda and perr
alpha_stddev_l = calc_alpha_sellmeier(l + l_stddev, fixed_lambda, A, B, C, D, E) - alpha
alpha_stddev_fixed_lambda = calc_alpha_sellmeier(l, fixed_lambda + fixed_lambda_stddev, A, B, C, D, E) - alpha
alpha_stddev_A = calc_alpha_sellmeier(l, fixed_lambda, A + perr[0], B, C, D, E) - alpha
alpha_stddev_B = calc_alpha_sellmeier(l, fixed_lambda, A, B + perr[1], C, D, E) - alpha
alpha_stddev_C = calc_alpha_sellmeier(l, fixed_lambda, A, B, C + perr[2], D, E) - alpha
alpha_stddev_D = calc_alpha_sellmeier(l, fixed_lambda, A, B, C, D + perr[3], E) - alpha
alpha_stddev_E = calc_alpha_sellmeier(l, fixed_lambda, A, B, C, D, E + perr[4]) - alpha
# Calculate the RMS stddev of alpha
alpha_stddev = numpy.sqrt(alpha_stddev_l**2 + alpha_stddev_fixed_lambda **2 + alpha_stddev_A **2 + alpha_stddev_B **2 + alpha_stddev_C **2 + alpha_stddev_D **2 + alpha_stddev_E **2)
return(alpha,alpha_stddev)
############################################################################
def calc_alpha_error(alpha1, alpha2, crtt_fixed_lambda):
......@@ -353,15 +395,18 @@ def group_delay(l, A, B, C, D, E):
############################################################################
def alpha_err_writeln(alpha_err_file, l1, l2, l_fix, l_idx, popt, crtt_fixed_lambda, alpha_3wl, time_err):
def alpha_err_writeln(alpha_err_file, l1, l_stddev, l2, l_fix, l_fix_stddev, l_idx, popt, perr, crtt_fixed_lambda, alpha_3wl, time_err):
"""
write a line to the alpha_err file in a pre formatted way.
alpha_err_file -- output file pointer
l1,l2,l_fix -- <float> 1, 2 and fixed lambda's in [nm].
l_stddev -- <float> std-dev for lambda [m]
l_fix_stddev -- <float> std-dev for fixed return channel lambda [m]
l_idx -- <tuple> (idx1,idx2) lambda index into two dimensional array of alpha's and errors
A..E -- <float> sellmeier constants
popt -- <float> sellmeier constants
perr -- <float> std-dev for sellmeier constants
crtt_fixed_lambda -- <float> round trip time @ fixed lambda in [ps]
alpha_3wl -- <2D array of float> of alpha values determined via 3-wavelength method at lambda's l1 and l2
(axes: idx1 => l1, idx2 => l2)
......@@ -374,13 +419,15 @@ def alpha_err_writeln(alpha_err_file, l1, l2, l_fix, l_idx, popt, crtt_fixed_lam
l2_str = "{0:.2f}".format(l2)
crtt_l1_str = "{0:.0f}".format(group_delay(l1, popt[0], popt[1], popt[2], popt[3], popt[4]))
crtt_l2_str = "{0:.0f}".format(group_delay(l2, popt[0], popt[1], popt[2], popt[3], popt[4]))
alpha_sel_str = "{0:.3e}".format(calc_alpha_sellmeier(l1/1e9, l_fix , popt))
alpha_sel, alpha_stddev_sel = calc_alpha_sellmeier_stddev(l1/1e9, l_stddev, l_fix, l_fix_stddev, popt, perr)
alpha_sel_str = "{0:.3e}".format(alpha_sel)
alpha_stddev_sel_str = "{0:.3e}".format(alpha_stddev_sel)
alpha_3wl_str = "{0:.3e}".format(alpha_3wl[l_idx])
alpha_spec_3wl_str = "{0:.0f}".format(spec_alpha(alpha_3wl[l_idx]))
time_err_str = "{0:.3e}".format(time_err[l_idx])
time_err_rel_str = "{0:.3e}".format(time_err[l_idx]/crtt_fixed_lambda)
alpha_err_file.write(l1_str + "\t" + l2_str + "\t" + crtt_l1_str + "\t" + crtt_l2_str + "\t" + alpha_sel_str + "\t" + alpha_3wl_str + "\t" + alpha_spec_3wl_str + "\t\t\t" + time_err_str + "\t" + time_err_rel_str + "\n")
alpha_err_file.write(l1_str + "\t" + l2_str + "\t" + crtt_l1_str + "\t" + crtt_l2_str + "\t" + alpha_sel_str + "\t" + alpha_stddev_sel_str + "\t" + alpha_3wl_str + "\t" + alpha_spec_3wl_str + "\t\t\t" + time_err_str + "\t" + time_err_rel_str + "\n")
return
......@@ -444,7 +491,7 @@ def average_insitu_files(files):
############################################################################
def analyze_plot(insitu_file, analyse_single, x, y, name, tolerance, use_itu_channels, ref_name, fixed_itu_channel, tangent_array, temperature_array):
def analyze_plot(insitu_file, analyse_single, x, y, name, tolerance, use_itu_channels, ref_name, fixed_itu_channel, l_stddev, fixed_lambda_stddev, tangent_array, temperature_array):
"""
Reads all insitu_files from a directory and calculates the average of all measurements.
......@@ -458,6 +505,8 @@ def analyze_plot(insitu_file, analyse_single, x, y, name, tolerance, use_itu_cha
use_itu_channels -- <boolean> use ITU channels instead of wavelength
ref_name -- <string> file or directory name containing crtt reference measurements
fixed_itu_channel -- <int> channel used for fixed Slave-to-Master communication
l_stddev -- <float> std-dev for lambda [m]
fixed_lambda_stddev -- <float> std-dev for fixed return channel lambda [m]
tangent_array -- <numpy.array> keep tangent values in an array in case multiple files are processed
temperature_array -- <numpy.array> keep temperature values in an array in case multiple files are processed
......@@ -574,7 +623,7 @@ def analyze_plot(insitu_file, analyse_single, x, y, name, tolerance, use_itu_cha
ly_clean.append(clean_y) # 5-t Sellmeier
alpha_tan = calc_alpha_tangent(i/1e9, clean_y, tangent*1e9, fixed_lambda)
alpha_sel = calc_alpha_sellmeier(i/1e9, fixed_lambda, popt)
alpha_sel, alpha_stddev_sel = calc_alpha_sellmeier_stddev(i/1e9, l_stddev, fixed_lambda, fixed_lambda_stddev, popt, perr)
alpha_tangent.append(alpha_tan)
#alpha_tangent.append(calc_alpha_tangent(i/1e9, clean_y, dispersion, fixed_lambda))
......@@ -706,7 +755,7 @@ def analyze_plot(insitu_file, analyse_single, x, y, name, tolerance, use_itu_cha
crtt_l1 = group_delay(l1,popt[0],popt[1],popt[2],popt[3],popt[4])
crtt_l2 = group_delay(l2,popt[0],popt[1],popt[2],popt[3],popt[4])
alpha_3wl = calc_alpha_l1(l1/1e9, l2/1e9, crtt_l1, crtt_l2, fixed_lambda)
alpha_3wlsel = calc_alpha_sellmeier(l1/1e9, fixed_lambda, popt)
alpha_3wlsel, alpha_stddev_3wlsel = calc_alpha_sellmeier_stddev(l1/1e9, l_stddev, fixed_lambda, fixed_lambda_stddev, popt, perr)
alpha_err = alpha_3wlsel - alpha_3wl
time_err = calc_alpha_error(alpha_3wlsel, alpha_3wl, crtt_fixed_lambda)
time_err_relative = time_err/crtt_fixed_lambda
......@@ -816,7 +865,7 @@ def analyze_plot(insitu_file, analyse_single, x, y, name, tolerance, use_itu_cha
#for i in rng_sel:
# alpha_err_file.write(str(l2[i[0]][0]) + "\t" + str(l2[i[1]][0]) + "\n")
alpha_err_file.write("l1[nm]\tl2[nm]\tcrtt_(l1)[ps]\tcrtt(l2)[ps]\talpha_sel(l1)\talpha_3wl(l1)\tSPEC_alpha_3wl(l1)\ttime_err[ps]\ttime_relative_err[t_err ps/link_length ps]\n")
alpha_err_file.write("l1[nm]\tl2[nm]\tcrtt_(l1)[ps]\tcrtt(l2)[ps]\talpha_sel(l1)\talpha_stddev_sel(l1)\talpha_3wl(l1)\tSPEC_alpha_3wl(l1)\ttime_err[ps]\ttime_relative_err[t_err ps/link_length ps]\n")
alpha_ddelay_file.write("# lambda-1[nm]\tITU-ch\tSPEC_alpha_3wl(l1)\tlambda-2[nm]\t\ttime_err[ps]\n")
for i in rng_sel:
......@@ -827,7 +876,7 @@ def analyze_plot(insitu_file, analyse_single, x, y, name, tolerance, use_itu_cha
itu_ch = itu_conv.wavelength_2_itu(lambda1/1e9)
time_err_str = "{0:.3e}".format(time_err[lambda_idx]*1e-12)
alpha_err_writeln(alpha_err_file, lambda1, lambda2, fixed_lambda, lambda_idx, popt, crtt_fixed_lambda, alpha_3wl, time_err)
alpha_err_writeln(alpha_err_file, lambda1, l_stddev, lambda2, fixed_lambda, fixed_lambda_stddev, lambda_idx, popt, perr, crtt_fixed_lambda, alpha_3wl, time_err)
alpha_ddelay_file.write(str(lambda1)+", "+str(itu_ch)+", "+str(alpha_spec_3wl)+", "+str(lambda2)+", "+time_err_str+"\n")
alpha_err_file.write("============================\n")
......@@ -847,7 +896,7 @@ def analyze_plot(insitu_file, analyse_single, x, y, name, tolerance, use_itu_cha
itu_ch = itu_conv.wavelength_2_itu(lambda1/1e9)
time_err_str = "{0:.3e}".format(time_err[lambda_idx]*1e-12)
alpha_err_writeln(alpha_err_file, lambda1, lambda2, fixed_lambda, lambda_idx, popt, crtt_fixed_lambda, alpha_3wl, time_err)
alpha_err_writeln(alpha_err_file, lambda1, l_stddev, lambda2, fixed_lambda, fixed_lambda_stddev, lambda_idx, popt, perr, crtt_fixed_lambda, alpha_3wl, time_err)
alpha_ddelay_file.write(str(lambda1)+", "+str(itu_ch)+", "+str(alpha_spec_3wl)+", "+str(lambda2)+", "+time_err_str+"\n")
alpha_err_file.write("============================\n")
......@@ -858,9 +907,8 @@ def analyze_plot(insitu_file, analyse_single, x, y, name, tolerance, use_itu_cha
lambda1 = l1[l1_idx][l2_idx]
lambda2 = l2[l1_idx][l2_idx]
lambda_idx = l1_idx,l2_idx
alpha_err_writeln(alpha_err_file, lambda1, lambda2, fixed_lambda, lambda_idx, popt, crtt_fixed_lambda, alpha_3wl, time_err)
alpha_err_writeln(alpha_err_file, lambda1, l_stddev, lambda2, fixed_lambda, fixed_lambda_stddev, lambda_idx, popt, perr, crtt_fixed_lambda, alpha_3wl, time_err)
#alpha_err_file.write(str(l1[l1_idx][l2_idx]) + "\t" + str(l2[l1_idx][l2_idx]) + "\t" + str(group_delay(l1[l1_idx][l2_idx],popt[0],popt[1],popt[2],popt[3],popt[4])) + "\t " + str(group_delay(l2[l1_idx][l2_idx],popt[0],popt[1],popt[2],popt[3],popt[4])) + "\t" + str(calc_alpha_sellmeier(l1[l1_idx][l2_idx]/1e9, fixed_lambda,popt[0],popt[1],popt[2],popt[3],popt[4])) + "\t " + str(alpha_3wl[l1_idx][l2_idx]) + "\t " + str(spec_alpha(alpha_3wl[l1_idx][l2_idx])) + "\t" + str(time_err[l1_idx][l2_idx]) + "\t" + str(time_err[l1_idx][l2_idx]/crtt_fixed_lambda) + "\n")
alpha_err_file.close()
alpha_ddelay_file.close()
......@@ -921,6 +969,11 @@ if __name__ == "__main__":
ref_name = args.r
fixed_itu_channel = 20
# Define 3 pico-meter standard deviations for wavelengths
# (= resolution of the ANDO AQ6140)
l_stddev = 3e-12 # 3 [pm]
fixed_lambda_stddev = 3e-12 # 3 [pm]
analyse_single = False
ref_analyse_single = False
files = []
......@@ -1003,6 +1056,7 @@ if __name__ == "__main__":
axs = fig_sell_fits.add_subplot(111)
axs.set_title("Sellmeier fits")
axs.set_xlabel('Wavelength [nm]')
axs.set_ylabel('CRTT [ps]')
if measurementfile == True:
print(insitu_file)
......@@ -1027,7 +1081,7 @@ if __name__ == "__main__":
t = crtt_array["temp"]
tangent_array, temperature_array, sellmeier_fit, result_line = analyze_plot(insitu_file, analyse_single, x, y, name, tolerance, use_itu_channels, ref_name, fixed_itu_channel, tangent_array, temperature_array)
tangent_array, temperature_array, sellmeier_fit, result_line = analyze_plot(insitu_file, analyse_single, x, y, name, tolerance, use_itu_channels, ref_name, fixed_itu_channel, l_stddev, fixed_lambda_stddev, tangent_array, temperature_array)
if not(analyse_single):
numb_of_meas = numb_of_meas + 1
result_file.write(result_line)
......@@ -1067,28 +1121,28 @@ if __name__ == "__main__":
else:
insitu_file = name + "/averaged"
# Process the averaged measurements as a single measurement
analyse_single = True
# Process the averaged measurements as a single measurement
analyse_single = True
if ref_name != '':
y = crtt_array - ref_crtt_array
y_stdev = numpy.sqrt(crtt_array_std**2 + ref_crtt_array_std**2)
else:
y = crtt_array
if ref_name != '':
y = crtt_array - ref_crtt_array
y_stdev = numpy.sqrt(crtt_array_std**2 + ref_crtt_array_std**2)
else:
y = crtt_array
t = temp_array
t = temp_array
tangent_array, temperature_array, sellmeier_fit, result_line = analyze_plot(insitu_file, analyse_single, x, y, name, tolerance, use_itu_channels, ref_name, fixed_itu_channel, tangent_array, temperature_array)
tangent_array, temperature_array, sellmeier_fit, result_line = analyze_plot(insitu_file, analyse_single, x, y, name, tolerance, use_itu_channels, ref_name, fixed_itu_channel, l_stddev, fixed_lambda_stddev, tangent_array, temperature_array)
result_file.write("\n")
result_file.write("Wavelength, Mean_CRTT, StdDev for CRTT per Averaged Wavelength for Reference setup\n")
for i in range(len(x)):
result_file.write(str(x[i]) + ", " + str(ref_crtt_array[i]) + ", " + str(ref_crtt_array_std[i]) + "\n")
result_file.write("\n")
result_file.write("Wavelength, Mean_CRTT, StdDev for CRTT per Averaged Wavelength for Reference setup\n")
for i in range(len(x)):
result_file.write(str(x[i]) + ", " + str(ref_crtt_array[i]) + ", " + str(ref_crtt_array_std[i]) + "\n")
result_file.write("\n")
result_file.write("Wavelength, Mean CRTT, StdDev for CRTT per Averaged Wavelength for Measurement setup\n")
for i in range(len(x)):
result_file.write(str(x[i]) + ", " + str(crtt_array[i]) + ", " + str(crtt_array_std[i]) + "\n")
result_file.write("\n")
result_file.write("Wavelength, Mean CRTT, StdDev for CRTT per Averaged Wavelength for Measurement setup\n")
for i in range(len(x)):
result_file.write(str(x[i]) + ", " + str(crtt_array[i]) + ", " + str(crtt_array_std[i]) + "\n")
sys.exit()
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment