lab-II/circuits/C3-RL2.py

265 lines
6.1 KiB
Python
Raw Permalink Normal View History

2018-03-18 18:02:21 +01:00
# 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 (II)
## (all SI units)
## measured quantities
R = ufloat(996, 4) # resistor
Rl = ufloat(64.8, 0.1) # internal resitance (inductor)
# frequency
nu = array(60, 100, 200, 500, 800, 1.2e3, 1.6e3, 2e3, 3e3, 5e3,
10e3, 15e3, 20e3, 40e3, 50e3)
# V input
Va = array(9.60, 9.60, 9.60, 9.60, 9.60, 9.80, 9.80,
10.00, 10.20, 10.20, 10.20, 10.20, 10.20, 10.00, 10.00)/2
# V output
Vb = array(9.00, 9.00, 9.00, 8.60, 8.00, 7.40, 6.60, 6.00,
4.40, 3.00, 1.60, 1.10, 0.80, 0.40, 0.32)/2
# time offset Va - Vb
Oab = array(200e-6, 40e-6, 80e-6, 90e-6, 90e-6, 84e-6, 70e-6,
69e-6, 58e-6, 38e-6, 22e-6, 15e-6, 12e-6, 6.2e-6, 4.8e-6)
# time offset I - V
Oiv = array(880e-6, 500e-6, 620e-6, 370e-6, 260e-6, 184e-6, 140e-6,
118e-6, 83e-6, 48e-6, 25e-6, 17e-6, 12.4e-6, 6.6e-6, 5.1e-6)
## 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.2), ufloat(Va[0], 0.2)
# uncertainties Oiv
sigmaOiv = array(2e-4, 5e-6, 5e-6, 5e-6, 5e-6, 5e-6, 5e-6, 5e-6,
1e-6, 1e-6, 1e-6, 1e-6, 1e-6, 1e-6, 1e-6)
# uncertainties Oab
sigmaOab = array(1e-5, 1e-6, 1e-6, 1e-6, 1e-6, 1e-6, 1e-6, 1e-6,
1e-6, 1e-6, 1e-6, 1e-6, 2e-7, 2e-7, 2e-7)
# 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, (Rlo1, 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, Rlo1)
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, Rlo1,
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()