Lab-I/viscosita.py

167 lines
4.7 KiB
Python
Raw Normal View History

2016-06-19 20:37:57 +02:00
#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/
η: ({})Pa*s
'''.format(rho_s, eta))