# coding: utf-8
from __future__ import division, print_function, unicode_literals
from lab import *
from math import pi, sqrt
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab

##
## measure focal length
##

# convert micrometer readings (mm)
length = lambda x,y,z: x+0.5*y+0.01*z

# object - lens distance (mm)
p = ufloat(4608, 10, tag='p')

# object size (mm)
h0 = array(10.00, 15.00, 18.60, 20.00, 25.00, 30.00, 35.00, 40.00, 45.00, 50.00)

# image size (mm)
h1 = array(length(11, 0, 11.5) - length(10, 1,  1.5),
           length(11, 0, 39.5) - length(10, 1,  1.5),
           length(11, 0, 27.5) - length(10, 0, 20.5),
           length(11, 0, 41.0) - length(10, 0, 27.0),
           length(11, 0, 34.5) - length( 9, 1, 42.0),
           length(11, 0, 46.5) - length( 9, 1, 28.0),
           length(11, 1,  3.0) - length( 9, 1,  3.0),
           length(11, 1, 32.0) - length( 9, 1,  1.5),
           length(11, 0, 23.5) - length( 8, 1, 14.5),
           length(12, 0,  5.0) - length( 9, 0, 22.0))


# linear fit
sigma_h = 0.03
x  = np.linspace(5, 55, 10)
k  = simple_linear(h0, h1, sigma_h)
k.tag = 'k'
f  = ufloat(252, 1)
f0 = p*k/(k+1)

# tests
alpha = chi_squared_fit(h0, h1, lambda x: k.n*x, sigma_h, s=1)
beta  = check_measures(f0, f)

print(mformat('''
# focal lenth measure
k: {}
fₒ: {}
χ² test: α={:.2f}
t-test fₒ-f: α={:.2f}
''', k, f0, alpha, beta))

plt.figure(1)
plt.title('focal length (f=252mm)')
plt.xlabel('object size (mm)')
plt.ylabel('image size (mm)')
plt.scatter(h0, h1, color='#206a99')
plt.errorbar(h0, h1, sigma_h, color='#1a567c',linestyle='')
plt.plot(x, k.n*x, color='#1a567c')
plt.show()

f = f0*1e-3


##
## measure Δs
##

ds = sample(
  length(12,0,5.5) - length(11,1,21.5),
  length(12,0,4.5) - length(11,1,20.0),
  length(12,0,5.5) - length(11,1,21.0),
  length(12,0,4.5) - length(11,1,20.5),
  length(12,0,6.5) - length(11,1,22.5),
  length(12,0,6.5) - length(11,1,21.0),
  length(12,0,6.5) - length(11,1,21.0),
  length(12,0,3.5) - length(11,1,20.0),
  length(12,0,5.5) - length(11,1,21.0),
  length(12,0,5.5) - length(11,1,20.5),
  length(12,0,6.5) - length(11,1,22.0),
  length(12,0,4.5) - length(11,1,21.0),
  length(12,0,3.5) - length(11,1,20.0),
  length(12,0,5.5) - length(11,1,21.5),
  length(12,0,7.5) - length(11,1,22.5),
  length(12,0,4.5) - length(11,1,21.0),
  length(12,0,4.0) - length(11,1,21.0),
  length(12,0,6.5) - length(11,1,22.0),
  length(12,0,6.0) - length(11,1,21.0),
  length(12,0,6.5) - length(11,1,22.0),
  length(12,0,4.5) - length(11,1,20.0),
  length(12,0,6.5) - length(11,1,21.5),
  length(12,0,6.5) - length(11,1,21.0),
  length(12,0,5.5) - length(11,1,22.0),
  length(12,0,5.5) - length(11,1,20.5),
  length(12,0,5.5) - length(11,1,22.0),
  length(12,0,7.5) - length(11,1,22.5),
  length(12,0,4.5) - length(11,1,20.0),
  length(12,0,7.5) - length(11,1,22.5),
  length(12,0,5.5) - length(11,1,21.5),
  length(12,0,4.5) - length(11,1,21.0),
  length(12,0,5.5) - length(11,1,21.0),
  length(12,0,5.5) - length(11,1,20.0),
  length(12,0,5.5) - length(11,1,20.0),
  length(12,0,8.5) - length(11,1,22.5),
  length(12,0,6.5) - length(11,1,20.0),
  length(12,0,6.5) - length(11,1,22.0),
  length(12,0,5.5) - length(11,1,20.5),
  length(12,0,4.5) - length(11,1,20.5),
  length(12,0,5.5) - length(11,1,20.0),
  length(12,0,6.5) - length(11,1,22.5),
  length(12,0,4.5) - length(11,1,20.5),
  length(12,0,4.5) - length(11,1,22.0),
  length(12,0,7.5) - length(11,1,22.5),
  length(12,0,4.0) - length(11,1,20.0),
  length(12,0,4.5) - length(11,1,20.0),
  length(12,0,5.5) - length(11,1,22.0),
  length(12,0,5.5) - length(11,1,21.0),
  length(12,0,4.5) - length(11,1,20.5),
  length(12,0,7.5) - length(11,1,22.5))

## Δs histogram
plt.figure(2)
plt.xlabel('displacement (mm)')
plt.ylabel('frequency')

# gaussian best fit
x = np.linspace(ds.min-1e-3, ds.max+1e-3, 300)
gauss = mlab.normpdf(x, ds.mean, ds.std)

# χ² test
O, B, _ = plt.hist(np.asarray(ds), bins=6, normed=1, facecolor='#506d72')
alpha   = chi_squared(ds, O, B)
print('''
# Δs histogram
χ² test: α={:.2f}
'''.format(alpha))

# plot
plt.plot(x, gauss, 'black', linewidth=1)
plt.show()

ds = sample(*1e-3*ds).val()
ds.tag = 'Δs'


##
## calculate c
##

B = ufloat(0.292, 0.001, tag='B') # lens-rotating mirror distance
D = ufloat(5.491, 0.005, tag='D') # rotating-fixed mirrors distance
w = ufloat(3000,  2,     tag='ω') # motor angular velocity

A = (B+D)*f/(B+D-f)
c0 = 8*np.pi * A*D**2 * w / ((D+B)*ds) # speed of light

print(mformat('''
# speed of light
c: {}
σc/c: {:.2f}%

error components:''',c0, c0.s/c0.n*100))

for p, sigma in c0.error_components().items():
  print('  {}: {:.2f}%'.format(p.tag, abs(sigma/c0.n)*100))

# t-test cₒ-c
c = ufloat(299792458,0)
alpha = check_measures(c0, c)

print(mformat('t-test cₒ-c: α={:.2f}', alpha))
exit()