#coding: utf8

from __future__ import print_function, division, unicode_literals

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

##
## Parte 1
##

# misure
M  = sample(0.5096,0.5100,0.5090,0.5095,0.5096,0.5102,0.5104,0.5103) # massa (g)
D  = sample(4.45,4.46,4.47,4.47,4.47,4.47,4.47,4.46)                 # diametro (mm)
T0 = sample(2.43,2.28,2.46,2.42,2.46,2.44,2.32,2.29)      # time at 20cm, 22°C
T1 = sample(0.39,0.46,0.42,0.41,0.41,0.43,0.37,0.44,0.39) # time at 5cm,  22°C
T  = sample(2.13,2.29,2.44,2.25,2.29,2.19,2.19,2.07,2.25, # time at 20cm, 24°C
            2.14,2.29,2.27,2.25,2.10,2.11,2.12,2.26,2.03,
            2.10,2.18,2.28,2.29,2.13,2.14,2.17,2.24,2.33,
            2.04,2.06,2.12,2.17,1.97,1.97,2.14,2.07,2.22,
            2.05,1.87,2.14,2.21,1.89,2.32,2.23,2.14,2.07,
            2.04,2.00,2.14,1.96,2.40,2.22,2.04,2.16,2.01,
            2.08,2.03,2.04,2.09,2.30,2.09,2.14,1.92,2.01,
            2.27,2.26,2.19,2.19,2.07,1.94,2.00,1.87,2.20,
            2.02,1.95,2.03,2.27,2.14,1.92,1.94,1.94,2.07,
            1.82,1.88,1.86,2.13,2.05,2.20,1.92,2.09,2.34,
            2.36,2.24,2.24,2.20,2.07,2.35,2.20,2.13,2.10,2.12)

# compare velocity 5cm - 20cm
V0 = 0.200/T0.to_ufloat()
V1 = 0.050/T1.to_ufloat()
alpha = check_measures(V0, V1)

print('''
v_20cm: ({})m/s
v_5cm: ({})m/s
Χ²: α={:.5f}, α>ε: {}
'''.format(V0, V1, alpha, alpha>epsilon))

# times histogram
plt.figure(1)
plt.xlabel('tempo di caduta da 20cm (s)')
plt.ylabel('frequenza\'')
plt.xlim(1.80, 2.46)

# gaussian best fit
t = np.linspace(1.80, 2.45, 100)
gauss = mlab.normpdf(t, T.mean, T.std)

plt.plot(t, gauss, 'black', linewidth=1)
O, B, _ = plt.hist(T, bins=6, normed=1, facecolor='#506d72')
plt.show()

# χ² test histogram
alpha = chi_squared(T, O, B)
print('χ²: α={:.3f}, α>ε: {}'.format(alpha, alpha>epsilon))

# velocity
V = 0.200/T.to_ufloat()

# results
print("""
m: ({})g
d: ({})mm
T: ({})s
v: ({})m/s
""".format(M, D, T, V))


## 
## Part 2
##

# masses
M = [
      sample(0.0323,0.0326,0.0325,0.0327,0.0328)                      #2mm
    , sample(0.1102,0.1101,0.1103,0.1102,0.1100)                      #3mm
    , sample(0.2611,0.2610,0.2610,0.2607,0.2612,0.2610,0.2610)        #4mm
    , sample(0.5096,0.5099,0.5101,0.5096,0.5089,0.5099,0.5093,0.5093) #5mm
    , sample(0.8808,0.8802,0.8802,0.8795)                             #6mm
    ]

# diameters
D = [
      sample(1.98,1.99,2.00,2.00,1.99)           #2mm
    , sample(2.99,3.00,3.00,2.99,3.00)           #3mm
    , sample(3.99,3.99,3.99,4.00,4.00,4.00,3.99) #4mm
    , sample(4.99,5.00,4.99,5.00,4.99,5.00,4.99) #5mm
    , sample(6.00,6.00,5.99,5.99)                #6mm
    ]

# time to fall (20cm)
T = [
      sample(16.69,17.07,16.83,16.74,16.02,16.83,17.07,16.40,16.77,16.43,
             16.82,16.93,16.98,16.81,16.54,16.57,16.54,16.41,16.48,16.62) #2mm
    , sample(7.87,7.70,7.85,7.79,7.76,7.93,7.90,7.72,8.09,7.73,
             7.65,7.86,7.96,7.99,8.04,7.76,7.86,7.75,7.91,7.86)           #3mm
    , sample(4.57,4.63,4.39,4.46,4.61,4.64,4.58,4.56,4.53,4.63,
             4.48,4.72,4.63,4.88,4.53,4.82,4.58,4.70,4.39,4.43)           #4mm
    , sample(3.25,3.31,3.43,3.30,3.35,3.31,3.30,3.35,3.37,3.26,
             3.31,3.18,3.19,3.27,3.21,3.48,3.22,3.13,3.25,3.29)           #5mm
    , sample(2.68,2.60,2.58,2.59,2.53,2.50,2.60,2.53,2.38,2.63,
             2.43,2.57,2.48,2.41,2.53,2.55,2.42,2.46,2.43,2.56,2.50)      #6mm
    ]

m = 0.001*ufloat(M[4].mean, M[4].std) # sample for density calculation
d = 0.001*ufloat(D[4].mean, D[4].std)

SV  = sample(0.200/i.mean*i.std for i in T) #uncertainty on V (≈all identical)
V   = sample(0.200/i.mean       for i in T) #limit velocity (m/s)
R   = sample(0.001*i.mean/2     for i in D) #radii (m)
Rmm = sample(i.mean/2           for i in D) #radii (mm)


k = simple_linear(Rmm**2, V, SV)
f = lambda x: k.n*x
g = lambda x: k.n*x**2

# plot R - V
plt.figure(2)
plt.xlabel('raggio (mm)')
plt.ylabel('velocita\' terminale (m/s)')
plt.xlim(0, 4)
plt.ylim(0, 0.10)

x = np.linspace(0, Rmm.max**2, 100)
plt.plot(x, map(g, x), color='#973838')
plt.scatter(Rmm, V, color='#92182b')
plt.show()

# plot R² - V
plt.figure(3)
plt.xlabel('raggio (mm)^2')
plt.ylabel('velocita\' terminale (m/s)')
plt.xlim(0, 10)
plt.ylim(0, 0.10)

plt.plot(x, map(f, x), color='#286899')
plt.scatter(Rmm**2, V)
plt.errorbar(Rmm**2, V, yerr=SV, color='#3a3838', linestyle='')
plt.show()

#linear fit R² - V
k = simple_linear(R**2, V, SV)
f = lambda x: k.n*x

# χ² test for linear fit
alpha = chi_squared_fit(R**2, V, f, SV)
print('''
χ²: α={}
k: ({})1/m*s
'''.format(alpha, k))

# calculate eta
rho_s = m/(4/3* np.pi * (d/2)**3)
rho_g = ufloat(1261, 0)
g     = ufloat(9.81, 0)
eta   = 2/9*g*(rho_s-rho_g)/k

print('''
ρ: ({})kg/m³
η: ({})Pa*s
'''.format(rho_s, eta))