Speed of Light Part 2

Here is an example program for fitting the data from the speed of light experiment. It reads in an associated file called "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.

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 = 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()

#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

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)