315 lines
9.5 KiB
Python
315 lines
9.5 KiB
Python
# coding: utf-8
|
||
from __future__ import print_function, division, unicode_literals
|
||
|
||
import numpy as np
|
||
import uncertainties.umath as um
|
||
import matplotlib.pyplot as plt
|
||
from lab import *
|
||
|
||
##
|
||
## part II
|
||
##
|
||
|
||
R = ufloat(9.93e2, 0.05e2) # resistance
|
||
L = ufloat(3.82e-3, 0.5e-3) # inductance
|
||
C = ufloat(8.32e-6, 0.5e-6) # capacitance
|
||
Rl = ufloat(19.9, 0.1) # inductor resistance
|
||
|
||
## Resistor
|
||
|
||
# frequency (Hz)
|
||
nu = array( 50, 100, 150, 200, 300, 400, 500, 600, 800, 1e3,
|
||
1.5e3, 2e3, 3e3, 5e3, 8e3, 10e3, 15e3, 20e3, 30e3, 40e3,
|
||
60e3, 80e3, 100e3, 200e3, 250e3)
|
||
|
||
# input tension (V)
|
||
Va = array(10.2, 10.2, 10.2, 10.0, 9.8, 9.8, 9.8, 9.8, 9.8, 9.8, 8.6,
|
||
10.0, 9.8, 10.0, 9.6, 10.0, 10.0, 10.2, 10.4, 10.4, 10.4, 10.4,
|
||
10.6, 10.6, 10.6)/2
|
||
|
||
# output tension (V)
|
||
Vb = array(2.56, 4.80, 6.40, 7.20, 8.40, 8.80, 9.20, 9.40, 9.50, 9.60, 8.50,
|
||
9.80, 9.70, 9.75, 9.20, 9.40, 9.00, 8.60, 7.60, 6.80, 5.40, 4.40,
|
||
3.44, 1.80, 1.44)/2
|
||
|
||
# phase difference Va - Vb (s)
|
||
Oab = -array(4.0e-3, 1.7e-3, 840e-6, 580e-6, 300e-6, 160e-6, 110e-6, 60e-6,
|
||
40e-6, 20e-6, 4e-6, 2e-6, -4e-6, -4e-6, -3.6e-6, -4.4e-6,
|
||
-4e-6, -4.4e-6, -3.8e-6, -3.3e-6, -2.8e-6, -2.2e-6, -1.88e-6, -1.12e-6,
|
||
-0.92e-6)
|
||
|
||
# phase error (s)
|
||
sigmaOab = array(0.1e-3, 0.1e-3, 20e-6, 20e-6, 20e-6, 20e-6, 20e-6, 20e-6,
|
||
2e-6, 2e-6, 2e-6, 2e-6, 2e-6, 2e-6, 0.2e-6, 0.2e-6,
|
||
0.2e-6, 0.2e-6, 0.2e-6, 0.2e-6, 0.2e-6, 0.2e-6, 0.2e-6, 0.05e-6,
|
||
0.05e-6)
|
||
|
||
# transfer function
|
||
H = Vb/Va * np.exp(2j*np.pi*nu * Oab)
|
||
|
||
# estimate uncertainties
|
||
sigmaHm = (ufloat(Vb[0], 0.1)/ufloat(Va[0], 0.1)).s
|
||
sigmaHp = array(*((2*np.pi*v * ufloat(o, so)).s for v, o, so in zip(nu, Oab, sigmaOab)))
|
||
|
||
# magnitude, plot and fit
|
||
x = np.arange(2, 500e3, 1)
|
||
|
||
def magn(nu, R, L, C):
|
||
w = 2*np.pi*nu
|
||
return R/np.sqrt(R**2+(w*L+1/(w*C))**2)
|
||
f, (Ro, Lo, Co) = curve(nu, abs(H), magn, sigmaHm, guess=[R.n, L.n, C.n])
|
||
|
||
# χ² and compatibility tests
|
||
alpha = chi_squared_fit(nu, abs(H), f, sigma=sigmaHm, s=3)
|
||
beta = [check_measures(*i) for i in zip([R, L, C], [Ro, Lo, Co])]
|
||
|
||
print(mformat('''
|
||
# |H| fit (resistor)
|
||
Rₒ: {} Ω - β={:.2f}
|
||
Lₒ: {} H - β={:.2f}
|
||
Cₒ: {} F - β={:.2f}
|
||
|
||
χ² test:
|
||
α={:.2f}, α>ε: {}
|
||
''', Ro, beta[0]
|
||
, Lo, beta[1]
|
||
, Co, beta[2]
|
||
, alpha, alpha>epsilon))
|
||
|
||
plt.figure(1)
|
||
plt.clf()
|
||
plt.subplot(2, 1, 1)
|
||
plt.title('transfer function (R)')
|
||
plt.xscale('log')
|
||
plt.ylabel('magnitude')
|
||
plt.scatter(nu, abs(H), color='#8c4f4a')
|
||
plt.plot(x, f(x), color='#845336')
|
||
|
||
|
||
# phase, fit and plot
|
||
def phase(nu, R, L, C):
|
||
w = 2*np.pi*nu
|
||
return np.arctan((w*L-1/(w*C))/R)
|
||
f, (Ro, Lo, Co) = curve(nu, np.angle(H), phase, sigmaHp, guess=[R.n, L.n, C.n])
|
||
|
||
# χ² and compatibility tests
|
||
alpha = chi_squared_fit(nu, np.angle(H), f, sigma=sigmaHp, s=3)
|
||
beta = [check_measures(*i) for i in zip([R, L, C], [Ro, Lo, Co])]
|
||
|
||
print(mformat('''
|
||
# ∠H fit (resistor)
|
||
Rₒ: {} Ω - β={:.2f}
|
||
Lₒ: {} H - β={:.2f}
|
||
Cₒ: {} F - β={:.2f}
|
||
|
||
χ² test:
|
||
α={:.2f}, α>ε: {}
|
||
''', Ro, beta[0]
|
||
, Lo, beta[1]
|
||
, Co, beta[2]
|
||
, alpha, alpha>epsilon))
|
||
|
||
plt.subplot(2, 1, 2)
|
||
plt.xscale('log')
|
||
plt.scatter(nu, np.angle(H), color='#43454f')
|
||
plt.plot(x, f(x), color='#65788f')
|
||
plt.xlabel('frequency (Hz)')
|
||
plt.ylabel('phase (rad)')
|
||
plt.show()
|
||
|
||
|
||
## Capacitor
|
||
|
||
# frequency (Hz)
|
||
nu = array( 10, 30, 50, 100, 150, 200, 300, 400, 500, 600, 800,
|
||
1e3, 1.2e3, 1.6e3, 2.5e3, 3e3, 5e3, 10e3, 20e3, 30e3, 50e3, 100e3)
|
||
|
||
# input tension (V)
|
||
Va = array(10.2, 10.2, 10.2, 10.0, 10.0, 9.8, 9.8, 9.8, 9.8, 9.8, 9.6,
|
||
9.8, 9.8, 9.8, 9.8, 10.0, 10.0, 10.2, 10.0, 10.4, 10.4, 10.4)
|
||
|
||
# output tension (V)
|
||
Vb = array(10.2, 9.8, 9.8, 8.8, 8.0, 6.8, 5.2, 4.2, 3.4, 3.0, 2.16, 1.8,
|
||
1.4, 1.2, 800e-3, 600e-3, 360e-3, 50.4e-3, 82.4e-3, 176e-3, 24e-3, 8.8e-3)
|
||
|
||
# phase difference Va - Vb (s)
|
||
Oab = array(1.6e-3, 1.1e-3, 1e-3, 840e-6, 760e-6, 720e-6, 540e-6, 460e-6,
|
||
390e-6, 350e-6, 272e-6, 230e-6, 200e-6, 152e-6, 104e-6, 86e-6,
|
||
53e-6, 29.6e-6, 16.9e-6, 12.3e-6, 8.2e-6, 4.6e-6)
|
||
|
||
# phase error (s)
|
||
sigmaOab = array(0.5e-3, 0.5e-3, 0.5e-3, 10e-6, 10e-6, 10e-6, 10e-6, 10e-6,
|
||
10e-6, 10e-6, 5e-6, 5e-6, 5e-6, 5e-6, 1e-6, 1e-6,
|
||
1e-6, 1e-6, 0.3e-6, 0.2e-6, 0.2e-6, 0.1e-6)
|
||
|
||
# transfer function
|
||
H = Vb/Va * np.exp(2j*np.pi*nu * Oab)
|
||
|
||
# estimate uncertainties
|
||
sigmaHm = (ufloat(Vb[0], 0.1)/ufloat(Va[0], 0.1)).s
|
||
sigmaHp = array(*((2*np.pi*v * ufloat(o, so)).s for v, o, so in zip(nu, Oab, sigmaOab)))
|
||
|
||
def magn(nu, R, L, C):
|
||
w = 2*np.pi*nu
|
||
return 1/(w*C*np.sqrt(R**2 + (w*L - 1/(w*C))**2))
|
||
f, (Ro, Lo, Co) = curve(nu, abs(H), magn, sigmaHm, guess=[R.n, L.n, C.n])
|
||
|
||
# χ² and compatibility tests
|
||
alpha = chi_squared_fit(nu, abs(H), f, sigma=sigmaHm, s=3)
|
||
beta = [check_measures(*i) for i in zip([R, L, C], [Ro, Lo, Co])]
|
||
|
||
print(mformat('''
|
||
# |H| fit (capacitor)
|
||
Rₒ: {} Ω - β={:.2f}
|
||
Lₒ: {} H - β={:.2f}
|
||
Cₒ: {} F - β={:.2f}
|
||
|
||
χ² test:
|
||
α={:.2f}, α>ε: {}
|
||
''', Ro, beta[0]
|
||
, Lo, beta[1]
|
||
, Co, beta[2]
|
||
, alpha, alpha>epsilon))
|
||
|
||
plt.figure(2)
|
||
plt.clf()
|
||
plt.subplot(2, 1, 1)
|
||
plt.title('transfer function (C)')
|
||
plt.xscale('log')
|
||
plt.ylabel('magnitude')
|
||
plt.scatter(nu, abs(H), color='#57553c')
|
||
plt.plot(x, f(x), color='#898471')
|
||
|
||
# phase, fit and plot
|
||
def phase(nu, R, L, C):
|
||
w = 2*np.pi*nu
|
||
return np.arctan((L*w - 1/(C*w))/R)+np.pi/2
|
||
f, (Ro, Lo, Co) = curve(nu, np.angle(H), phase, sigmaHp, guess=[R.n, L.n, C.n])
|
||
|
||
# χ² and compatibility tests
|
||
alpha = chi_squared_fit(nu, np.angle(H), f, sigma=sigmaHp, s=3)
|
||
beta = [check_measures(*i) for i in zip([R, L, C], [Ro, Lo, Co])]
|
||
|
||
print(mformat('''
|
||
# ∠H fit (capacitor)
|
||
Rₒ: {} Ω - β={:.2f}
|
||
Lₒ: {} H - β={:.2f}
|
||
Cₒ: {} F - β={:.2f}
|
||
|
||
χ² test:
|
||
α={:.2f}, α>ε: {}
|
||
''', Ro, beta[0]
|
||
, Lo, beta[1]
|
||
, Co, beta[2]
|
||
, alpha, alpha>epsilon))
|
||
|
||
plt.subplot(2, 1, 2)
|
||
plt.xscale('log')
|
||
plt.scatter(nu, np.angle(H), color='#a17e3e')
|
||
plt.plot(x, f(x), color='#c8b491')
|
||
plt.xlabel('frequency (Hz)')
|
||
plt.ylabel('phase (rad)')
|
||
plt.show()
|
||
|
||
|
||
## Inductor
|
||
|
||
# frequency (Hz)
|
||
nu = array(10, 50, 100, 200, 300, 500, 800, 1e3, 1.5e3, 2e3, 3e3, 5e3,
|
||
8e3, 10e3, 15e3, 20e3, 30e3, 50e3, 80e3, 100e3, 120e3, 150e3, 200e3, 300e3, 500e3)
|
||
|
||
# input tension (V)
|
||
Va = array( 6.4, 10.4, 10.0, 10.2, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0,
|
||
10.2, 10.2, 10.4, 10.2, 10.4, 10.4, 10.6, 10.4, 10.4, 10.4, 10.4, 10.4, 10.8)/2
|
||
|
||
# output tension (V)
|
||
Vb = array(6e-3, 54e-3, 98e-3, 148e-3, 154e-3, 232e-3, 292e-3, 336e-3, 448e-3, 560e-3, 816e-3,
|
||
1.26, 2.06, 2.52, 3.56, 4.64, 6.16, 7.92, 9.4, 10.0, 10.2, 10.4, 10.4, 10.4, 10.8)/2
|
||
|
||
# output tension error (V)
|
||
sigmaVb = array(1e-3, 2e-3, 2e-3, 5e-3, 5e-3, 10e-3, 10e-3, 10e-3, 15e-3, 15e-3, 20e-3,
|
||
0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.2, 0.2, 0.2, 0.2,
|
||
0.2, 0.2, 0.2)
|
||
|
||
# phase difference Va - Vb (s)
|
||
Oab = -array(25.2e-3, 4.3e-3, 1.98e-3, 820e-6, 512e-6, 290e-6, 208e-6,
|
||
172e-6, 120e-6, 98e-6, 71e-6, 42e-6, 26e-6, 19.6e-6,
|
||
11.8e-6, 8.2e-6, 4.7e-6, 1.9e-6, 840e-9, 560e-9, 320e-9,
|
||
200e-9, 120e-9, 50e-9, 20e-9)
|
||
|
||
# phase error (s)
|
||
sigmaOab = array(1e-3, 0.5e-3, 0.5e-3, 20e-6, 15e-6, 10e-6, 10e-6,
|
||
10e-6, 5e-6, 5e-6, 5e-6, 2e-6, 2e-6, 1e-6,
|
||
1e-6, 0.5e-6, 0.2e-6, 0.2e-6, 40e-9, 20e-9, 20e-9,
|
||
10e-9, 10e-9, 5e-9, 2e-9)
|
||
|
||
# transfer function
|
||
H = Vb/Va * np.exp(2j*np.pi*nu * Oab)
|
||
|
||
# estimate uncertainties
|
||
sigmaHm = array(*((ufloat(vb, svb)/ufloat(va, 0.1)).s for va, vb, svb in zip(Va, Vb, sigmaVb)))
|
||
sigmaHp = array(*((2*np.pi*v * ufloat(o, so)).s for v, o, so in zip(nu, Oab, sigmaOab)))
|
||
exit()
|
||
# magnitude, plot and fit
|
||
def magn(nu, R, Rl, L, C):
|
||
w = 2*np.pi*nu
|
||
return C*w*np.sqrt(w**2*(C*L**2*w**2 + C*R*Rl + C*Rl**2 - L)**2 + (C*L*R*w**2 + Rl)**2)/(C**2*w**2*(R + Rl)**2 + (C*L*w**2 - 1)**2)
|
||
f, (Ro, Rlo, Lo, Co) = curve(nu, abs(H), magn, guess=[R.n, Rl.n, L.n, C.n])
|
||
|
||
# χ² and compatibility tests
|
||
alpha = chi_squared_fit(nu, abs(H), f, sigma=sigmaHm, s=4)
|
||
beta = [check_measures(*i) for i in zip([R, Rl, L, C], [Ro, Rlo, Lo, Co])]
|
||
|
||
print(mformat('''
|
||
# |H| fit (inductor)
|
||
Rₒ: {} Ω - β={:.2f}
|
||
Lₒ: {} H - β={:.2f}
|
||
Cₒ: {} F - β={:.2f}
|
||
|
||
χ² test:
|
||
α={:.2f}, α>ε: {}
|
||
''', Ro, beta[0]
|
||
, Lo, beta[1]
|
||
, Co, beta[2]
|
||
, alpha, alpha>epsilon))
|
||
|
||
plt.figure(3)
|
||
plt.subplot(2, 1, 1)
|
||
plt.title('transfer function (L)')
|
||
plt.xscale('log')
|
||
plt.ylabel('magnitude')
|
||
plt.scatter(nu, abs(H), color='#a18b62')
|
||
plt.plot(x, f(x), color='#bc9d66')
|
||
|
||
# phase, fit and plot
|
||
def phase(nu, R, Rl, L, C):
|
||
w = 2*np.pi*nu
|
||
return np.arctan2(C*w**2*(C*Rl*(R + Rl) + L*(C*L*w**2 - 1))/(C**2*w**2*(R + Rl)**2 + (C*L*w**2 - 1)**2),
|
||
C*w*(C*L*w**2*(R + Rl) - Rl*(C*L*w**2 - 1))/(C**2*w**2*(R + Rl)**2 + (C*L*w**2 - 1)**2))
|
||
f, (Ro, Rlo, Lo, Co) = curve(nu, np.angle(H), phase, sigmaHp, guess=[R.n, Rl.n, L.n, C.n])
|
||
|
||
# χ² and compatibility tests
|
||
alpha = chi_squared_fit(nu, np.angle(H), f, sigma=sigmaHp, s=4)
|
||
beta = [check_measures(*i) for i in zip([R, Rl, L, C], [Ro, Rlo, Lo, Co])]
|
||
|
||
print(mformat('''
|
||
# ∠H fit (inductor)
|
||
Rₒ: {} Ω - β={:.2f}
|
||
Lₒ: {} H - β={:.2f}
|
||
Cₒ: {} F - β={:.2f}
|
||
|
||
χ² test:
|
||
α={:.2f}, α>ε: {}
|
||
''', Ro, beta[0]
|
||
, Lo, beta[1]
|
||
, Co, beta[2]
|
||
, alpha, alpha>epsilon))
|
||
|
||
plt.subplot(2, 1, 2)
|
||
plt.xscale('log')
|
||
plt.scatter(nu, np.angle(H), color = '#43454f')
|
||
plt.plot(x, f(x), color='#65788f')
|
||
plt.xlabel('frequency (Hz)')
|
||
plt.ylabel('phase (rad)')
|
||
plt.show()
|