Speed of Light Part 2

Back to 2014.

Here is an example program for fitting the data from the speed of light experiment. It reads in an associated file media:averages.txt that contains the average values for the rotational velocity of the mirror and the displacement of the reflected beams focus position, as measured with the microscope. The error on the displacements is assumed to be common among the data points and is stored in the variable disp_var (displacement variance). The error in the rotational velocities is taken to be small, but this assumption is also checked in the code. If everything goes well, the program will output best fit parameters etc as well as a PNG plot of the data and fit: media:SL_data_with_fit.png.

import numpy as np
import pylab as plt
from scipy.stats.stats import chisqprob

vel_rev, disp = np.loadtxt('averages.txt')
print 'velocities', vel_rev
print 'displacements', disp
vel_rad = vel_rev*2*np.pi # Convert from rev/sec to radians/sec
disp_var = 0.00174356**2 #mm^2

#Plot data
fs = 16
plt.errorbar(vel_rev, disp, yerr=disp_var**.5, fmt='b.')
plt.xlabel("Rotational Velocity (rev/s)", fontsize=fs)
plt.ylabel("Displacement (mm)", fontsize=fs)
#plt.show()

#Compute linear least squares solution
def linear_fit(x,y,y_var): #function definition!
    """Assumes same variance on all y"""
    N=len(x)
    det = N*(x**2).sum() - (x.sum())**2 # determinant
    a = ((x**2).sum()*y.sum()-x.sum()*(x*y).sum())/det
    b = (N*(x*y).sum()-x.sum()*y.sum())/det
    a_var = (y**2).sum()*y_var/det
    b_var = N*y_var/det
    return a, b, a_var, b_var

a, b, a_var, b_var = linear_fit(vel_rad, disp, disp_var)

print "a =", a, '+/-', a_var**.5, 'mm'
print "b =", b, '+/-', b_var**.5, 'mm/(rad/s)'

#Plot line
plt.plot(vel_rev, a+b*vel_rad, 'blue')
#plt.show()

#Compute Chi-sq
nu = len(disp)-2 # number of deg of freedom
chisq = np.sum((disp-(a+b*vel_rad))**2/disp_var)
print "Chi-sq =", chisq, 'for', nu, 'degrees of freedom'
print "PTE = ", chisqprob(chisq,nu)

# PTE is a bit small...suspect error bars underestimated
# Compute correction to errors
disp_var_new = disp_var * chisq / nu
a, b, a_var, b_var = linear_fit(vel_rad, disp, disp_var_new)

print "Adjusted quantities"
print "New displacement error", disp_var_new**.5, 'mm'
print "a =", a, '+/-', a_var**.5, 'mm'
print "b =", b, '+/-', b_var**.5, 'mm/(rad/s)'

#Replot
plt.clf()
plt.errorbar(vel_rev, disp, yerr=disp_var_new**.5, fmt='b.')
plt.xlabel("Rotational Velocity (rev/s)", fontsize=fs)
plt.ylabel("Displacement (mm)", fontsize=fs)
plt.plot(vel_rev, a+b*vel_rad, 'blue')
#plt.show()
plt.savefig('data_with_fit.png')

#Check that velocity errors are safe to neglect
x_err = 6*np.pi #estimated error in rad/s
y_err_indirect = b*x_err
print 'displacement error due to velocity error', y_err_indirect, 'mm'

#Estimate Speed of Light
L1=930.
L2=613.
fL1=48.
MR=132. #error 5.
D=4870.  #error 20
A=L1-L2-fL1
B=L2-MR
print 'A =', A, 'mm'
print 'B =', B, 'mm'

c= 4*A*D**2/b/(D+B)
print "c = %.5f x 10^8 m/s" % (c/1e11)

#error
dc_b = c/b*b_var**.5 
print "sigma_c due to b error = %.5f x 10^8 m/s" % (dc_b/1e11)

Expected Output:

velocities [  100.5   351.    600.5   849.   1099.5   -99.5  -349.5  -600.   -850.
 -1100. ]
displacements [ 11.925   11.9535  11.981   12.0025  12.03    11.9075  11.8815  11.8575
  11.838   11.8095]
a = 11.9185850459 +/- 1.50210913346e-06 mm
b = 1.58667545642e-05 +/- 1.26028529539e-07 mm/(rad/s)
Chi-sq = 16.0709045274 for 8 degrees of freedom
PTE =  0.0413764210529
Adjusted quantities
New displacement error 0.00247122372129 mm
a = 11.9185850459 +/- 2.1290048651e-06 mm
b = 1.58667545642e-05 +/- 1.78625738005e-07 mm/(rad/s)
displacement error due to velocity error 0.000299081277452 mm
A = 269.0 mm
B = 481.0 mm
c = 3.00571 x 10^8 m/s
sigma_c due to b error = 0.03384 x 10^8 m/s