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()
|