rnhmjoj/Lab-I
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

initial commit

Browse Source
This commit is contained in:
rnhmjoj 2016-06-19 20:37:57 +02:00
commit b5b8c1ffee
No known key found for this signature in database
GPG Key ID: 362BB82B7E496B7C
11 changed files with 1384 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

99
calorimetro.py Normal file
View File

@ -0,0 +1,99 @@
# coding: utf-8
from __future__ import print_function
from lab import *
import matplotlib.pyplot as plt
import numpy as np
##
## Part 1
## equivalent mass
mt = 353.36
mi = np.array([592.67, 576.99, 512.66, 479.75, 498.35])
mf = np.array([682.53, 678.08, 634.87, 626.51, 626.81])
T1 = np.array([20.2, 35.1, 21.1, 38.9, 22.5])
T2 = np.array([83.1, 87.0, 89.9, 87.2, 95.0])
Te = np.array([38.0, 51.0, 48.9, 61.9, 53.9])
m1 = mi-mt
m2 = mf-mi
me = m2*(T2-Te)/(Te-T1) - m1
me = me[2:].mean()
print('m_e: ', me)
##
## Part 2
## specific heat
mt = 359.67
mi = np.array([713.98, 717.14, 693.60])
ms = np.array([130.02, 121.85, 38.95])
m1 = mi-mt
Ts = 100
T1 = np.array([21.2, 23.2, 21.2])
Te = np.array([23.6, 25.3, 23.0])
cs = (Te-T1)*(m1+me)/((Ts-Te)*ms)
print('''
rame: {:.1f}
ottone: {:.1f}
alluminio: {:.1f}
'''.format(*cs*4186))
##
## Part 3
## Joule constant
mt = 354.71
m = 760.53 - mt
V = 14
I = 3
T = sample(22.0, 23.4, 24.8, 26.4, 27.6, 29.0, 30.4)
t = sample(0, 1, 2, 3, 4, 5, 6)*60
# linear interpolation
x = np.linspace(0,370)
a,b = linear(t, T, 0.1)
f = lambda x: a.n+b.n*x
# plot T - t
plt.xlabel('tempo (s)')
plt.ylabel('temperatura (s)')
plt.xlim(0,400)
plt.scatter(t, T, color='#135964')
plt.plot(t, f(t), '#317883')
plt.show()
# χ² test
alpha = chi_squared_fit(t, T, f, 0.1)
print('χ²: α={:.3f}, α>ε: {}'.format(alpha, alpha>epsilon))
# find J from b
J = I*V/(b*(m+me))
print('J={}'.format(J))
##
## Part 4
## latent heat
mt = 355.12 # calorimeter tare
mi = 386.61 # calorimeter + ice
mf = 561.62 # calorimeter + ice + water
m1 = mf-mi # water
m2 = mi-mt # ice
t1 = 91.1
t2 = -17
te = 66.8
ca = 4.2045 # water specific heat at 91°C
cg = 2.0461 # ice specific heat at -17°C
l = ca*(m1+me)/m2*(t1-te) + cg*t2 - ca*te
print(l)

145
corda.py Normal file
View File

@ -0,0 +1,145 @@
# coding: utf-8
from __future__ import division, unicode_literals
from lab import *
import numpy as np
import matplotlib.pyplot as plt
import uncertainties.umath as u
mass = lambda x: ufloat(x, 0.01)*0.001
length = lambda x: ufloat(x, 0.5)*0.01
##
## Costants
##
mL = [ (16.47, 251.0) # violet
, (5.49, 252.7) # black
, (6.16, 287.7) # green
, (1.48, 188.8) # white
]
mu = [mass(m)/length(l) for m,l in mL] # linear density (kg/m)
pp = mass(50.03) # tare mass (kg)
g = 9.80665 # acceleration due to gravity (m/s^2)
sigma = 0.8 # frequency uncertainty
## green string
## T,L,μ constant, plot ν(n)
m = mass(400.21) + pp
L = length(98.5)
T = m*g
n = np.arange(1,8)
nu = np.array([23.97, 48.20, 72.54, 96.97, 122.20, 147.10, 171.30])
# linear regression y=kx where
# y=ν, x=n, k=1/(2L)√(T/μ)
x = np.linspace(0,8,100)
k = simple_linear(n, nu, sigma)
f = lambda x: k.n*x
plt.figure(1)
plt.xlim(0,8)
plt.ylim(0,200)
plt.xlabel('n armonica')
plt.ylabel('frequenza di risonanza (Hz)')
plt.scatter(n, nu, color='#468174')
plt.plot(x, f(x), '#5ca897')
plt.show()
# χ² test
alpha = chi_squared_fit(n, nu, f, sigma, s=1)
v0 = 2*L*k # actual propagation velocity
v = u.sqrt(T/mu[2]) # theoretical propagation velocity
print '''
χ²: α={:.3f}, α>ε: {}
k: ({})1/s
v°: ({})m/s
v: ({})m/s
'''.format(alpha, alpha>epsilon, k, v0, v)
##
## green string
## T,μ,n constant, plot ν(T), n=2
M = [mass(i)+pp for i in [700.45, 499.95, 450.16, 400.19, 200.16, 140.92]]
T = np.array([g*i.n for i in M])
nu = np.array([60.26, 52.20, 49.70, 47.06, 34.96, 31.22])
# linear regression y=kx where
# y=ν, x=√T, k=1/(L√μ)
x = np.linspace(1, 8, 100)
k = simple_linear(np.sqrt(T), nu, sigma)
f = lambda x: k.n*np.sqrt(x)
plt.figure(2)
plt.xlabel('Tensione (N)')
plt.ylabel('Frequenza 2ª armonica (Hz)')
plt.scatter(T, nu, color='#673f73')
plt.plot(x, f(x), '#975ca8')
plt.show()
# χ² test
alpha = chi_squared_fit(T, nu, f, sigma, s=1)
print "χ²: α={}, α>ε: {}".format(alpha, alpha>epsilon)
##
## green string
## T,μ,n constant, plot ν(L), n=2
m = mass(400.21) + pp
T = m*g
L = np.array([116.5, 104.5, 95.4, 82.2, 75.0, 64.5, 54.7])*0.01
nu = np.array([44.58, 49.43, 54.60, 63.58, 69.34, 80.75, 96.40])
# linear regression y=kx where
# y=ν, x=1/L, k=n/2 √(T/μ)
x = np.linspace(0.5, 1.25, 100)
k = simple_linear(1/L, nu, sigma)
f = lambda x: k.n/x
plt.figure(3)
plt.xlabel('Lunghezza corda (m)')
plt.ylabel('Frequenza 2ª armonica (Hz)')
plt.scatter(L, nu, color='#844855')
plt.plot(x, f(x), '#a85c6c')
plt.show()
# χ² test
alpha = chi_squared_fit(L, nu, f, sigma, s=1)
print "χ^2: α={}, α>ε: {}".format(alpha, alpha>epsilon)
##
## every string
## T,n,L constant, plot ν(μ), n=2
m = mass(400.21) + pp
L = length(98.5)
n = 2
mu = np.array([i.n for i in mu])
nu = np.array([26.90, 47.23, 47.60, 78.20])
# linear regression y=kx where
# y=ν, x=1/√μ, k=√T/L
x = np.linspace(0.0005,0.007,100)
k = simple_linear(1/np.sqrt(mu), nu, sigma)
f = lambda x: k.n/np.sqrt(x)
plt.figure(4)
plt.xlabel('Densità lineare (kg/m)')
plt.ylabel('Frequenza 2ª armonica')
plt.scatter(mu, nu, color='#44507c')
plt.plot(x, f(x), '#5c6ca8')
plt.show()
# χ² test
alpha = chi_squared_fit(mu, nu, f, sigma, s=1)
print "χ²: α={}, α>ε: {}".format(alpha, alpha>epsilon)

137
kater.py Normal file
View File

@ -0,0 +1,137 @@
# coding: utf-8
from __future__ import division, print_function, unicode_literals
from lab import *
import numpy as np
import matplotlib.pyplot as plt
import uncertainties.umath as u
# known constants
D = ufloat(0.994, 0.001)
Ma = ufloat(1.000, 0.001)
Mb = ufloat(1.400, 0.001)
##
## Part 1
##
X = sample(0.067,0.080,0.100,0.135,0.150,0.170,0.196,0.227) # distance from 1
T1 = [ sample(2.0899,2.0842,2.0828,2.0823) # period 1
, sample(2.0604,2.0575,2.0574,2.0555)
, sample(2.0203,2.0186,2.0162,2.0170)
, sample(1.9509,1.9493,1.9485,1.9470)
, sample(1.9283,1.9260,1.9253,1.9208)
, sample(1.9038,1.9009,1.8994,1.9007)
, sample(1.8748,1.8731,1.8722,1.8733)
, sample(1.8463,1.8485,1.8458,1.8428)
]
T2 = [ sample(2.0196,2.0177,2.0174,2.0179) # period 2
, sample(2.0165,2.0159,2.0155,2.0151)
, sample(2.0073,2.0064,2.0066,2.0070)
, sample(1.9892,1.9899,1.9901,1.9913)
, sample(1.9880,1.9872,1.9865,1.9865)
, sample(1.9821,1.9815,1.9808,1.9828)
, sample(1.9747,1.9751,1.9740,1.9730)
, sample(1.9693,1.9683,1.9678,1.9672)
]
t1 = [i.mean for i in T1]
t2 = [i.mean for i in T2]
s1 = [i.stdm for i in T1]
s2 = [i.stdm for i in T1]
# quadratic fit
a1,b1,c1 = polynomial(X, t1, 2, sigma_Y=s1)
a2,b2,c2 = polynomial(X, t2, 2, sigma_Y=s2)
f = lambda x: a1*x**2+b1*x+c1
g = lambda x: a2*x**2+b2*x+c2
x = np.linspace(X.min, X.max, 100)
# plot T1, T2 - X
plt.figure(1)
plt.xlabel('distanza da 1 (m)')
plt.ylabel('periodo (s)')
plt.xlim(0.05,0.24)
plt.scatter(X, t1, color='#21812b')
plt.scatter(X, t2, color='#7755ff')
plt.plot(x, [f(i).n for i in x], color='#21812b')
plt.plot(x, [g(i).n for i in x], color='#7755ff')
plt.show()
# find intersection and calculate g
x = ((b2-b1)-(u.sqrt)((b1-b2)**2-4*(a1-a2)*(c1-c2)))/(2*(a1-a2)) # intersection
T = f(x) # isochronic period
g1 = D/T**2*4*np.pi**2 # acceleration due to mass
print('''
Kater pendulum:
T: ({})s
g₁: ({})m/s²
'''.format(T, g1))
##
## Part 2
## free fall - direct measure of g
H = sample(0.300,0.400,0.500,0.550,0.650,0.750) # drop height
T = [ sample(0.241,0.242,0.241,0.248,0.251,0.252,0.249,0.243,0.248,0.244) # time to fall
, sample(0.285,0.285,0.284,0.284,0.284,0.282,0.284,0.282,0.281,0.281)
, sample(0.318,0.319,0.318,0.316,0.318,0.316,0.316,0.318,0.317,0.317)
, sample(0.336,0.333,0.333,0.332,0.340,0.330,0.330,0.332,0.335,0.330)
, sample(0.362,0.361,0.365,0.367,0.359,0.361,0.357,0.366,0.361,0.359)
, sample(0.385,0.387,0.390,0.385,0.393,0.386,0.386,0.392,0.394,0.387)
]
t = sample(i.mean for i in T)
St = sample(i.std for i in T)
# linear fit y=kx where y=h, k=g/2, x=t^2
x = np.linspace(0.2, 0.4, 100)
k = simple_linear(t**2, H, 0.02)
f = lambda x: x**2 * k.n
# plot H - T
plt.figure(2)
plt.ylabel('altezza (m)')
plt.xlabel('tempo di caduta (s)')
plt.scatter(t, H)
plt.plot(x, f(x))
plt.show()
# calculate g
g2 = 2*k
# χ² test
alpha = chi_squared_fit(t, H, f, St)
print('''
Free fall:
g₂: ({})m/s²
χ²: α={:.3f}, α>ε: {}
'''.format(g2, alpha, alpha>epsilon))
##
## Part 3
## compare g1, g2
alpha = check_measures(g1, g2)
g = combine_measures(g1, g2)
print('''
compare g₁ g₂:
α: {:.3f}, α>ε: {}
combine:
g_best: ({})m/s²
'''.format(alpha, alpha>epsilon, g))

233
lab.py Normal file
View File

@ -0,0 +1,233 @@
# coding: utf-8
from __future__ import division, unicode_literals
from sympy.functions.special.gamma_functions import lowergamma, gamma
from sympy.functions.special.error_functions import erf
from uncertainties.core import UFloat
from uncertainties import ufloat, wrap, correlated_values
from numpy.polynomial.polynomial import polyfit
import collections
import types
import numpy as np
import string
##
## Variables
##
class sample(np.ndarray):
"""
Sample type (ndarray subclass)
Given a data sample outputs many statistical properties:
n: number of measures
min: minimum value
max: maximum value
mean: sample mean
med: median value
var: sample variance
std: sample standard deviation
stdm: standard deviation of the mean
"""
__array_priority__ = 2
def __new__(cls, *sample):
if isinstance(sample[0], types.GeneratorType):
sample = list(sample[0])
s = np.asarray(sample).view(cls)
return s
def __str__(self):
return format_measure(self.mean, self.stdm)
def __array_finalize__(self, obj):
if obj is None:
return
obj = np.asarray(obj)
self.n = len(obj)
self.min = np.min(obj)
self.max = np.max(obj)
self.mean = np.mean(obj)
self.med = np.median(obj)
self.var = np.var(obj, ddof=1)
self.std = np.sqrt(self.var)
self.stdm = self.std / np.sqrt(self.n)
def to_ufloat(self):
return ufloat(self.mean, self.stdm)
##
## Normal distribution
##
def phi(x, mu, sigma):
"""
Normal CDF
computes the probability P(x<a) given X is a normally distributed
random variable with mean μ and standard deviation σ.
"""
return float(1 + erf((x-mu) / (sigma*np.sqrt(2))))/2
def p(mu, sigma, a, b=None):
"""
Normal CDF shorthand
given a normally distributed random variable X, with mean μ
and standard deviation σ, computes the probability P(a<X<b)
or P(X<a) whether b is given.
"""
if b:
return phi(b, mu, sigma) - phi(a, mu, sigma) # P(a<X<b)
return phi(a, mu, sigma) # P(x<a)
def chi(x, d):
"""
χ² CDF
given a χ² distribution with degrees of freedom
computes the probability P(X^2x)
"""
return float(lowergamma(d/2, x/2)/gamma(d/2))
def q(x, d):
"""
χ² 1-CDF
given a χ² distribution with d degrees of freedom
computes the probability P(X^2>x)
"""
return 1 - chi(x, d)
def chi_squared(X, O, B, s=2):
"""
χ² test for a histogram
given
X: a normally distributed variable X
O: ndarray of normalized absolute frequencies
B: ndarray of bins delimiters
s: number of constraints on X
computes the probability P(χ²>χ₀²)
"""
N, M = X.n, len(O)
d = M-1-s
delta = (X.max - X.min)/M
E = [N*p(X.mean, X.std, B[k], B[k+1]) for k in range(M)]
Chi = sum((O[k]/delta - E[k])**2/E[k] for k in range(M))
return q(Chi, d)
def chi_squared_fit(X, Y, f, sigma=None, s=2):
"""
χ² test for fitted data
given
X: independent variable
Y: dependent variable
f: best fit function
s: number of constraints on the data (optional)
sigma: uncertainty on Y, number or ndarray (optional)
"""
if sigma is None:
sigma_Y = np.mean([i.std for i in Y])
else:
sigma_Y = np.asarray(sigma).mean()
Chi = sum(((Y - f(X))/sigma_Y)**2)
return q(Chi, Y.size-s)
##
## Regressions
##
def simple_linear(X, Y, sigma_Y):
"""
Linear regression of line Y=kX
"""
sigma = np.asarray(sigma_Y).mean()
k = sum(X*Y) / sum(X**2)
sigma_k = sigma / np.sqrt(sum(X**2))
return ufloat(k, sigma_k)
def linear(X, Y, sigma_Y):
"""
Linear regression of line Y=A+BX
"""
sigma = np.asarray(sigma_Y).mean()
N = len(Y)
D = N*sum(X**2) - sum(X)**2
A = (sum(X**2)*sum(Y) - sum(X)*sum(X*Y))/D
B = (N*sum(X*Y) - sum(X)*sum(Y))/D
sigma_A = sigma * np.sqrt(sum(X**2)/D)
sigma_B = sigma * np.sqrt(N/D)
return (ufloat(A, sigma_A),
ufloat(B, sigma_B))
def polynomial(X, Y, d=2, sigma_Y=None):
"""
d-th degree polynomial fit
"""
weights = None
if sigma_Y:
if isinstance(sigma_Y, collections.Iterable):
weights = 1/np.asarray(sigma_Y)
else:
weights = np.repeat(1/sigma_Y, len(X))
coeff, cov = np.polyfit(X, Y, d, w=weights,
cov=True, full=False)
print cov
return correlated_values(coeff, cov)
##
## Misc
##
def check_measures(X1, X2):
"""
Checks whether the results of two set of measures
are compatible with each other. Gives the α-value
α=1-P(X within tσ) where t is the weighted difference
of X₁ and X₂
"""
t = np.abs(X1.n - X2.n)/np.sqrt(X1.s**2 + X2.s**2)
print t
return float(1 - erf(t/np.sqrt(2)))
def combine_measures(*variables):
"""
Combines different (compatible) measures of the same
quantity producing a new value with a smaller uncertainty
than the individually taken ones
"""
W = np.array([1/i.s**2 for i in variables])
X = np.array([i.n for i in variables])
best = sum(W*X)/sum(W)
sigma = 1/np.sqrt(sum(W))
return ufloat(best, sigma)
def format_measure(mean, sigma):
"""
Formats a measure in the standard declaration
"""
prec = int(np.log10(abs(sigma)))
limit = 2-prec if prec < 0 else 2
sigma = round(sigma, limit)
mean = round(mean, limit)
return "{}±{}".format(mean, sigma)
sin = wrap(lambda x: np.sin(np.radians(x)))
epsilon = 0.05

100
oscillatore.py Normal file
View File

@ -0,0 +1,100 @@
# coding: utf-8
from __future__ import unicode_literals, division
from lab import *
import uncertainties.umath as u
import matplotlib.pyplot as plt
import numpy as np
##
## almost free oscillation
## θ(t) = Ae^(-γt)sin(ωt+φ)+θ0
G = sample(0.0421, 0.0514, 0.0421, 0.0392)
W = sample(2.6800, 2.6600, 2.6709, 2.6700)
g = ufloat(G.mean, G.std) # damping ratio γ
w = ufloat(W.mean, W.std) # pulsation ω
w0 = u.sqrt(w**2-g**2) # corrected value ω0
T = 2*np.pi/w # oscillation period
print('ω0: {}, γ: {}, T: {}'.format(w0, g, T))
##
## damped oscillation
##
G = sample(0.2280, 0.2130, 0.2180, 0.2150)
W = sample(2.6602, 2.6500, 2.6600, 2.6600)
g = ufloat(G.mean, G.std) # damping ratio γ
w = ufloat(W.mean, W.std) # pulsation ω
w0 = u.sqrt(w**2-g**2) # corrected value ω0
T = 2*np.pi/w # oscillation period
print('ω0: {}, γ: {}, T: {}'.format(w0, g, T))
##
## even more damped
##
G = sample(0.6060, 0.6150, 0.6057, 0.6060)
W = sample(2.6000, 2.6200, 2.6100, 2.5900)
g = ufloat(G.mean, G.std) # damping ratio γ
w = ufloat(W.mean, W.std) # pulsation ω
w0 = u.sqrt(w**2-g**2) # corrected value ω0
T = 2*np.pi/w # oscillation period
print('ω0: {}, γ: {}, T: {}'.format(w0, g, T))
##
## Driven oscillation
##
V = np.array([1.8000, 1.9000, 2.0000, 2.1000, 2.2000, # generator tension (V)
2.3000, 2.3500, 3.0000, 3.2000, 3.4000])
A = np.array([2.3900, 2.5400, 4.0000, 3.2900, 5.4500, # oscillation amplitude (rad)
6.2500, 4.5200, 3.7500, 2.4200, 1.8200])
T = np.array([4.3100, 4.0600, 3.1000, 3.4000, 2.8300, # oscillation period (s)
2.7500, 2.9800, 1.9700, 1.8300, 1.7000])
Tf = np.array([4.3300, 4.1048, 3.0883, 3.4069, 2.8270, # rotor period (s)
2.7229, 2.9592, 1.9721, 1.8240, 1.7210])
w = 2*np.pi/T # disc angular velocity
wf = 2*np.pi/Tf # rotor angular velocity
# polynomial fit Y=aX^2+bX+c
Y = 1/A**2
X = w**2
y = 1/ufloat(A.mean(), 0.08)**2 # error estimation
a,b,c = polynomial(X, Y, 2, y.s) # polynomial coeffs
# plot Α - ω
x = np.linspace(1.4,3.8,250)
f = lambda w: 1/np.sqrt(a.n*w**4 + b.n*w**2 + c.n)
plt.xlabel('pulsazione (rad/s)')
plt.ylabel('ampiezza (rad)')
plt.scatter(w, A, color='#b99100')
plt.plot(x, f(x), '#c89d00')
plt.show()
# find curve parameters
M = 1/u.sqrt(a)
w0 = (M**2 * c)**(1/4)
g = 1/2*u.sqrt(b*M**2 + 2*w0**2)
print '''
M: {}
ω0: {}
γ: {}
'''.format(M,w0,g)

144
ottica.py Normal file
View File

@ -0,0 +1,144 @@
# 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_)))

23
shell.nix Normal file
View File

@ -0,0 +1,23 @@
{ pkgs ? import <nixpkgs> {}, mode ? "shell" }:
with pkgs;
let
matplotlib = pythonPackages.matplotlib.override {
enableGtk3 = true;
};
in stdenv.mkDerivation rec {
name = "lab-data";
source = ".";
buildInputs = with pythonPackages; [
matplotlib numpy
sympy uncertainties
];
shellHook = ''
#exec 2> /dev/null
exec fish
'';
}

154
termocoppia.pickle Normal file
View File

File diff suppressed because one or more lines are too long

68
termocoppia.py Normal file
View File

@ -0,0 +1,68 @@
# coding: utf-8
from lab import *
import pickle
import numpy as np
import matplotlib.pyplot as plt
def plot((name, points), n):
plt.figure(n+1)
if name == 'azoto':
plt.ylim(-5,-5.5)
elif name == 'acqua':
plt.ylim(4, 4.5)
elif name =='stearico':
name = 'acido stearico'
elif name == 'glicole':
plt.ylim(-0.7,-0.4)
plt.title(name.capitalize())
plt.xlabel('tempo')
plt.ylabel('differenza di potenziale (mV)')
plt.plot(filter(lambda x: x!=0, points), '#44b34e')
plt.show()
set = pickle.load(open('termocoppia.pickle'))
for i, curve in enumerate(set.iteritems()):
#plot(curve, i)
pass
V = map(np.mean,
[ set['stagno'][150:250] # Tf
, set['indio'][125:250] # Tf
, set['acqua'][150:225] # Te
, set['stearico'][80:120] # Tf
, [0] # miscela acqua-ghiaccio
, set['glicole'][130:170] # Tf
, set['etere'][50:150] # Tf
, set['azoto'] # Te
])
T = [ 505.08 # Tf Sn
, 429.75 # Tf In
, 373.15 # Te H2O
, 342.40 # Tf CH3-(CH2)16-COOH
, 273.15 # Tf H2O
, 260.20 # Tf HO-(CH2)2-OH
, 156.80 # Tf CH3CH2-O-CH2CH3
, 77.36 # Te N2
]
x = np.linspace(-6, 6, 100)
a,b,c,d,e = polynomial(V, T, 4)
f = lambda x: a.mean*x**4 + b.mean*x**3 + c.mean*x**2 + d.mean*x + e.mean
g = lambda n,x: sum(a*x**i for i, a in enumerate(polyfit(V,T,n)))
# q(np.sum((np.polyval(np.polyfit(V,T,7), V) - T) ** 2), 1)*100
plt.figure(8)
plt.xlabel('Differenza di potenziale (mV)')
plt.ylabel('Temperatura (K)')
plt.scatter(V, T, color='#6972cc')
plt.plot(x, g(5,x), '#6972cc')
plt.show()
print '''
T(V)={a.n:.3f}V^4+{b.n:.3f}V^3{c.n:.3f}V^2+{d.n:.3f}V+{e.n:.3f}
'''.format(**dict(a=a,b=b,c=c,d=d,e=e))

115
urti.py Normal file
View File

@ -0,0 +1,115 @@
# coding: utf-8
from __future__ import division
from lab import *
import numpy as np
# helpers
length = lambda x: x/100
mass = lambda x: x/1000
# constants
m = mass(260.04) # cart mass (kg)
g = 9.80665 # acceleration due to gravity
##
## Collisions
##
# 3xmagnet, 2xrubber bumper, 2xputty
I = [0.23, 0.28, 0.36, 0.21, 0.27, 0.11, 0.15] # force impulse
# cart velocity before-after colliding
v = [ (0.46,-0.44), (0.55,-0.54), (0.70,-0.71) # elastic
, (0.52,-0.29), (0.68,-0.35) # partially anelastic
, (0.41, 0.00), (0.59, 0.00) # anelastic
]
for i,_ in enumerate(I):
K0, K1 = 1/2*m*v[i][0]**2, 1/2*m*v[i][1] # kinetic energy
p0, p1 = m*v[i][0], m*v[i][1] # momentuum
dp = round(p1-p0, 3)
dk = round(K1-K0, 3)
print '''
urto {}:
Δp={:.3f}, I={}, ΔK={:.2f} {:.2f}%
'''.format(i+1, dp, I[i], dk, abs((dp+I[i])/dp*100))
##
## Springs
##
x0 = 0.8675 # spring at rest offset
# force impulse
I = [ 0.32, 0.26, 0.35, 0.21, 0.34 # spring A
, 0.42, 0.47, 0.50, 0.45, 0.68 # spring B
]
# cart velocity before-after colliding
v = [ (0.61,-0.58), (0.51,-0.48), (0.66,-0.66), (0.41,-0.38), (0.67,-0.63) # spring A
, (0.81,-0.77), (0.89,-0.87), (0.94,-0.93), (0.85,-0.83), (1.32,-1.26) # spring B
]
# maximum position reached (spring compressed)
x = [ 0.888, 0.885, 0.890, 0.882, 0.887 # spring A
, 0.881, 0.884, 0.885, 0.885, 0.888 # spring B
]
# maximum force
F = [ 5.55, 4.58, 6.29, 3.66, 6.10 # spring A
, 18.61, 20.97, 22.25, 19.17, 31.25 # spring B
]
ka = []
kb = []
for i,_ in enumerate(I):
k = F[i]/(x[i]-x0) # spring elastic constant
Ep = 1/2*k*(x[i] -x0)**2 # theoretical value (assuming ΔEm=0)
K0, K1 = 1/2*m*v[i][0]**2, 1/2*m*v[i][1] # kinetic energy
p0, p1 = m*v[i][0], m*v[i][1] # momentuum
dp = round(p1-p0, 3)
dk = round(K1-K0, 3)
dEm = round(K0-Ep, 3)
if i<5:
ka.append(k)
else:
kb.append(k)
print '''
urto {}:
Δp={:.3f}, I={:.2f}, ΔK={:.2f} {:.2f}%
k_{}={:.3f}, ΔΕm={:.3f}
'''.format(i+1, dp, I[i], dk, abs((dp+I[i])/dp*100),('a' if i<5 else 'b'), k, dEm)
##
## Attrito
##
## cart
# a=0 => μ=tanθ=h/l
mu_s = length(11.5)/length(128.7)
# a=g(sinθ-μcosθ)
a = sample(1.54,1.54,1.50,1.48,1.54,1.53)
h, l = length(52.0), length(117.1)
mu_d = (h/l) - (a.mean/g)*np.sqrt(1+(h/l)**2)
print 'μs={:.2f} μd={:.2f}'.format(mu_s, mu_d)
## teflon block
# a=0 => μ=tanθ=h/l
mu_s = length(19.8)/length(127.0)
# a=g(sinθ-μcosθ)
a = sample(0.67,0.79,0.76,0.75,0.78,0.73)
h, l = length(34.1), length(125.5)
mu_d = (h/l) - (a.mean/g)*np.sqrt(1+(h/l)**2)
print 'μs={:.2f} μd={:.2f}'.format(mu_s, mu_d)

166
viscosita.py Normal file
View File

@ -0,0 +1,166 @@
#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))