265 lines
6.3 KiB
Python
265 lines
6.3 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 *
|
|||
|
|
|
|||
|
|
##
|
|||
|
|
## Impedence of an inductor (I)
|
|||
|
|
## (all SI units)
|
|||
|
|
|
|||
|
|
|
|||
|
|
## measured quantities
|
|||
|
|
R = ufloat(996, 4) # resistor
|
|||
|
|
Rl = ufloat(19.9, 0.1) # internal resistance (inductor)
|
|||
|
|
|
|||
|
|
# frequency
|
|||
|
|
nu = array(60, 100, 200, 500, 800, 1.2e3, 1.6e3, 2e3, 3e3, 5e3,
|
|||
|
|
10e3, 15e3, 20e3, 30e3, 40e3, 50e3, 80e3, 100e3, 200e3, 300e3)
|
|||
|
|
|
|||
|
|
# V input
|
|||
|
|
Va = array(9.60, 9.60, 9.60, 9.60, 9.60, 9.60, 9.60, 9.80, 9.80, 9.80,
|
|||
|
|
9.80, 9.60, 9.60, 9.80, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0)/2
|
|||
|
|
|
|||
|
|
# V output
|
|||
|
|
Vb = array(9.40, 9.40, 9.20, 9.20, 9.20, 9.20, 9.20, 9.40, 9.20, 9.20,
|
|||
|
|
8.30, 7.45, 6.81, 5.77, 5.20, 4.63, 3.45, 2.96, 1.73, 1.21)/2
|
|||
|
|
|
|||
|
|
# time offset Va - Vb
|
|||
|
|
Oab = -array(100e-6, 80e-6, 20e-6, 20e-6, 12e-6, 8e-6, 8e-6, 8e-6, 6e-6, 3e-6,
|
|||
|
|
4.5e-6, 4.1e-6, 3.9e-6, 3.5e-6, 3.1e-6, 2.8e-6, 2.1e-6, 1.8e-6, 1.1e-6, 780e-9)
|
|||
|
|
|
|||
|
|
# time offset I - V
|
|||
|
|
Oiv = array( 3e-4, 2.1e-4, 1e-4, 65e-6, 88e-6, 80e-6, 78e-6, 72e-6, 56e-6, 40e-6,
|
|||
|
|
22e-6, 16.1e-6, 12.2e-6, 8.2e-6, 6.1e-6, 4.9e-6, 3.1e-6, 2.5e-6, 1.26e-6, 840e-9)
|
|||
|
|
|
|||
|
|
|
|||
|
|
## derived quantities
|
|||
|
|
om = 2*np.pi*nu # angular frequency
|
|||
|
|
I = Vb/R.n # output current
|
|||
|
|
Vab = Va - Vb # tension drop
|
|||
|
|
Fab = om * Oab # phase difference Vᵢ - Vₒ
|
|||
|
|
Fiv = om * Oiv # phase difference I - Vₒ
|
|||
|
|
|
|||
|
|
Z = Vab/I * np.exp(1j*Fiv) # impedance
|
|||
|
|
H1 = Vab/Vb * np.exp(1j*Fab) # transfer function ΔV→Vb
|
|||
|
|
H2 = Vb/Va * np.exp(1j*Fab) # transfer function Va→Vb
|
|||
|
|
|
|||
|
|
|
|||
|
|
## estimate uncertainties
|
|||
|
|
Evb, Eva = ufloat(Vb[0], 0.15), ufloat(Va[0], 0.15)
|
|||
|
|
|
|||
|
|
# uncertainties Oiv
|
|||
|
|
sigmaOiv = array(5e-4, 5e-4, 5e-4, 1e-06, 1e-06, 1e-06, 1e-06, 1e-06, 1e-06, 1e-06,
|
|||
|
|
1e-06, 1e-07, 1e-07, 1e-07, 1e-07, 1e-07, 1e-07, 1e-07, 5e-08, 5e-08)
|
|||
|
|
|
|||
|
|
# uncertainties Oab
|
|||
|
|
sigmaOab = array(5e-06, 5e-06, 5e-06, 5e-06, 2e-06, 2e-06, 2e-06, 2e-06, 2e-06, 2e-06,
|
|||
|
|
2e-07, 2e-07, 2e-07, 2e-07, 2e-07, 2e-07, 2e-07, 2e-07, 2e-07, 2e-08)
|
|||
|
|
|
|||
|
|
# phase uncertainties
|
|||
|
|
sigmaFiv = sigmaOiv * 2*np.pi*nu
|
|||
|
|
sigmaFab = sigmaOab * 2*np.pi*nu
|
|||
|
|
|
|||
|
|
# magnitude ucnertainties
|
|||
|
|
sigmaZ = (R*(Eva-Evb)/Evb).s
|
|||
|
|
sigmaH1 = ((Eva-Evb)/Evb).s
|
|||
|
|
sigmaH2 = ((Evb-Eva)/Evb).s
|
|||
|
|
|
|||
|
|
|
|||
|
|
## plot and fit Z(ν)
|
|||
|
|
plt.figure(4)
|
|||
|
|
plt.clf()
|
|||
|
|
|
|||
|
|
# magnitude
|
|||
|
|
plt.subplot(2, 1, 1)
|
|||
|
|
plt.title('impedance (RL circuit)')
|
|||
|
|
plt.xscale('log')
|
|||
|
|
plt.ylabel('magnitude (kΩ)')
|
|||
|
|
plt.scatter(nu, abs(Z)/1e3, color="#36913d")
|
|||
|
|
|
|||
|
|
x = np.arange(nu.min()-10, nu.max()+2e3, 10)
|
|||
|
|
|
|||
|
|
# fit |Z|(ω) = √(Rl² + (iωL)²)
|
|||
|
|
f = lambda x, Rl, L: np.sqrt(Rl**2 + (2*np.pi*L * x)**2)
|
|||
|
|
f, (Rlo, L1) = curve(nu, abs(Z), f, sigmaZ, guess=[Rl.n, 0.04], method='trf')
|
|||
|
|
plt.plot(x, f(x)/1e3, color="#36913d")
|
|||
|
|
|
|||
|
|
alpha = check_measures(Rl, Rlo)
|
|||
|
|
beta = chi_squared_fit(nu, abs(Z), f, sigmaZ)
|
|||
|
|
|
|||
|
|
print(mformat('''
|
|||
|
|
# fit magnitude Z
|
|||
|
|
L₁: {} H
|
|||
|
|
Rlₒ₁: {} Ω
|
|||
|
|
|
|||
|
|
compatibility test Rl/Rlₒ₁:
|
|||
|
|
α={:.2f}, α>ε: {}
|
|||
|
|
|
|||
|
|
χ² test:
|
|||
|
|
β={:.2f}, β>ε: {}
|
|||
|
|
''', L1, Rlo,
|
|||
|
|
alpha, alpha>epsilon,
|
|||
|
|
beta, beta>epsilon))
|
|||
|
|
|
|||
|
|
# phase
|
|||
|
|
plt.subplot(2, 1, 2)
|
|||
|
|
plt.xscale('log')
|
|||
|
|
plt.xlabel('frequency (Hz)')
|
|||
|
|
plt.ylabel('phase (rad)')
|
|||
|
|
plt.scatter(nu, Fiv, color="#33859d")
|
|||
|
|
|
|||
|
|
# fit ∠Z(ω) = arctan(ωL/Rl)
|
|||
|
|
f = lambda x, Rl, L: np.arctan(2*np.pi*x*L/Rl)
|
|||
|
|
f, (Rlo2, L2) = curve(nu, Fiv, f, sigmaFiv, guess=[Rl.n, L1.n])
|
|||
|
|
plt.plot(x, f(x), color="#245361")
|
|||
|
|
plt.show()
|
|||
|
|
|
|||
|
|
alpha = check_measures(Rl, Rlo2)
|
|||
|
|
beta = chi_squared_fit(nu, Fiv, f, sigmaFiv)
|
|||
|
|
|
|||
|
|
print(mformat('''
|
|||
|
|
# fit phase Z
|
|||
|
|
L₂: {} H
|
|||
|
|
Rlₒ₂: {} Ω
|
|||
|
|
|
|||
|
|
compatibility test Rl/Rlₒ₂:
|
|||
|
|
α={:.2f}, α>ε: {}
|
|||
|
|
|
|||
|
|
χ² test:
|
|||
|
|
β={:.2f}, β>ε: {}
|
|||
|
|
''', L2, Rlo2,
|
|||
|
|
alpha, alpha>epsilon,
|
|||
|
|
beta, beta>epsilon))
|
|||
|
|
|
|||
|
|
|
|||
|
|
## plot, fit H₁(ν)
|
|||
|
|
plt.figure(5)
|
|||
|
|
plt.clf()
|
|||
|
|
|
|||
|
|
# magnitude
|
|||
|
|
plt.subplot(2,1,1)
|
|||
|
|
plt.title('transmission function 1')
|
|||
|
|
plt.xscale('log')
|
|||
|
|
plt.ylabel('amplitude (Vout-Vin / Vout)')
|
|||
|
|
plt.scatter(nu, abs(H1), color="#2e3340")
|
|||
|
|
|
|||
|
|
# fit |H₁|(ω) = √(Rl² + (iωL)²)/R
|
|||
|
|
f = lambda x, R, L: np.sqrt(Rl.n**2 + (2*np.pi*L*x)**2)/R
|
|||
|
|
f, (Ro1, L1) = curve(nu, abs(H1), f, sigmaH1, guess=[R.n, L1.n]) # NB. does not converge adding Rl
|
|||
|
|
plt.plot(x, f(x), color="#61778d")
|
|||
|
|
|
|||
|
|
alpha = check_measures(R, Ro1)
|
|||
|
|
beta = chi_squared_fit(nu, abs(H1), f, sigmaH1)
|
|||
|
|
|
|||
|
|
print(mformat('''
|
|||
|
|
# fit magnitude H₁
|
|||
|
|
L₁: {} H
|
|||
|
|
Rₒ₁: {} Ω
|
|||
|
|
|
|||
|
|
compatibility test R/Rₒ₁:
|
|||
|
|
α={:.2f}, α>ε: {}
|
|||
|
|
|
|||
|
|
χ² test:
|
|||
|
|
γ={:.2f}, γ>ε: {}
|
|||
|
|
''', L1, Ro1,
|
|||
|
|
alpha, alpha>epsilon,
|
|||
|
|
beta, beta>epsilon))
|
|||
|
|
|
|||
|
|
# phase
|
|||
|
|
plt.subplot(2,1,2)
|
|||
|
|
plt.xscale('log')
|
|||
|
|
plt.xlabel('frequency (Hz)')
|
|||
|
|
plt.ylabel('phase (rad)')
|
|||
|
|
plt.scatter(nu, Fab, color="#a54242")
|
|||
|
|
|
|||
|
|
# fit ∠H₁(ω) = arctan(-ωL/(R+Rl))
|
|||
|
|
f = lambda x, Re, L: np.arctan(-2*np.pi*L/Re * x)
|
|||
|
|
f, (Re, L2) = curve(nu, Fab, f, sigmaFab, guess=[Rl.n+R.n, L1.n])
|
|||
|
|
plt.plot(x, f(x), color='#cc6666')
|
|||
|
|
plt.show()
|
|||
|
|
|
|||
|
|
alpha = check_measures(Rl+R, Re)
|
|||
|
|
beta = chi_squared_fit(nu, Fab, f, 2*sigmaFab)
|
|||
|
|
|
|||
|
|
print(mformat('''
|
|||
|
|
# fit phase H₁
|
|||
|
|
L₂: {} H
|
|||
|
|
Re: {} Ω
|
|||
|
|
|
|||
|
|
compatibility test (Rl+R)/Re:
|
|||
|
|
α={:.2f}, α>ε: {}
|
|||
|
|
|
|||
|
|
χ² test:
|
|||
|
|
β={:.2f}, γ>ε: {}
|
|||
|
|
''', L2, Re,
|
|||
|
|
alpha, alpha>epsilon,
|
|||
|
|
beta, beta>epsilon))
|
|||
|
|
|
|||
|
|
|
|||
|
|
## plot, fit H₂(ν)
|
|||
|
|
plt.figure(6)
|
|||
|
|
plt.clf()
|
|||
|
|
|
|||
|
|
# magnitude
|
|||
|
|
plt.subplot(2,1,1)
|
|||
|
|
plt.title('transmission function 2')
|
|||
|
|
plt.xscale('log')
|
|||
|
|
plt.ylabel('amplitude (Vout / Vin)')
|
|||
|
|
plt.scatter(nu, abs(H2), color='#845336')
|
|||
|
|
|
|||
|
|
# fit |H₂|(ω) = R/√((Rl+R)² + (ωL)²) ≈ 1/√(1 + (iωL)²)
|
|||
|
|
f = lambda x, R, L: 1/np.sqrt(1 + (2*np.pi*L*x/R)**2)
|
|||
|
|
f, (Ro1, L1) = curve(nu, abs(H2), f, sigmaH2, guess=[R.n, L1.n])
|
|||
|
|
plt.plot(x, f(x), color="#8c4f4a")
|
|||
|
|
|
|||
|
|
alpha = check_measures(R, Ro1)
|
|||
|
|
beta = chi_squared_fit(nu, abs(H2), f, sigmaH2)
|
|||
|
|
|
|||
|
|
print(mformat('''
|
|||
|
|
# fit magnitude H₂
|
|||
|
|
L₁: {} H
|
|||
|
|
Rₒ₁: {} Ω
|
|||
|
|
|
|||
|
|
compatibility test R/Rₒ₁:
|
|||
|
|
α={:.2f}, α>ε: {}
|
|||
|
|
|
|||
|
|
χ² test:
|
|||
|
|
γ={:.2f}, γ>ε: {}
|
|||
|
|
''', L1, Ro1,
|
|||
|
|
alpha, alpha>epsilon,
|
|||
|
|
beta, beta>epsilon))
|
|||
|
|
|
|||
|
|
# phase
|
|||
|
|
plt.subplot(2,1,2)
|
|||
|
|
plt.xscale('log')
|
|||
|
|
plt.xlabel('frequency (Hz)')
|
|||
|
|
plt.ylabel('phase (rad)')
|
|||
|
|
plt.scatter(nu, Fab, color='#5c6652')
|
|||
|
|
|
|||
|
|
# fit ∠H₁(ω) = arctan(-ωL/(R+Rl))
|
|||
|
|
f = lambda x, Re, L: np.arctan(-2*np.pi*L/Re * x)
|
|||
|
|
f, (Re, L2) = curve(nu, Fab, f, sigmaFab, guess=[Rl.n+R.n, L1.n])
|
|||
|
|
plt.plot(x, f(x), color='#718062')
|
|||
|
|
plt.show()
|
|||
|
|
|
|||
|
|
alpha = check_measures(Rl+R, Re)
|
|||
|
|
beta = chi_squared_fit(nu, Fab, f, 2*sigmaFab)
|
|||
|
|
|
|||
|
|
print(mformat('''
|
|||
|
|
# fit phase H₂
|
|||
|
|
L₂: {} H
|
|||
|
|
Re: {} Ω
|
|||
|
|
|
|||
|
|
compatibility test (Rl+R)/Re:
|
|||
|
|
α={:.2f}, α>ε: {}
|
|||
|
|
|
|||
|
|
χ² test:
|
|||
|
|
β={:.2f}, γ>ε: {}
|
|||
|
|
''', L2, Re,
|
|||
|
|
alpha, alpha>epsilon,
|
|||
|
|
beta, beta>epsilon))
|
|||
|
|
plt.show()
|