# coding: utf-8
from __future__ import print_function, division, unicode_literals

from lab import *
import numpy as np
import matplotlib.pyplot as plt

np.set_printoptions(precision=2)

epsilon = 0.05

##
## Parte 1
## test Snell Law and calculate n

i_a = sample(5,10,15,20,25,30,35) # direct path
r_a = sample(8,15,23,30,40,48,59)

i_b = sample(8,15,23,30,40,48,59) # reversed path
r_b = sample(6,10,16,20,26,30,36)

# linear regression y=kx, where y=sin(r^) x=sin(i^)
x = np.linspace(0, 1, 100)
ka = simple_linear(sin(i_a), sin(r_a), sin(1))
kb = simple_linear(sin(i_b), sin(r_b), sin(1))

# plot r^ - i^
plt.figure(1)
plt.xlim(0,1.0)
plt.ylim(0,1.0)
plt.xlabel('sin(angolo incidenza)')
plt.ylabel('sin(angolo rifrazione)')
plt.scatter(sin(i_a), sin(r_a), color='gray')
plt.plot(x, ka.n*x, '#72aa3d')
plt.scatter(sin(i_b), sin(r_b), color='gray')
plt.plot(x, kb.n*x, '#3d72aa')
plt.show()

# χ^2 test
alpha_a = chi_squared_fit(sin(i_a), sin(r_a), lambda x: ka.n*x, sin(1))
alpha_b = chi_squared_fit(sin(i_b), sin(r_b), lambda x: kb.n*x, sin(1))

print('''
χ^2 tests:
 α={:.3f}, α>ε: {}
 α_rev={:.3f}, α>ε: {}
'''.format(alpha_a, alpha_a>epsilon,
           alpha_b, alpha_b>epsilon))

#compare n, 1/n_rev
alpha = check_measures(ka, 1/kb)
n_best = combine_measures(ka, 1/kb)
print('''
n: ({})
n_rev: ({})
α: {:.2f}
n_best: ({})
'''.format(ka, kb, alpha, n_best))


##
## Parte 2
## calculate n from minimum deviation δ

delta_min = ufloat(38,1) # minimum deviation angle
alpha     = ufloat(60,0) # equilater triangle prism 
n = sin(delta_min/2+alpha/2)/sin(alpha/2)

print('''
n: ({})
n_rev: ({})
n_δ: ({})
'''.format(ka, kb, n))


##
## Part 3
## test thin lens equation, find focal lengths

## lens with f=200mm

d = 1.5 # transversal object length (cm)

s  = sample(110.0,100.0,90.0,80.0) # screen distance (p+q)
p  = sample(24.9,26.3,28.1,32.6)   # object - lens distance
q  = sample(84.7,73.7,61.5,46.8)   # lens - sceen dinstance
t  = sample(5.1,4.2,3.3,2.1)       # transversal image length
t_ = d*q/p                         # theoretical length   

# linear fit y=a+bx, where x=1/p, y=1/q
x = np.linspace(0.030, 0.041, 100)
a, b = linear(1/p, 1/q, 2e-4) # 1/f, -1
f = 1/a                       # focal length

# plot y=a+bx, where x=1/p, y=1/q
plt.figure(2)
plt.xlabel('1/p (cm)')
plt.ylabel('1/q (cm)')
plt.scatter(1/p, 1/q, color='#ba6b88')
plt.plot(x, a.n + b.n*x, '#ba6b88')
plt.show()

# χ^2 test
alpha = chi_squared_fit(1/p, 1/q, lambda x: a.n + b.n*x, 2e-4)
print('χ^2: α={:.3f}, α>ε: {}'.format(alpha, alpha>epsilon))

print('''
lente 200mm
f: ({})cm
t:  {}cm
t°: {}cm
'''.format(f, np.asarray(t), np.asarray(t_)))

## lens with f=100mm

s  = sample(80,70,60,50)         # see above
p  = sample(12.5,13.0,13.7,15.4) #
q  = sample(67.2,57.0,46.3,34.4) #
t  = sample(7.9,6.7,5.5,3.3)     #
t_ = d*q/p                       #

# linear fit y=a+bx
x    = np.linspace(0.064, 0.081, 100)
a, b = linear(1/p, 1/q, 2e-4) # 1/f, -1
f    = 1/a                    # focal length

# plot y=a+bx, where x=1/p, y=1/q
plt.figure(3)
plt.xlabel('1/p (cm)')
plt.ylabel('1/q (cm)')
plt.scatter(1/p, 1/q, color='#6b88ba')
plt.plot(x, a.n + b.n*x, '#6b88ba')
plt.show()

# χ^2 test
alpha = chi_squared_fit(1/p, 1/q, lambda x: a.n + b.n*x, 2e-4)
print('χ^2: α={:.3f}, α>ε: {}'.format(alpha, alpha>epsilon))

print('''
lente 100mm
f: ({})cm
t:  {}cm
t°: {}cm
'''.format(f, np.asarray(t), np.asarray(t_)))