#!/usr/bin/env python # ----------------------------------------------------------------------------------- # # Python macro for reading text (ASCII) file from LabView for ball rolling down incline # and calculating the gravitational accceleration out of it # # Author: Troels C. Petersen (petersen@nbi.dk), # Vojtech Pacik (vojtech.pacik@nbi.ku.dk) # Date: 16th of November 2017 # ----------------------------------------------------------------------------------- # from __future__ import print_function, division # import numpy as np # Matlab like syntax for linear algebra and functions import matplotlib.pyplot as plt # Plots and figures like you know them from Matlab from iminuit import Minuit, describe, Struct from probfit import BinnedLH, Chi2Regression, Extended # Helper tool for fitting from math import sqrt, cos, sin, tan, atan, pi from scipy import stats plt.close('all') # to close all open figures r = np.random r.seed(42) # ----------------------------------------------------------------------------------- # def sqr( a ) : """ Return squared value """ return a*a # --------------------------------------------------------- # def nice_string_output(names, values, extra_spacing = 2): """ Function to create a nice string output for figures """ max_values = len(max(values, key=len)) max_names = len(max(names, key=len)) string = "" for name, value in zip(names, values) : string += "{0:s} {1:>{spacing}} \n".format(name, value, spacing = extra_spacing + max_values + max_names - len(name)) return string[:-2] # --------------------------------------------------------- # def GetMeanAndErrorOnMeanAndRMS( inlist = [] ) : """ Return list with mean and its error out of list of input values """ if len( inlist ) == 0 : print("WARNING: Called function with an empty list") return [-9999999999.9, -99999999999.9] elif len( inlist ) == 1 : print("WARNING: Called function with a list of length one") return [inlist[0], -99999999999.9] # Find the mean: mu = 0.0 for value in inlist : mu += value mu = float(mu) / len( inlist ) # Find the standard deviation (rms) and error on mean (emu): rms = 0.0 for value in inlist : rms += (value - mu)*(value - mu) rms = sqrt(float(rms) / (len(inlist) - 1)) emu = float(rms) / sqrt(len(inlist)) return [mu, emu, rms] # --------------------------------------------------------- # def GetMeanAndChi2( inlist = [] ) : """ Return a list containing weighted mean, its error, chi2, ndf and chi2 probability out of input lists in form of [ [value1, err1 ], [value2, err2], ... ] """ if len( inlist ) == 0 : print("WARNING: Called function with an empty list") return [-9999999999.9, -99999999999.9, -99999999999.9, -99999999999.9, -99999999999.9] elif len( inlist ) == 1 : print("WARNING: Called function with a list of length one") return [inlist[0][0], inlist[0][1], -99999999999.9, -99999999999.9, -99999999999.9] # Calculate the weighted mean and standard deviation sums: wmu = 0.0 sigmu = 0.0 for value in inlist : # You put the code! # Obtain the weighted mean and standard deviation results: # Calculate the Chi-Square and print: chi2_val = 0.0 for value in inlist : # You put the code! prob_chi2 = stats.chi2.sf(chi2_val, len(inlist)-1) return [float(wmu), float(sigmu), float(chi2_val), float(len(inlist)-1), float(prob_chi2)] # ----------------------------------------------------------------------------------- # # Input values & parameters & settings # ----------------------------------------------------------------------------------- # Blinding = True # Determining if random offset is added to result for blinding SavePlots = False # For now, don't save plots (once you trust your code, switch on) Verbose = True # For now, print a lot of output (once you trust your code, switch off) Nverbose = 5 # But only print a lot for the first 10 random numbers # ----------------------------------------------------------------------------------- # # Measured angles and estimated uncertianties (with goniometer / "angle-meter") [deg] # ----------------------------------------------------------------------------------- # # NOTE: Normal and reverse refer to the orientation of the EXPERIMENT. # Reg(ular) and Inv(erted) refer to the orientation of the GONIOMETER. # The difference between regular and inverted values comes from two orientation of goniometer. # Thus, the ~1 degree change is due to lacking calibration of the measuring device (goniometer). # For exactly this reason, we measure twice swithcing the goniometer direction. # Example results: angle_normal_reg = [ [ 12.70, 0.3 ], [ 12.60, 0.3 ], [ 12.75, 0.2 ], [ 12.50, 0.3 ] ] # NORMAL direction: calc_angle_normal_reg = GetMeanAndChi2(angle_normal_reg) mean_angle_normal_reg = calc_angle_normal_reg[0] emean_angle_normal_reg = max(0.1, calc_angle_normal_reg[1]) # I don't trust an error below 0.1 degree print(" Angle (normal, regular) = {0:7.2f} +- {1:4.2f} deg".format(mean_angle_normal_reg, emean_angle_normal_reg)) print(" Agreement between {0:1d} angle measurements: Prob(Chi2={1:4.1f}, Ndof={2:2d}) = {3:5.3f}".format(int(len(angle_normal_reg)), calc_angle_normal_reg[2], int(calc_angle_normal_reg[3]), calc_angle_normal_reg[4])) if (calc_angle_normal_reg[4] < 0.01) : print(" WARNING: Chi2 probability is alarmingly low! p = {0:8.6f}".format(calc_angle_normal_reg[4])) # ----------------------------------------------------------------------------------- # # Measured angle using trigonometry of the incline [mm] # ----------------------------------------------------------------------------------- # # Get average height: # ----------------------------------------------------------------------------------- # # Measured ball diameter [mm] (no calculation needed!): # ----------------------------------------------------------------------------------- # ball_diam_small = [ 10.00, 0.05 ] # small ball # Measured rail distance [mm] (no calculation needed!): # ----------------------------------------------------------------------------------- # rails = [ [ 6.05, 0.05 ], [ 6.00, 0.05 ], [ 6.10, 0.05 ], [ 6.05, 0.05 ], [ 6.05, 0.05 ] ] calc_rails = GetMeanAndChi2(rails) mean_rails = calc_rails[0] emean_rails = calc_rails[1] print("\n Average rail distance = {0:7.2f} +- {1:4.2f} mm".format(calc_rails[0], calc_rails[1])) print(" Agreement between {0:1d} rail measurements: Prob(Chi2={1:4.1f}, Ndof={2:2d}) = {3:5.3f}".format(int(len(rails)), calc_rails[2], int(calc_rails[3]), calc_rails[4])) if (calc_rails[4] < 0.01) : print(" WARNING: Chi2 probability is alarmingly low! p = {0:8.6f}".format(calc_rails[4])) # ----------------------------------------------------------------------------------- # # Measured positions of gates [mm] # ----------------------------------------------------------------------------------- # # Put d_gate1 = [] err_d_gate1 = [] d_gate1_final = 0.0 # You do the math... #---------------------------------------------------------------------------------- # Fit of acceleration: #---------------------------------------------------------------------------------- # List of input file names with voltages and times files = [ "data_Example/data_NormDir_MedBall1_SummaryTimes.dat", "data_Example/data_RevDir_MedBall1_SummaryTimes.dat" ] # Loop over different data files and perform acceleration fit: for filename in files : time = [] etime = [] # Read the passing times from file: print("\n ========== Reading passing times from file '{0:s}' ===========\n".format(filename)) with open( filename, 'r' ) as file : for line in file : line = line.strip() line = line.split() time.append(float(line[0])) etime.append(float(line[1])) file.close() # Given the passing times and distances, fit for the acceleration: etime = [ 0.0, 0.0, 0.0, 0.0, 0.0 ] # If one wants to disregard the uncertainty in the timing measurements. # acc, eacc = EstimateAcceleration(time, dist, etime, edist, filename) # This function you have to write! # The filename is included in order for the routine to give the plot file a logical name! #---------------------------------------------------------------------------------- # Possibly estimate systematic uncertainty between runs and delta angle correction # (from table mis-alignment). #---------------------------------------------------------------------------------- # The difference between the two accelerations is due to the table, which is not level. # The angle of the table can be infered from this constraint: # a_normal / sin(theta + dtheta) = a_reversed / sin(theta - dtheta) d2r = 2.0 * pi / 360.0 theta = 13.0 * d2r err_theta = 0.5 * d2r a_normal = 1.301 err_a_normal = 0.014 a_reverse = 1.423 err_a_reverse = 0.014 # In order to determine dtheta, a binary search is performed, as this the analytical # solution is "hard". # In order to get the uncertainty on dtheta, it is calculated many times for choices # of the input parameters (varied without their uncertainties). list_dtheta = [] for i in xrange( 1000 ) : # Pick a set of numbers consistent with the accelerations and the angle: an_used = r.normal( a_normal, err_a_normal ) ar_used = r.normal( a_reverse, err_a_reverse ) theta_used = r.normal( theta, err_theta ) # Start with a value for dtheta, and through a binary search get a good value: dist = 0.01 dtheta = 0.02 ddtheta = 0.01 N = 0 while (abs(dist) > 0.00000001 and N < 50) : dist = (an_used / sin(theta_used - dtheta)) - (ar_used / sin(theta_used + dtheta)) if (dist > 0.0) : dtheta -= ddtheta else : dtheta += ddtheta ddtheta = ddtheta / 2.0 N += 1 # print(" {0:3d}: Nite = {1:2d} dist = {2:12.9f} ddtheta = {3:11.9f} dtheta = {4:7.5f}".format(i, N, dist, ddtheta, dtheta)) list_dtheta.append(dtheta) # Make a histogram of the results: Nbins = 100 dtheta_min = 0.01 dtheta_max = 0.02 binwidth = (dtheta_max - dtheta_min) / Nbins fig, ax = plt.subplots(figsize=(12, 6)) hist_dtheta = ax.hist(list_dtheta, bins=Nbins, range=(dtheta_min, dtheta_max), histtype='step', linewidth=2, label=r"Values for $\Delta \theta$", color='blue') ax.set_xlabel(r"$\Delta \theta$ [rad]") ax.set_ylabel("Frequency / 0.0001 rad") ax.set_title(r"Values for $\Delta \theta$ (obtained numerically with variations of input)") ax.legend(loc='upper right') # The distribution is Gaussian (as expected), and in order to extract a value and uncertainty for dtheta, # one could just use the mean and RMS of this distribution: dtheta = np.array(list_dtheta).mean() err_dtheta = np.array(list_dtheta).std() print(" The final result for dtheta = {0:.5f} +- {1:.5f} rad ({2:.4f})".format(dtheta, err_dtheta, err_dtheta/dtheta)) # However, one can also fit it... def gauss_extended(x, N, mu, sigma): """Non-normalized Gaussian""" return N / np.sqrt(2 * np.pi) / abs(sigma) * np.exp(-(x - mu) ** 2 / 2.0 / sigma ** 2) # .... #---------------------------------------------------------------------------------- # Finally, calculate g (for medium ball size): #---------------------------------------------------------------------------------- # First print the input parameters and their errors (to check): #print(" Input parameters for: NORMAL direction REVERSE direction:") #print(" Acceleration (a): = {0:6.4f} +- {1:6.4f} ({2:6.4f}) = {3:6.4f} +- {4:6.4f} ({5:6.4f})".format(mean_a_normal, emean_a_normal, emean_a_normal/mean_a_normal, mean_a_reverse, emean_a_reverse, emean_a_reverse/mean_a_reverse)) #print(" Angle (theta'): = {0:6.4f} +- {1:6.4f} ({2:6.4f}) = {3:6.4f} +- {1:6.4f} ({4:6.4f})".format(thetap_normal, err_thetap, err_thetap/thetap_normal, thetap_reverse, err_thetap/thetap_reverse)) #print(" Ball radius (R): = {0:6.4f} +- {1:6.4f} ({2:6.4f}) = {0:6.4f} +- {1:6.4f} ({2:6.4f})".format(Rball, err_Rball, err_Rball/Rball)) #print(" Rail distance (d): = {0:6.4f} +- {1:6.4f} ({2:6.4f}) = {0:6.4f} +- {1:6.4f} ({2:6.4f})".format(mean_rails, emean_rails, emean_rails/mean_rails)) # Defining the rails-ball-factor as const (for convinience) # const = 1.0 + 2.0/5.0... const = 1.67 # Calculation of g g = a_reverse / sin(theta) * const # Blind the result (if wanted): if (Blinding) : g_offset_blinding = np.random.normal(0.0, 1.0) print("NOTE: Result is blinded!!!") else : g_offset_blinding = 0.0 # Use the SAME blinding offset, so that values can still be used for cross-check: g = g + g_offset_blinding err_g = 0.01 # Just a guess - you have all the tools to get this one yourself! print(" The final value for g is: g = {0:7.4f} +- {1:6.4f} m/s^2 ({2:6.4f})".format(g, err_g, err_g/g)) # ----------------------------------------------------------------------------------- # Comparison to top notch measurement (VERY PRECISE!) and models (not as precise): # ----------------------------------------------------------------------------------- g_cph = 9.81546301 # +- 0.000000001 (from DTU 2005) lat_cph = 55.67 * pi / 180.0 g_poles = 9.832 g_45 = 9.806 g_equator = 9.780 g_cph_model1 = g_45 + (g_equator - g_poles) * cos(2.0 * lat_cph) g_cph_model2 = 9.780327 * (1.0 + 0.0053024 * sqr(sin(lat_cph)) + 0.0000058*sqr(sin(2.0*lat_cph))) print(" The offical value for Copenhagen is: {0:7.5f} (compared to models: {1:7.5f} and {2:7.5f})".format(g_cph, g_cph_model1, g_cph_model2)) print("") print(" Distance from result to official value: {0:6.4f} +- {1:6.4f} ({2:3.1f} sigma)".format(g - g_cph, err_g, (g - g_cph) / err_g)) if (Blinding) : print(" NOTE: result is blinded!!!") # Hold script until you want to exit raw_input( ' ... Press enter to exit ... ' )