added calc_alpha_l1_stddev

d20edd2f · Peter Jansweijer · d1726297 · d20edd2f
Commit d20edd2f authored Apr 05, 2019 by Peter Jansweijer
Hide whitespace changes
Inline Side-by-side

Showing with 84 additions and 5 deletions

analyze_sellmeier.py sw/insitu_alpha/analyze_sellmeier.py +84 -5

No files found.
--- a/sw/insitu_alpha/analyze_sellmeier.py
+++ b/sw/insitu_alpha/analyze_sellmeier.py
@@ -222,6 +222,43 @@ def calc_alpha_l1(lambda_1, lambda_2, crtt_lambda_1, crtt_lambda_2, fixed_lambda

 ############################################################################

+def calc_alpha_l1_stddev(lambda_1, lambda_1_stddev, lambda_2, lambda_2_stddev, crtt_lambda_1, crtt_lambda_1_stddev, crtt_lambda_2, crtt_lambda_2_stddev, fixed_lambda, fixed_lambda_stddev):
+  """
+  calc_alpha_l1 calculates the alpha factor for wavelength lambda l1.
+  For this caluculation we need two cable round trip times that were measured
+  using lamdba l1 and lambda l2 while both used a common fixed_lambda for the
+  return channel.
+
+  lambda_1, lambda_2, fixed_lambda  -- <float> lambda's in [m]
+  lambda_1_stddev, lambda_2_std_dev,
+  fixed_lambda_stddev               -- <float> std-dev for lambda's in [m]
+  crtt_l1, crtt_l2                  -- <float> cable round trip times [ps]measured
+                                               at lambda_1 and lambda_2
+  crtt_l1_stddev, crtt_l2_stddev    -- <float> std-dev for CRTT [ps]measured
+
+  returns:
+  alpha_lambda_1                    -- <float>
+  alpha_lambda_1_stddev             -- <float>
+
+  """
+
+  # Calculate the (mean) value of alpha_lambda_1
+  alpha_lambda_1 = calc_alpha_l1(lambda_1, lambda_2, crtt_lambda_1, crtt_lambda_2, fixed_lambda)
+
+  # Calculate the stddev contribution to alpha_lambda_1 for each of the parameters lambda_1, lambda_2, crtt_lambda_1, crtt_lambda_2 and fixed_lambda
+  alpha_stddev_lambda1      = calc_alpha_l1(lambda_1 + lambda_1_stddev, lambda_2, crtt_lambda_1, crtt_lambda_2, fixed_lambda) - alpha_lambda_1
+  alpha_stddev_lambda2      = calc_alpha_l1(lambda_1, lambda_2 + lambda_2_stddev, crtt_lambda_1, crtt_lambda_2, fixed_lambda) - alpha_lambda_1
+  alpha_stddev_crtt_lambda1 = calc_alpha_l1(lambda_1, lambda_2, crtt_lambda_1 + crtt_lambda_1_stddev, crtt_lambda_2, fixed_lambda) - alpha_lambda_1
+  alpha_stddev_crtt_lambda2 = calc_alpha_l1(lambda_1, lambda_2, crtt_lambda_1, crtt_lambda_2 + crtt_lambda_2_stddev, fixed_lambda) - alpha_lambda_1
+  alpha_stddev_fixed_lambda = calc_alpha_l1(lambda_1, lambda_2, crtt_lambda_1, crtt_lambda_2, fixed_lambda + fixed_lambda_stddev) - alpha_lambda_1
+
+  # Calculate the RMS stddev of alpha
+  alpha_lambda_1_stddev = numpy.sqrt(alpha_stddev_lambda1**2 + alpha_stddev_lambda2 **2 + alpha_stddev_crtt_lambda1 **2 + alpha_stddev_crtt_lambda2 **2 + alpha_stddev_fixed_lambda)
+
+  return (alpha_lambda_1, alpha_lambda_1_stddev)
+
+############################################################################
+
 def calc_alpha_tangent(l, crtt_l, tangent, fixed_lambda):
  """
  calc_alpha_tangent calculates the alpha factor for the wavelength lambda l.
@@ -421,6 +458,49 @@ def group_delay(l, A, B, C, D, E):
  
 ############################################################################

+def group_delay_stddev(l, l_stddev, popt, perr):
+  """
+  group_delay calculates the group delay for the wavelength lambda l based on
+  the 5-term Sellmeier equation with coefficients A..E.
+  
+  5-term sellmeier group delay:
+  tau(lambda) = A + B * lambda^2 + C * lambda^-2 + D * lambda^4 + E * lambda^-4
+
+  popt              -- <numpy.float64> Sellmeier coefficients
+  perr              -- <numpy.float64> std-dev from covariance matix pcov
+  l                 -- <float> lambda is in [nm].
+  l_stddev          -- <float> std-dev for lambda is in [nm].
+                    -- Note! Convert to [m] before executing sellmeier fit
+
+  returns:
+  delay             -- <float> group delay value based on selm5 and l
+  delay_stddev      -- <float> std_dev for delay
+  """
+
+  A = popt[0]
+  B = popt[1]
+  C = popt[2]
+  D = popt[3]
+  E = popt[4]
+
+  # Calculate the (mean) value of sellmeier_5
+  sellmeier_5 = group_delay(l, A, B, C, D, E)
+
+  # Calculate the stddev contribution to t_err of each of the parameters l, A to E
+  sellmeier_5_l_stddev = group_delay(l + l_stddev, A, B, C, D, E) - sellmeier_5
+  sellmeier_5_A        = group_delay(l, A + perr[0], B, C, D, E) - sellmeier_5
+  sellmeier_5_B        = group_delay(l, A, B + perr[1], C, D, E) - sellmeier_5
+  sellmeier_5_C        = group_delay(l, A, B, C + perr[2], D, E) - sellmeier_5
+  sellmeier_5_D        = group_delay(l, A, B, C, D + perr[3], E) - sellmeier_5
+  sellmeier_5_E        = group_delay(l, A, B, C, D, E + perr[4]) - sellmeier_5
+
+  # Calculate the RMS stddev of alpha
+  sellmeier_5_stddev = numpy.sqrt(sellmeier_5_l_stddev**2 + sellmeier_5_A **2 + sellmeier_5_B **2 + sellmeier_5_C **2 + sellmeier_5_D **2 + sellmeier_5_E **2)
+ 
+  return(sellmeier_5, sellmeier_5_stddev)
+  
+############################################################################
+
 # 5-term sellmeier 
 # Dispersion coefficient (is the derivative of group_delay equation)
 # D(lambda) = 2 * B * lambda - 2 * C * lambda^-3 + 4 * D * lambda^3 -4 * E * lambda^-5
@@ -793,13 +873,12 @@ def analyze_plot(insitu_file, analyse_single, x, y, name, tolerance, use_itu_cha

      # crtt_l1 and crtt_l2 are 2D arrays of <float> crtt in [ps]
      l1,l2 = numpy.meshgrid(x_clean, x_clean)
-      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)
+      crtt_l1, crtt_l1_stddev = group_delay_stddev(l1, l_stddev, popt, perr)
+      crtt_l2, crtt_l2_stddev = group_delay_stddev(l2, l_stddev, popt, perr)
+      alpha_3wl, alpha_3wl_stddev = calc_alpha_l1_stddev(l1/1e9, l_stddev, l2/1e9, l_stddev, crtt_l1, crtt_l1_stddev, crtt_l2, crtt_l2_stddev, fixed_lambda, l_stddev)
      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, time_err_stddev = calc_alpha_error_stddev(alpha_3wlsel, alpha_stddev_3wlsel, alpha_3wl, 0, crtt_fixed_lambda, l_stddev)
+      time_err, time_err_stddev = calc_alpha_error_stddev(alpha_3wlsel, alpha_stddev_3wlsel, alpha_3wl, alpha_3wl_stddev, crtt_fixed_lambda, l_stddev)
      time_err_relative = time_err/crtt_fixed_lambda
      
      #pdb.set_trace()