From f0fa70042718d4085764486f9f135a9f30ae9afe Mon Sep 17 00:00:00 2001 From: rnhmjoj Date: Sun, 18 Mar 2018 18:02:21 +0100 Subject: [PATCH] initial commit --- README.md | 27 ++ circuits/C1.py | 179 ++++++++++++ circuits/C2.py | 274 ++++++++++++++++++ circuits/C3-?.py | 188 +++++++++++++ circuits/C3-RC1.py | 205 ++++++++++++++ circuits/C3-RC2.py | 204 ++++++++++++++ circuits/C3-RL1.py | 264 ++++++++++++++++++ circuits/C3-RL2.py | 264 ++++++++++++++++++ circuits/C4-1C.py | 263 ++++++++++++++++++ circuits/C4-1D.py | 258 +++++++++++++++++ circuits/C4-2.py | 314 +++++++++++++++++++++ circuits/__pycache__/lab.cpython-36.pyc | Bin 0 -> 10552 bytes circuits/lab.py | 1 + circuits/lab.pyc | Bin 0 -> 12721 bytes circuits/transfer.py | 59 ++++ lab.py | 351 ++++++++++++++++++++++++ optics/__pycache__/lab.cpython-36.pyc | Bin 0 -> 10550 bytes optics/interf.py | 148 ++++++++++ optics/lab.py | 1 + optics/lab.pyc | Bin 0 -> 12657 bytes optics/light.py | 173 ++++++++++++ optics/micro.py | 339 +++++++++++++++++++++++ optics/spectra.py | 192 +++++++++++++ optics/spectra/ar | 24 ++ optics/spectra/he | 10 + optics/spectra/hg | 10 + optics/spectra/kr | 48 ++++ optics/spectra/n | 5 + optics/spectra/na | 9 + optics/spectra/ne | 63 +++++ optics/spectra/xe | 66 +++++ shell.nix | 22 ++ 32 files changed, 3961 insertions(+) create mode 100644 README.md create mode 100644 circuits/C1.py create mode 100644 circuits/C2.py create mode 100644 circuits/C3-?.py create mode 100644 circuits/C3-RC1.py create mode 100644 circuits/C3-RC2.py create mode 100644 circuits/C3-RL1.py create mode 100644 circuits/C3-RL2.py create mode 100644 circuits/C4-1C.py create mode 100644 circuits/C4-1D.py create mode 100644 circuits/C4-2.py create mode 100644 circuits/__pycache__/lab.cpython-36.pyc create mode 120000 circuits/lab.py create mode 100644 circuits/lab.pyc create mode 100644 circuits/transfer.py create mode 100644 lab.py create mode 100644 optics/__pycache__/lab.cpython-36.pyc create mode 100644 optics/interf.py create mode 120000 optics/lab.py create mode 100644 optics/lab.pyc create mode 100644 optics/light.py create mode 100644 optics/micro.py create mode 100644 optics/spectra.py create mode 100644 optics/spectra/ar create mode 100644 optics/spectra/he create mode 100644 optics/spectra/hg create mode 100644 optics/spectra/kr create mode 100644 optics/spectra/n create mode 100644 optics/spectra/na create mode 100644 optics/spectra/ne create mode 100644 optics/spectra/xe create mode 100644 shell.nix diff --git a/README.md b/README.md new file mode 100644 index 0000000..3d6861a --- /dev/null +++ b/README.md @@ -0,0 +1,27 @@ +# Esperimenti di Laboratorio II +## Programmi usati per svolgere i calcoli degli esperimenti + +### Info +Ogni modulo rappresenta un esperimento e produce come output i risultati +richiesti (grafici, valori, tabelle...). +Il modulo [lab](https://github.com/rnhmjoj/lab-II/blob/master/lab.py) esporta +funzioni statistiche usate usate nei vari programmi. Usare `help()` per ottenere +la documentazione specifica ad ogni funzione. + +### Dipendenze +I programmi sono in python 2 (ma facilmente portabile a 3) e richiedono i seguenti +pacchetti di [PyPi](https://pypi.python.org/pypi): + + 1. [numpy](https://pypi.python.org/pypi/numpy) + 2. [matplotlib](https://pypi.python.org/pypi/matplotlib) + 3. [sympy](https://pypi.python.org/pypi/sympy) + 4. [scipy](https://pypi.python.org/pypi/scipy) + 4. [uncertainties](https://pypi.python.org/pypi/uncertainties) + +È inclusa una nix expression ([shell.nix](https://github.com/rnhmjoj/lab-II/blob/master/shell.nix)) +che definisce l'ambiente minimo per eseguire il codice. + +### License +Dual licensed under the MIT and GPL licenses: +http://www.opensource.org/licenses/mit-license.php +http://www.gnu.org/licenses/gpl.html diff --git a/circuits/C1.py b/circuits/C1.py new file mode 100644 index 0000000..88cf453 --- /dev/null +++ b/circuits/C1.py @@ -0,0 +1,179 @@ +# coding: utf-8 +from __future__ import print_function, division, unicode_literals + +from lab import * +import numpy as np +import uncertainties.umath as um +import matplotlib.pyplot as plt + + +## +## Internal resistance +## + +Req = array(2.021, 3.048) / \ + array(2.25e-6, 1.19e-6) # potential/current +Rn = array(0.98e6, 3.322e6) # nominal value +Rv = sample(Rn*Req/(Rn - Req)).val() # internal resistance + +Req = array(1.101, 1.104, 1.108) / \ + array(70.08e-3, 34.545e-3, 23.555e-3) # potential/current +Rn = array(9.88, 21.66, 32.64) # nominal value +Ra = sample(Req - Rn).val() # internal resistance + +print(mformat(""" +Internal resistance: +{} Ω (voltmeter) +{} Ω (amperometer) +""", Ra, Rv)) + + +## +## Ohm Law +## + +# 10Ω resistor +I = array(44.10, 52.65, 61.60, 66.00, 70.60, 75.07, 79.12, 83.66, + 88.39, 92.41, 96.44, 101.28, 105.69, 109.83, 114.23, 118.67, + 123.35, 127.40, 132.13, 140.52)/1e3 +V = array(0.432, 0.515, 0.604, 0.646, 0.691, 0.734, 0.774, 0.820, 0.864, 0.904, + 0.944, 0.991, 1.034, 1.074, 1.117, 1.161, 1.207, 1.247, 1.293, 1.375) + +# linear fit y=a+bx where y=V and x=I +x = np.linspace(40, 145, 100)/1e3 +a, b = linear(I, V, 1e-3) +f = lambda x: a.n + b.n*x +R = b + +alpha = chi_squared_fit(I, V, f, sigma=1e-3) + +# plot y=a+bx +plt.figure(1) +plt.xlabel('I (A)') +plt.ylabel('dV (V)') +plt.scatter(I, V, color='#ba6b88') +plt.plot(x, f(x), '#ba6b88') +plt.show() + +print(mformat(''' +R={} Ω + +χ² test: + α={:.3f}, α>ε: {} +''', R, alpha, alpha>epsilon)) + + +# 1kΩ resistor +V = array(0.738, 0.828, 0.922, 1.011, 1.103, 1.198, 1.288, 1.377, 1.472, 1.568, + 1.658, 1.751, 1.846, 1.936, 2.028, 2.114, 2.213, 2.300, 2.391, 2.487) +I = array(0.7445, 0.8350, 0.9304, 1.0194, 1.1127, 1.2089, 1.2993, 1.3897, 1.4850, + 1.5824, 1.6732, 1.7670, 1.8625, 1.9537, 2.0474, 2.1339, 2.2334, 2.3213, + 2.4134, 2.5099)/1e3 + +# linear fit y=a+bx where y=V and x=I +x = np.linspace(0.7, 2.6, 100)/1e3 +a, b = linear(I, V, 1e-3) +f = lambda x: a.n + b.n*x +R = b + +alpha = chi_squared_fit(I, V, f, sigma=1e-3) + +# plot y=a+bx +plt.figure(2) +plt.xlabel('I (mA)') +plt.ylabel('dV (V)') +plt.scatter(I*1e3, V, color='#3e49ba') +plt.plot(x*1e3, f(x), '#3e49ba') +plt.show() + +print(mformat(''' +R={} Ω + +χ² test: + α={:.3f}, α>ε: {} +''', R, alpha, alpha>epsilon)) + + +# 1MΩ resistor + +V = array(1.012, 1.118, 1.214, 1.317, 1.421, 1.524, 1.621, 1.722, 1.823, 1.924, + 2.030, 2.130, 2.230, 2.330, 2.434, 2.534, 2.640, 2.739, 2.845, 2.940) +I = array(1.12, 1.24, 1.34, 1.46, 1.58, 1.69, 1.80, 1.91, 2.03, 2.14, 2.26, 2.37, + 2.47, 2.58, 2.71, 2.82, 2.94, 3.04, 3.16, 3.26)/1e6 + +# linear fit y=a+bx where y=V and x=I +x = np.linspace(1.10, 3.30, 100)/1e6 +a, b = linear(I, V, 5e-3) +f = lambda x: a.n + b.n*x +R = b + +alpha = chi_squared_fit(I, V, f, sigma=5e-3) + +# plot y=a+bx +plt.figure(3) +plt.xlabel('I (muA)') +plt.ylabel('dV (V)') +plt.scatter(I*1e6, V, color='#d5b647') +plt.plot(x*1e6, f(x), '#d5b647') +plt.show() + +print(mformat(''' +R={} Ω + +χ² test: + α={:.3f}, α>ε: {} +''', R, alpha, alpha>epsilon)) + + +## +## Diode +## + +V = array(3.172, 3.293, 3.393, 3.471, 3.529, 3.579, 3.619, 3.651) +I = array(1.22, 86.77, 323.60, 749.20, 1260.20, 1880.90, 2597.5, 3351.45)/1e6 + +sigmaI = 3e-5 + +# fit +x = np.linspace(3.0, 3.7, 100) +shockley = lambda V, a, b, c: a*V + b*(np.exp(c*V) - 1) +f, (a, b, c) = curve(V, I, shockley, sigmaI, guess=[1, 1e-10, 10]) + +y = np.linspace(0, 5000, 100)/1e6 +threshold = lambda I, Vth: Vth +h, (Vth,) = curve(I[4:], V[4:], threshold, sigmaI) + +# parameters +e = 1.602e-19 # elementary charge +k = 1.380e-23 # boltzmann constant +T = 298 # temperature +g = e/(k*T*c) + +# χ² test +alpha = chi_squared_fit(V, I, f, sigmaI, s=3) + +print(mformat(''' +a: {} Ω⁻¹ +I₀: {} A +g: {} +Vth: {} V + +χ² test: + α={:.3f}, α>ε: {} +''', a, b, g, Vth, alpha, alpha>epsilon)) + +# plot I - V +plt.figure(4) +plt.title('Shockley law') +plt.xlabel('ΔV (V)') +plt.ylabel('I (mA)') + +# shockley curve +plt.scatter(V, I*1e3, color='#3363aa') +plt.plot(x, f(x)*1e3, color='#2a518c') + +# threshold line +plt.plot(np.repeat(h(y), len(y)), y*1e3, color='#91375b') +plt.text(Vth.n+0.01, 0.01, 'Vth') +plt.show() + diff --git a/circuits/C2.py b/circuits/C2.py new file mode 100644 index 0000000..d39e81a --- /dev/null +++ b/circuits/C2.py @@ -0,0 +1,274 @@ +# 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 * + +## +## Impedance of a resistor +## (all SI units) + +R = ufloat(996, 4) # resistor (for oscilloscope) + +# frequency +nu = array(60, 100, 150, 200, 300, 500, 750, 1e3, 1.5e3, 2e3, 3e3, + 5e3, 8e3, 12e3, 15e3, 20e3, 25e3, 40e3, 50e3, 60e3, 75e3, + 100e3, 150e3, 200e3, 250e3, 300e3, 350e3) + +# tension drop (via voltmeter) +Vm = array(1.564, 1.564, 1.562, 1.558, 1.563, 1.552, 1.536, 1.518, + 1.470, 1.417, 1.303, 1.096, 0.878, 0.715, 0.644, 0.575, + 0.533, 0.470, 0.445, 0.423, 0.395, 0.349, 0.257, 0.169, + 0.085, 0.038, 0.020)*np.sqrt(2) + +# current (via amperometer) +Im = array(1.663, 1.664, 1.663, 1.664, 1.664, 1.664, 1.664, 1.664, + 1.664, 1.664, 1.664, 1.666, 1.671, 1.679, 1.687, 1.702, + 1.721, 1.798, 1.865, 1.943, 2.076, 2.327, 2.817, 3.143, + 3.255, 3.294, 3.301)*1e-3*np.sqrt(2) + +# tension input - ground (via oscilloscope) +Va = array(9.6, 9.6, 9.6, 9.6, 9.6, 9.6, 9.6, 9.6, 9.6, 9.6, + 9.6, 9.6, 9.6, 9.6, 9.6, 9.6, 9.6, 9.6, 9.6, 9.6, + 9.6, 9.6, 9.6, 9.6, 9.6, 9.6, 9.6)/2 + +# tension output - ground (via oscilloscope) +Vb = array(4.8, 4.8, 4.8, 4.8, 4.8, 4.8, 4.8, 4.8, 4.8, 4.8, + 4.8, 4.8, 4.8, 4.8, 4.8, 4.8, 4.8, 4.8, 4.8, 4.8, + 4.8, 4.8, 4.8, 4.8, 4.8, 4.8, 4.8)/2 + +# add uncertainties +Vb, Va = ufloat(Vb[0], 0.01), ufloat(Va[0], 0.01) + +Io = Vb/R # current (via oscilloscope) +Vo = Va - Vb # tension drop (via oscilloscope) +Z = Vo/Io # impedance + +plt.figure(1) +plt.clf() +plt.xlabel("frequency (Hz)") +plt.ylabel("voltage voltmeter / oscilloscope") +plt.semilogx(nu, Vm/Vo.n, 'o-', color='#bb3e5f', markersize=4.5) +plt.show() + +plt.figure(2) +plt.clf() +plt.xlabel("frequency (Hz)") +plt.ylabel("current amperometer / oscilloscope") +plt.semilogx(nu, Io.n/Im, 'o-', color="#56ad5e", markersize=4.5) +plt.show() + +plt.figure(3) +plt.clf() +plt.title('impedance (purely resistive circuit)') +plt.xlabel('frequency (Hz)') +plt.ylabel('magnitude (kΩ)') +plt.semilogx(nu, nu/nu*Z.n/1e3, 'o-', color="#d0aa23", markersize=4.5) +plt.show() + +alpha = check_measures(Z, R) + +print(mformat(''' +Z: {:.3}Ω +R: {:.3}Ω + +compatibility test: + α={:.2f}, α>ε: {} +''', Z, R, alpha, alpha>epsilon)) + + +## +## Impedance of a capacitor +## (all SI units) + +C = ufloat(10.74e-9, 1e-9) # capacitor +R = ufloat(996, 4) # resistor + +# frequency (Hz) +nu = array(60, 75, 100, 150, 200, 300, 400, 500, 600, 700, 800, 1000, 1200, + 1300, 1400, 1500, 1800, 2500, 5000, 10e3, 50e3, 100e3, 150e3, + 200e3, 300e3, 350e3) + +# V input +Va = array(10.2, 10.2, 10.2, 10.2, 10.2, 10.2, 10.2, 10.2, 10.2, 10.2, 10.2, + 10.2, 10.2, 10.2, 10.2, 10.2, 10.2, 10.2, 10.2, 10.2, 10.2, 10.2, + 10.2, 10.2, 10.2, 10.2)/2 + +# V output +Vb = array(45.0e-3, 55.2e-3, 74e-3, 106e-3, 140e-3, 212e-3, 276e-3, + 344e-3, 460e-3, 520e-3, 580e-3, 700e-3, 860e-3, 920e-3, + 1.05, 1.16, 1.41, 1.88, 3.60, 6.15, 9.89, 10.09, 10.16, + 10.16, 10.17, 10.17)/2 + +# time offset Va - Vb +Oab = array(4.123, 3.300, 2.480, 1.650, 1.230, 0.820, 0.611, 0.490, + 0.404, 0.344, 0.302, 0.240, 0.196, 0.181, 0.167, 0.155, + 0.128, 0.089, 0.040, 0.016, 0.95e-3, 0.23e-3, 0.10e-3, + 0.05e-3, 0.025e-3, 0.020e-3)*(-1e-3) + +# time offset I - V +Oiv = array(4.200e-3, 3.320e-3, 2.500e-3, 1.645e-3, 1.250e-3, 830e-6, 635e-6, + 498e-6, 420e-6, 350e-6, 310e-6, 250e-6, 208e-6, 192e-6, + 176e-6, 166e-6, 138e-6, 100e-6, 50e-6, 25e-6, 4.91e-6, + 2.49e-6, 1.68e-6, 1230e-9, 826e-9, 720e-9)*(-1) + +# estimate uncertainties +Evb, Eva = ufloat(Vb[0], 1e-3), ufloat(Va[0], 0.1) +sigmaZ = (R*(Eva-Evb)/Evb).s +sigmaH1 = ((Eva-Evb)/Evb).s +sigmaH2 = ((Evb-Eva)/Evb).s + +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 + +# plot and fit Z(ν) +plt.figure(4) +plt.clf() + +# magnitude +plt.subplot(2, 1, 1) +plt.title('impedance (capacitive circuit)') +plt.ylabel('magnitude (kΩ)') +plt.semilogx(nu, abs(Z)/1e3, 'o', color="#36913d", markersize=4.5) + +# fit Y=kX where Y=|Z|, X=1/ν, k=1/2πC +k = simple_linear(1/nu, abs(Z), sigmaZ) +Co = 1/(2*np.pi*k) +f = lambda x: k.n/x + +x = np.arange(nu.min()-10, nu.max(), 10) +plt.semilogx(x, f(x)/1e3, color='#589f22') + +# phase +plt.subplot(2, 1, 2) +plt.xlabel('frequency (Hz)') +plt.ylabel('phase (rad)') +plt.ylim(-1.8, -1) +plt.semilogx(nu, Fiv, 'o', color="#36913d", markersize=4.5) +plt.semilogx(x, x/x * -np.pi/2, color="#589f22") +plt.show() + +phi = sample(Fiv).val() +alpha = check_measures(phi, ufloat(-np.pi/2, 0)) + +print(mformat(''' +φ: {} rad +-π/2: {:.4} rad + +compatibility test: + α={:.2f}, α>ε: {} +''', phi, -np.pi/2, alpha, alpha>epsilon)) + +alpha = check_measures(C, Co) +beta = chi_squared_fit(nu, abs(Z), f, sigmaZ) + +print(mformat(''' +k: {} +C: {} F +Cₒ: {} F + +compatibility test: + α={:.2f}, α>ε: {} + +χ² test: + β={:.2f}, β>ε: {} +''', k, C, Co, + 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.ylabel('amplitude (Vout-Vin / Vout)') +plt.semilogx(nu, abs(H1), 'o', color="#9b2e83", markersize=4.5) + +# fit Y=kX where Y=|H1|, X=1/ν, k=1/2πRC +k = simple_linear(1/nu, abs(H1), sigmaH1) +RCo = 1/(2*np.pi*k) +f = lambda x: k.n/x + +x = np.arange(nu.min()-10, nu.max(), 10) +plt.semilogx(x, f(x), color='#9b2e83') + +# phase +plt.subplot(2,1,2) +plt.xlabel('frequency (Hz)') +plt.ylabel('phase (rad)') +plt.semilogx(nu, Fab, 'o', color="#3a44ad", markersize=4.5) +plt.semilogx(x, -np.pi/2+np.arctan(2*np.pi*x*R.n*C.n)) +plt.show() + +alpha = check_measures(R*C, RCo) +beta = chi_squared_fit(nu, abs(H1), f, sigmaH1) + +print(mformat(''' +k: {} Hz +RC: {} s +RCₒ: {} s + +compatibility test: + α={:.2f}, α>ε: {} + +χ² test: + β={:.2f}, β>ε: {} +''', k, R*C, RCo, + 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.ylabel('amplitude (Vout / Vin)') +plt.semilogx(nu, abs(H2), 'o', color="#9b2e83", markersize=4.5) + +# fit Y=a+bX where Y=1/|H₂|², X=1/ν², a=1, b=1/(2πRC)² +a,b = linear(1/nu**2, 1/abs(H2)**2, 0.01) +RCo = 1/(2*np.pi*um.sqrt(b)) +f = lambda x: 1/np.sqrt(1 + b.n/x**2) + +x = np.arange(nu.min()-10, nu.max(), 10) +plt.semilogx(x, f(x), color='#9b2e83') + +# phase +plt.subplot(2,1,2) +plt.xlabel('frequency (Hz)') +plt.ylabel('phase (rad)') +plt.semilogx(nu, Fab, 'o', color="#3a44ad", markersize=4.5) +plt.semilogx(x, -np.pi/2+np.arctan(2*np.pi*x*R.n*C.n)) +plt.show() + +alpha = check_measures(R*C, RCo) +beta = chi_squared_fit(nu, abs(H2), f, sigmaH2) + +print(mformat(''' +b: {} +RC: {} s +RCₒ: {} s + +compatibility test: + α={:.2f}, α>ε: {} + +χ² test: + β={:.2f}, β>ε: {} +''', b, R*C, RCo, + alpha, alpha>epsilon, + beta, beta>epsilon)) diff --git a/circuits/C3-?.py b/circuits/C3-?.py new file mode 100644 index 0000000..ae43d97 --- /dev/null +++ b/circuits/C3-?.py @@ -0,0 +1,188 @@ +# 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) + +R = ufloat(996, 4) # resistor +L = ufloat(0.014561, 0.000009) +C = ufloat(10e-12, 1e-12) + +# frequency +nu = array( 2e3, 10e3, 25e3, 50e3, 75e3, 100e3, 125e3, 150e3, 175e3, 200e3, + 225e3, 250e3, 260e3, 270e3, 275e3, 285e3, 300e3, 350e3, 375e3) + +# V input +Va = array(10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, + 10.0, 10.0, 8.2, 10.0, 10.0, 8.20, 10.0, 10.0, 10.0)/2 + +# V output +Vb = array( 10.0, 9.4, 7.8, 5.0, 3.28, 2.20, 1.44, 0.96, 0.56, 0.32, 0.144, + 0.085, 0.0396, 0.069, 75.4e-3, 0.144, 0.194, 0.320, 0.360)/2 + +# time offset Va - Vb +Oab = array( 12e-6, 6.0e-6, 5.2e-6, 4.2e-6, 3.5e-6, 3.08e-6, 2.64e-6, 2.36e-6, + 2.08e-6, 1.88e-6, 1.71e-6, 260e-9, 260e-9, 280e-9, 280e-9, 300e-9, + 250e-9, 220e-9, 180e-9) + +# time offset I - V +Oiv = array( 110e-6, 26.0e-6, 11.2e-6, 5.9e-6, 4.2e-6, 3.36e-6, 2.80e-6, 2.40e-6, + 2.12e-6, 1.90e-6, 1.74e-6, 280e-9, 260e-9, 270e-9, 280e-9, 280e-9, + 260e-9, 220e-9, 184e-9) + + +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], 1e-3), ufloat(Va[0], 0.1) +sigmaF = (om[0]*ufloat(Oiv[0], 3e-5)).s +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.ylabel('magnitude (kΩ)') +plt.semilogx(nu, abs(Z)/1e3, 'o', color="#36913d", markersize=4.5) + +# fit Y=kX where Y=|Z|, X=ν, k=2πL +k = simple_linear(nu, abs(Z), sigmaZ) +L = k/(2*np.pi) +f = lambda x: k.n*x + +x = np.arange(nu.min()-10, nu.max(), 10) +plt.semilogx(x, f(x)/1e3, color='#589f22') + +# phase +plt.subplot(2, 1, 2) +plt.xlabel('frequency (Hz)') +plt.ylabel('phase (rad)') +plt.semilogx(nu, Fiv, 'o', color="#36913d", markersize=4.5) +plt.show() + + +alpha = chi_squared_fit(nu, abs(Z), f, sigmaZ) + +print(mformat(''' +k: {} +L: {} H + +χ² test: + α={:.2f}, α>ε: {} +''', k, L, + alpha, alpha>epsilon)) + + + +# plot, fit H₁(ν) +plt.figure(5) +plt.clf() + +# magnitude +plt.subplot(2,1,1) +plt.title('transfer function 1') +plt.ylabel('magnitude (Vout-Vin / Vout)') +plt.semilogx(nu, abs(H1), 'o', color="#9b2e83", markersize=4.5) + +# fit Y=kX where Y=|H1|, X=ν, k=2πL/R +k = simple_linear(nu, abs(H1), sigmaH1) +L = k*R/(2*np.pi) +f = lambda x: k.n/x + +x = np.arange(nu.min()-10, nu.max(), 10) +plt.semilogx(x, f(x), color='#9b2e83') + +# phase +plt.subplot(2,1,2) +plt.xlabel('frequency (Hz)') +plt.ylabel('phase (rad)') +plt.semilogx(nu, Fab, 'o', color="#3a44ad", markersize=4.5) +plt.show() + +alpha = check_measures(R*C, RCo) +beta = chi_squared_fit(nu, abs(H1), f, sigmaH1) + +print(mformat(''' +k: {} Hz +RC: {} s +RCₒ: {} s + +compatibility test: + α={:.2f}, α>ε: {} + +χ² test: + β={:.2f}, β>ε: {} +''', k, R*C, RCo, + alpha, alpha>epsilon, + beta, beta>epsilon)) + + +# plot, fit H₂(ν) +plt.figure(6) +plt.clf() + +# magnitude +plt.subplot(2,1,1) +plt.title('transfer function 2') +plt.ylabel('magnitude (Vout / Vin)') +plt.semilogx(nu, abs(H2), 'o', color="#9b2e83", markersize=4.5) + +# fit Y=a+bX where Y=1/|H₂|², X=1/ν², a=1, b=1/(2πRC)² +a,b = linear(1/nu**2, 1/abs(H2)**2, 0.01) +RCo = 1/(2*np.pi*um.sqrt(b)) +f = lambda x: 1/np.sqrt(1 + b.n/x**2) # magnitude +g = lambda x: -np.pi/2 + np.arctan(2*np.pi*x*R.n*C.n) # phase + +x = np.arange(nu.min()-10, nu.max(), 10) +plt.semilogx(x, f(x), color='#9b2e83') + +# phase +plt.subplot(2,1,2) +plt.xlabel('frequency (Hz)') +plt.ylabel('phase (rad)') +plt.semilogx(nu, Fab, 'o', color="#3a44ad", markersize=4.5) +plt.semilogx(x, g(x)) +plt.show() + +alpha = check_measures(R*C, RCo) +beta = chi_squared_fit(nu, abs(H2), f, sigmaH2) +gamma = chi_squared_fit(nu, Fab, g, sigmaF) + +print(mformat(''' +b: {} +RC: {} s +RCₒ: {} s + +compatibility test: + α={:.2f}, α>ε: {} + +χ² test (magnitude): + β={:.2f}, β>ε: {} + +χ² test (phase): +γ={:.2f}, γ>ε: {} +''', b, R*C, RCo, + alpha, alpha>epsilon, + beta, beta>epsilon, + gamma, gamma>epsilon)) + diff --git a/circuits/C3-RC1.py b/circuits/C3-RC1.py new file mode 100644 index 0000000..69c62ab --- /dev/null +++ b/circuits/C3-RC1.py @@ -0,0 +1,205 @@ +# 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 a capacitor (I) +## (all SI units) + +C = ufloat(10.78e-9, 1e-9) # capacitor +R = ufloat(996, 4) # resistor + +# frequency +nu = array(60, 150, 300, 500, 800, 1e3, 2e3, 3e3, 5e3, + 10e3, 25e3, 50e3, 100e3, 125e3, 175e3, 250e3, 300e3) + +# V input +Va = array(10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, + 10.0, 9.80, 9.60, 9.60, 9.44, 9.60, 9.44, 9.44)/2 + +# V output +Vb = array(42.0e-3, 102e-3, 198e-3, 336e-3, 540e-3, 680e-3, 1.36, 2.0, 3.2, + 5.48, 8.20, 9.0, 9.2, 9.28, 9.40, 9.28, 9.28)/2 + +# time offset Va - Vb +Oab = array(4.200e-3, 1.660e-3, 820e-6, 490e-6, 300e-6, 240e-6, 114e-6, 73e-6, 39e-6, + 15.60e-6, 3.4e-6, 900e-9, 240e-9, 130e-9, 90e-9, 40e-9, 20e-9)*(-1) + +# time offset I - V +Oiv = array(4.200e-3, 1.645e-3, 830e-6, 498e-6, 310e-6, 250e-6, 124e-6, 84e-6, 50e-6, + 25.30e-6, 10.01e-6, 5.0e-6, 2.49e-6, 2.01e-6, 1.426e-6, 997e-9, 826e-9)*(-1) + + +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], 1e-3), ufloat(Va[0], 0.1) +sigmaF = (om[0]*ufloat(Oiv[0], 3e-5)).s +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 (RC circuit)') +plt.ylabel('magnitude (kΩ)') +plt.semilogx(nu, abs(Z)/1e3, 'o', color="#36913d", markersize=4.5) + +# fit Y=kX where Y=|Z|, X=1/ν, k=1/2πC +k = simple_linear(1/nu, abs(Z), sigmaZ) +Co = 1/(2*np.pi*k) +f = lambda x: k.n/x + +x = np.arange(nu.min()-10, nu.max(), 10) +plt.semilogx(x, f(x)/1e3, color='#589f22') + +# phase +plt.subplot(2, 1, 2) +plt.xlabel('frequency (Hz)') +plt.ylabel('phase (rad)') +plt.ylim(-1.8,-1) +plt.semilogx(nu, Fiv, 'o', color="#36913d", markersize=4.5) +plt.semilogx(x, x/x * -np.pi/2, color="#589f22") +plt.show() + +phi = sample(Fiv).val() +alpha = check_measures(phi, ufloat(-np.pi/2, 0)) + +print(mformat(''' +φ: {} rad +-π/2: {:.4} rad + +compatibility test: + α={:.2f}, α>ε: {} +''', phi, -np.pi/2, alpha, alpha>epsilon)) + +alpha = check_measures(C, Co) +beta = chi_squared_fit(nu, abs(Z), f, sigmaZ) + +print(mformat(''' +k: {} +C: {} F +Cₒ: {} F + +compatibility test: + α={:.2f}, α>ε: {} + +χ² test: + β={:.2f}, β>ε: {} +''', k, C, Co, + alpha, alpha>epsilon, + beta, beta>epsilon)) + + +# plot, fit H₁(ν) +plt.figure(5) +plt.clf() + +# magnitude +plt.subplot(2,1,1) +plt.title('transfer function 1') +plt.ylabel('magnitude (Vout-Vin / Vout)') +plt.semilogx(nu, abs(H1), 'o', color="#9b2e83", markersize=4.5) + +# fit Y=kX where Y=|H1|, X=1/ν, k=1/2πRC +k = simple_linear(1/nu, abs(H1), sigmaH1) +RCo = 1/(2*np.pi*k) +f = lambda x: k.n/x + +x = np.arange(nu.min()-10, nu.max(), 10) +plt.semilogx(x, f(x), color='#9b2e83') + +# phase +plt.subplot(2,1,2) +plt.xlabel('frequency (Hz)') +plt.ylabel('phase (rad)') +plt.semilogx(nu, Fab, 'o', color="#3a44ad", markersize=4.5) +plt.semilogx(x, -np.pi/2+np.arctan(2*np.pi*x*R.n*C.n)) +plt.show() + +alpha = check_measures(R*C, RCo) +beta = chi_squared_fit(nu, abs(H1), f, sigmaH1) + +print(mformat(''' +k: {} Hz +RC: {} s +RCₒ: {} s + +compatibility test: + α={:.2f}, α>ε: {} + +χ² test: + β={:.2f}, β>ε: {} +''', k, R*C, RCo, + alpha, alpha>epsilon, + beta, beta>epsilon)) + + +# plot, fit H₂(ν) +plt.figure(6) +plt.clf() + +# magnitude +plt.subplot(2,1,1) +plt.title('transfer function 2') +plt.ylabel('magnitude (Vout / Vin)') +plt.semilogx(nu, abs(H2), 'o', color="#9b2e83", markersize=4.5) + +# fit Y=a+bX where Y=1/|H₂|², X=1/ν², a=1, b=1/(2πRC)² +a,b = linear(1/nu**2, 1/abs(H2)**2, 0.01) +RCo = 1/(2*np.pi*um.sqrt(b)) +f = lambda x: 1/np.sqrt(1 + b.n/x**2) # magnitude +g = lambda x: -np.pi/2 + np.arctan(2*np.pi*x*R.n*C.n) # phase + +x = np.arange(nu.min()-10, nu.max(), 10) +plt.semilogx(x, f(x), color='#9b2e83') + +# phase +plt.subplot(2,1,2) +plt.xlabel('frequency (Hz)') +plt.ylabel('phase (rad)') +plt.semilogx(nu, Fab, 'o', color="#3a44ad", markersize=4.5) +plt.semilogx(x, -np.pi/2+np.arctan(2*np.pi*x*R.n*C.n)) +plt.show() + +alpha = check_measures(R*C, RCo) +beta = chi_squared_fit(nu, abs(H2), f, sigmaH2) +gamma = chi_squared_fit(nu, Fab, g, sigmaF) + +print(mformat(''' +b: {} +RC: {} s +RCₒ: {} s + +compatibility test: + α={:.2f}, α>ε: {} + +χ² test (magnitude): + β={:.2f}, β>ε: {} + +χ² test (phase): +γ={:.2f}, γ>ε: {} +''', b, R*C, RCo, + alpha, alpha>epsilon, + beta, beta>epsilon, + gamma, gamma>epsilon)) + + diff --git a/circuits/C3-RC2.py b/circuits/C3-RC2.py new file mode 100644 index 0000000..787fdd4 --- /dev/null +++ b/circuits/C3-RC2.py @@ -0,0 +1,204 @@ +# 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 * + +## +## Impedance of a capacitor (II) +## (all SI units) + +C = ufloat(8.43e-7, 1e-8) # capacitor +R = ufloat(996, 4) # resistor + +# frequency +nu = array(60, 150, 200, 300, 400, 500, 700, 800, + 1e3, 1.5e3, 2e3, 2.5e3, 3e3, 4e3, 5e3, 10e3) + +# V input +Va = array(10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, + 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 9.92)/2 + +# V output +Vb = array(2.40, 4.32, 5.02, 6.12, 6.86, 7.30, 7.90, 8.28, + 8.55, 9.24, 9.30, 9.40, 9.60, 9.62, 9.65, 9.70)/2 + +# time offset Va - Vb +Oab = array(3.4e-3, 980e-6, 602e-6, 295e-6, 176e-6, 116e-6, 63e-6, 46e-6, + 31e-6, 12e-6, 8e-6, 4.8e-6, 4.1e-6, 2.1e-6, 1.2e-6, 0.4e-6)*(-1) + +# time offset I - V +Oiv = array(4.2e-3, 1.71e-3, 1.26e-3, 820e-6, 630e-6, 500e-6, 348e-6, 318e-6, + 248e-6, 164e-6, 125e-6, 101e-6, 83e-6, 62e-6, 50e-6, 25e-6)*(-1) + + +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], 1e-3), ufloat(Va[0], 0.1) +sigmaF = (om[0]*ufloat(Oiv[0], 3e-5)).s +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 (RC circuit)') +plt.ylabel('magnitude (kΩ)') +plt.semilogx(nu, abs(Z)/1e3, 'o', color="#36913d", markersize=4.5) + +# fit Y=kX where Y=|Z|, X=1/ν, k=1/2πC +k = simple_linear(1/nu, abs(Z), sigmaZ) +Co = 1/(2*np.pi*k) +f = lambda x: k.n/x + +x = np.arange(nu.min()-10, nu.max(), 10) +plt.semilogx(x, f(x)/1e3, color='#589f22') + +# phase +plt.subplot(2, 1, 2) +plt.xlabel('frequency (Hz)') +plt.ylabel('phase (rad)') +plt.ylim(-1.8,-1) +plt.semilogx(nu, Fiv, 'o', color="#36913d", markersize=4.5) +plt.semilogx(x, x/x * -np.pi/2, color="#589f22") +plt.show() + +phi = sample(Fiv).val() +alpha = check_measures(phi, ufloat(-np.pi/2, 0)) + +print(mformat(''' +φ: {} rad +-π/2: {:.4} rad + +compatibility test: + α={:.2f}, α>ε: {} +''', phi, -np.pi/2, alpha, alpha>epsilon)) + +alpha = check_measures(C, Co) +beta = chi_squared_fit(nu, abs(Z), f, sigmaZ) + +print(mformat(''' +k: {} +C: {} F +Cₒ: {} F + +compatibility test: + α={:.2f}, α>ε: {} + +χ² test: + β={:.2f}, β>ε: {} +''', k, C, Co, + alpha, alpha>epsilon, + beta, beta>epsilon)) + + +# plot, fit H₁(ν) +plt.figure(5) +plt.clf() + +# magnitude +plt.subplot(2,1,1) +plt.title('transfer function 1') +plt.ylabel('magnitude (Vout-Vin / Vout)') +plt.semilogx(nu, abs(H1), 'o', color="#9b2e83", markersize=4.5) + +# fit Y=kX where Y=|H1|, X=1/ν, k=1/2πRC +k = simple_linear(1/nu, abs(H1), sigmaH1) +RCo = 1/(2*np.pi*k) +f = lambda x: k.n/x + +x = np.arange(nu.min()-10, nu.max(), 10) +plt.semilogx(x, f(x), color='#9b2e83') + +# phase +plt.subplot(2,1,2) +plt.xlabel('frequency (Hz)') +plt.ylabel('phase (rad)') +plt.semilogx(nu, Fab, 'o', color="#3a44ad", markersize=4.5) +plt.semilogx(x, -np.pi/2+np.arctan(2*np.pi*x*R.n*Co.n)) +plt.show() + +alpha = check_measures(R*C, RCo) +beta = chi_squared_fit(nu, abs(H1), f, sigmaH1) + +print(mformat(''' +k: {} Hz +RC: {} s +RCₒ: {} s + +compatibility test: + α={:.2f}, α>ε: {} + +χ² test: + β={:.2f}, β>ε: {} +''', k, R*C, RCo, + alpha, alpha>epsilon, + beta, beta>epsilon)) + + +# plot, fit H₂(ν) +plt.figure(6) +plt.clf() + +# magnitude +plt.subplot(2,1,1) +plt.title('transfer function 2') +plt.ylabel('magnitude (Vout / Vin)') +plt.semilogx(nu, abs(H2), 'o', color="#9b2e83", markersize=4.5) + +# fit Y=a+bX where Y=1/|H₂|², X=1/ν², a=1, b=1/(2πRC)² +a,b = linear(1/nu**2, 1/abs(H2)**2, 0.01) +RCo = 1/(2*np.pi*um.sqrt(b)) +f = lambda x: 1/np.sqrt(1 + b.n/x**2) # magnitude +g = lambda x: -np.pi/2 + np.arctan(2*np.pi*x*R.n*C.n) # phase + +x = np.arange(nu.min()-10, nu.max(), 10) +plt.semilogx(x, f(x), color='#9b2e83') + +# phase +plt.subplot(2,1,2) +plt.xlabel('frequency (Hz)') +plt.ylabel('phase (rad)') +plt.semilogx(nu, Fab, 'o', color="#3a44ad", markersize=4.5) +plt.semilogx(x, g(x)) +plt.show() + +alpha = check_measures(R*C, RCo) +beta = chi_squared_fit(nu, abs(H2), f, sigmaH2) +gamma = chi_squared_fit(nu, Fab, g, sigmaF) + +print(mformat(''' +b: {} +RC: {} s +RCₒ: {} s + +compatibility test: + α={:.2f}, α>ε: {} + +χ² test (magnitude): + β={:.2f}, β>ε: {} + +χ² test (phase): +γ={:.2f}, γ>ε: {} +''', b, R*C, RCo, + alpha, alpha>epsilon, + beta, beta>epsilon, + gamma, gamma>epsilon)) + diff --git a/circuits/C3-RL1.py b/circuits/C3-RL1.py new file mode 100644 index 0000000..382d9c0 --- /dev/null +++ b/circuits/C3-RL1.py @@ -0,0 +1,264 @@ +# 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() diff --git a/circuits/C3-RL2.py b/circuits/C3-RL2.py new file mode 100644 index 0000000..44c18c5 --- /dev/null +++ b/circuits/C3-RL2.py @@ -0,0 +1,264 @@ +# 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() diff --git a/circuits/C4-1C.py b/circuits/C4-1C.py new file mode 100644 index 0000000..2f076ad --- /dev/null +++ b/circuits/C4-1C.py @@ -0,0 +1,263 @@ +# 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 I (charge) +## + +## +## overdamped oscillation +## + +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 +V0 = ufloat(10, 0) # generator peak tension +nu = ufloat(2, 0) # generator frequency + +# time +t = array(0.2, 0.4, 0.6, 0.8, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, + 12, 13, 14, 15, 16, 17, 18, 19, 20)*1e-3 + +# tension +V = array(9.8, 9.56, 9.32, 9.09, 8.89, 7.87, 6.97, 6.17, 5.47, + 4.85, 4.3, 3.81, 3.37, 2.99, 2.65, 2.34, 2.09, 1.83, + 1.64, 1.46, 1.27, 1.14, 1.01, 0.9) + +# fit V(t) +w0 = 1/um.sqrt(L*C) +y = R/(2*L) +a = V0 +b = um.sqrt(y**2 - w0**2) + +def over(t, a, b, y): + return a * y/b * (np.exp(-(y-b)*t) - np.exp(-(y+b)*t)) + +x = np.linspace(t.min(), t.max(), 500) +f, p = curve(t, V, over, 0.01, guess=[V0.n, b.n, y.n]) + +# χ² and compatibility tests +alpha = chi_squared_fit(t, V, f, sigma=0.01, s=3) +beta = [check_measures(*i) for i in zip([a, b, y], p)] + +print(mformat(''' +overdamped oscillation fit +α: {}V, β₀={:.2f} +β: {}Hz, β₁={:.2f} +γ: {}Hz, β₂={:.2f} + +χ²: α={:.2f}, α>ε: {} +''', p[0], beta[0] + , p[1], beta[1] + , p[2], beta[2] + , alpha, alpha>epsilon)) + +# plot V(t) +plt.figure(1) +plt.title('overdamped oscillations (charge)') +plt.xlabel('time (ms)') +plt.ylabel('tension (V)') +plt.scatter(t*1e3, V, color='#5c6652') +plt.plot(x*1e3, f(x), color='#718062') +plt.show() + + +## +## underdamped oscillations +## + +R = ufloat(100, 4) # resistance +L = ufloat(3.82e-3, 0.5e-03) # inductance +C = ufloat(44.5e-9, 0.5e-9) # capacitance +V0 = ufloat(5, 0) # generator peak tension +nu = ufloat(1, 0) # generator frequency + +# time +t = np.append(np.arange(0, 100, 2), np.arange(104, 250, 4))*1e-6 + +# tension +V = array(0.285, 0.522, 0.712, 0.879, 1.040, 1.158, 1.229, 1.286, + 1.292, 1.256, 1.214, 1.117, 1.026, 0.900, 0.777, 0.605, + 0.437, 0.271, 0.121, -0.045, -0.196, -0.335, -0.453, -0.564, + -0.636, -0.701, -0.758, -0.777, -0.770, -0.766, -0.715, -0.673, + -0.602, -0.524, -0.450, -0.351, -0.258, -0.142, -0.062, 0.041, + 0.122, 0.211, 0.268, 0.341, 0.382, 0.426, 0.454, 0.466, + 0.471, 0.453, 0.366, 0.259, 0.135, 0.017, -0.081, -0.165, + -0.233, -0.270, -0.287, -0.249, -0.217, -0.164, -0.078, -0.003, + 0.052, 0.109, 0.146, 0.147, 0.159, 0.146, 0.122, 0.100, + 0.050, -0.002, -0.034, -0.070, -0.073, -0.100, -0.115, -0.093, + -0.061, -0.054, -0.023, -0.014, 0.023, 0.051, 0.059) + + +# fit V(t) +a = R*C*V0 +y = R/(2*L) +w0 = 1/um.sqrt(L*C) +w = um.sqrt(w0**2 - y**2) + +def under(t, a, w, y): + return a * np.exp(-y*t) * (y*np.cos(w*t) + w*np.sin(w*t)) + +x = np.linspace(0, t.max(), 500) +f, p = curve(t, V, under, 0.01, guess=[a.n, w.n, y.n]) + +# χ² and compatibility tests +alpha = chi_squared_fit(t, V, f, sigma=0.01, s=3) +beta = [check_measures(*i) for i in zip([a, w, y], p)] + +print(mformat(''' +underdamped oscillation fit +α: {}Wb, β₀={:.2f} +ω: {}Hz, β₁={:.2f} +γ: {}Hz, β₂={:.2f} + +χ²: α={:.2f}, α>ε: {} +''', p[0], beta[0] + , p[1], beta[1] + , p[2], beta[2] + , alpha, alpha>epsilon)) + +plt.figure(2) +plt.title('underdamped oscillations (charge)') +plt.xlabel('time (μs)') +plt.ylabel('tension (V)') +plt.scatter(t*1e6, V, s=12, color='#a17e3e') +plt.plot(x*1e6, f(x), color='#c8b491') +plt.show() + + +## +## critical dampening +## + +R = ufloat(586, 4) # resistance +L = ufloat(3.82e-3, 0.5e-3) # inductance +C = ufloat(4.45e-8, 0.05e-8) # capacitance +V0 = ufloat(2.5, 0) # generator peak tension + +# time +t = array(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 14, 18, 22, 26, 30, 34, + 38, 42, 46, 50, 54, 58, 62, 66, 70, 74, 78, 82, 86, 90, + 94, 98, 102, 106, 110, 114, 118, 122, 126, 130, 134, 138, + 142, 146, 150)*1e-6 +# tension +V = array(0.007, 0.355, 0.654, 0.907, 1.120, 1.300, 1.437, 1.560, + 1.639, 1.726, 1.772, 1.821, 1.726, 1.557, 1.355, 1.145, + 0.954, 0.783, 0.642, 0.526, 0.411, 0.328, 0.254, 0.213, + 0.157, 0.125, 0.095, 0.069, 0.057, 0.040, 0.034, 0.026, + 0.025, 0.024, 0.013, 0.014, 0.003, 0.010, 0.010, 0.007, + -0.001, 0.002, -0.004, -0.006, -0.008, 0.002) + + +## fit as critical + +y = R/(2*L) +a = R*C*V0 + +def critic(t, a, y): + return a * y**2 * t*np.exp(-y*t) + +x = np.linspace(0, t.max(), 500) +f, p = curve(t, V, critic, 0.005, guess=[a.n, y.n]) + +# χ² and compatibility tests +alpha = chi_squared_fit(t, V, f, sigma=0.005, s=2) +beta = [check_measures(*i) for i in zip([a, y], p)] + +print(mformat(''' +critically damped oscillation fit +α: {}Wb, β₀={:.2f} +γ: {}Hz, β₂={:.2f} + +χ²: α={:.2f}, α>ε: {} +''', p[0], beta[0] + , p[1], beta[1] + , alpha, alpha>epsilon)) + +# plot V(t) +plt.figure(3) +plt.subplot(3, 1, 1) +plt.title('critically damped oscillations (charge)') +plt.xlabel('time (μs)') +plt.ylabel('tension (V)') +plt.scatter(t*1e6, V, color='#43454f') +plt.plot(x*1e6, f(x), color='#65788f') + + +## fit as overdamped + +w0 = 1/um.sqrt(L*C) +y = R/(2*L) +a = V0 +b = um.sqrt(y**2 - w0**2) + +def over(t, a, b, y): + return a * y/b * (np.exp(-(y-b)*t) - np.exp(-(y+b)*t)) + +x = np.linspace(t.min(), t.max(), 500) +f, p = curve(t, V, over, 0.01, guess=[V0.n, b.n, y.n]) + +# χ² and compatibility tests +alpha = chi_squared_fit(t, V, f, sigma=0.01, s=3) +beta = [check_measures(*i) for i in zip([a, b, y], p)] + +print(mformat(''' +critically damped oscillation fit (as overdamped) +α: {}V, β₀={:.2f} +β: {}Hz, β₁={:.2f} +γ: {}Hz, β₂={:.2f} + +χ²: α={:.2f}, α>ε: {} +''', p[0], beta[0] + , p[1], beta[1] + , p[2], beta[2] + , alpha, alpha>epsilon)) + +# plot V(t) +plt.subplot(3, 1, 2) +plt.xlabel('time (ms)') +plt.ylabel('tension (V)') +plt.scatter(t*1e3, V, color='#5c6652') +plt.plot(x*1e3, f(x), color='#718062') + + +## fit as underdamped + +a = R*C*V0 +y = R/(2*L) +w0 = 1/um.sqrt(L*C) +w = um.sqrt(abs(w0**2 - y**2)) + +def under(t, a, w, y): + return a * np.exp(-y*t) * (y*np.cos(w*t) + w*np.sin(w*t)) + +x = np.linspace(0, t.max(), 500) +f, p = curve(t, V, under, 0.01, guess=[a.n, w.n, y.n]) + +# χ² and compatibility tests +alpha = chi_squared_fit(t, V, f, sigma=0.01, s=3) +beta = [check_measures(*i) for i in zip([a, w, y], p)] + +print(mformat(''' +critically damped oscillation fit (as underdamped) +α: {}Wb, β₀={:.2f} +ω: {}Hz, β₁={:.2f} +γ: {}Hz, β₂={:.2f} + +χ²: α={:.2f}, α>ε: {} +''', p[0], beta[0] + , p[1], beta[1] + , p[2], beta[2] + , alpha, alpha>epsilon)) + +plt.subplot(3,1,3) +plt.xlabel('time (μs)') +plt.ylabel('tension (V)') +plt.scatter(t*1e6, V, s=12, color='#a17e3e') +plt.plot(x*1e6, f(x), color='#c8b491') +plt.show() diff --git a/circuits/C4-1D.py b/circuits/C4-1D.py new file mode 100644 index 0000000..744b793 --- /dev/null +++ b/circuits/C4-1D.py @@ -0,0 +1,258 @@ +# 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 I (discharge) +## + +## +## overdamped oscillation +## + +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 +V0 = ufloat(5, 0) # generator peak tension +nu = ufloat(2, 0) # generator frequency + +# time +t = array(0, 40e-6, 80e-6, 120e-6, 160e-6, 200e-6, 240e-6, 280e-6, 320e-6, + 360e-6, 400e-6, 440e-6, 480e-6, 520e-6, 560e-6, 600e-6, 640e-6, 680e-6, + 720e-6, 760e-6, 800e-6, 840e-6, 880e-6, 920e-6, 960e-6, 1e-3, 1.04e-3, + 1.08e-3, 1.12e-3, 1.16e-3, 1.24e-3, 1.32e-3, 1.48e-3, 1.64e-3, 1.88e-3, 2.12e-3, + 2.32e-3, 2.68e-3, 3.00e-3, 3.56e-3) + +# tension +V = array(0.0, -9.0, -8.6, -8.2, -7.9, -7.6, -7.2, -6.9, -6.6, -6.3, -6.0, -5.7, -5.5, + -5.2, -5.0, -4.8, -4.6, -4.4, -4.2, -4.0, -3.8, -3.6, -3.5, -3.3, -3.2, -3.0, + -2.9, -2.8, -2.7, -2.5, -2.3, -2.1, -1.8, -1.5, -1.1, -0.8, -0.7, -0.5, -0.3, + -0.18) + +# fit V(t) +w0 = 1/um.sqrt(L*C) +y = R/(2*L) +a = V0 +b = um.sqrt(y**2 - w0**2) + +def over(t, a, b, y): + return a * y/b * (np.exp(-(y+b)*t) - np.exp(-(y-b)*t)) + +x = np.linspace(t.min(), t.max(), 500) +f, p = curve(t, V, over, 0.01, guess=[V0.n, b.n, y.n]) + +# χ² and compatibility tests +alpha = chi_squared_fit(t, V, f, sigma=0.01, s=3) +beta = [check_measures(*i) for i in zip([a, b, y], p)] + +print(mformat(''' +overdamped oscillation fit +α: {}V, β₀={:.2f} +β: {}Hz, β₁={:.2f} +γ: {}Hz, β₂={:.2f} + +χ²: α={:.2f}, α>ε: {} +''', p[0], beta[0] + , p[1], beta[1] + , p[2], beta[2] + , alpha, alpha>epsilon)) + +# plot V(t) +plt.figure(1) +plt.title('overdamped oscillations (discharge)') +plt.xlabel('time (ms)') +plt.ylabel('tension (V)') +plt.scatter(t*1e3, V, color='#8c4f4a') +plt.plot(x*1e3, f(x), color='#845336') +plt.show() + + +## +## underdamped oscillations +## + +R = ufloat(200, 8) # resistance +L = ufloat(3.82e-3, 0.5e-03) # inductance +C = ufloat(44.5e-9, 0.5e-9) # capacitance +Vo = ufloat(5, 0) # generator peak tension +nu = ufloat(1, 0) # generator frequency + +# time +t = np.append(np.arange(0, 100, 2), np.arange(104, 250, 4))*1e-6 + +# tension +V = array(-0.318, -0.571, -0.761, -0.939, -1.088, -1.210, -1.290, -1.348, -1.379, + -1.350, -1.299, -1.253, -1.130, -0.993, -0.873, -0.721, -0.558, -0.391, + -0.238, -0.081, 0.075, 0.221, 0.359, 0.473, 0.573, 0.650, 0.726, + 0.760, 0.789, 0.790, 0.778, 0.740, 0.684, 0.646, 0.561, 0.484, + 0.390, 0.301, 0.201, 0.117, 0.023, -0.060, -0.149, -0.223, -0.304, + -0.349, -0.387, -0.431, -0.437, -0.464, -0.421, -0.355, -0.260, -0.160, + -0.052, 0.048, 0.148, 0.206, 0.249, 0.257, 0.266, 0.231, 0.168, + 0.109, 0.060, -0.009, -0.068, -0.098, -0.141, -0.147, -0.154, -0.138, + -0.110, -0.089, -0.043, 0.002, 0.029, 0.056, 0.073, 0.079, 0.094, + 0.084, 0.064, 0.055, 0.026, 0.001, -0.014) + +# fit V(t) +a = R*C*V0 +y = R/(2*L) +w0 = 1/um.sqrt(L*C) +w = um.sqrt(w0**2 - y**2) + +def under(t, a, w, y): + return a * np.exp(-y*t) * (- y*np.cos(w*t) - w*np.sin(w*t)) + +x = np.linspace(0, t.max(), 500) +f, p = curve(t, V, under, 0.01, guess=[a.n, w.n, y.n]) + +# χ² and compatibility tests +alpha = chi_squared_fit(t, V, f, sigma=0.01, s=3) +beta = [check_measures(*i) for i in zip([a, w, y], p)] + +print(mformat(''' +underdamped oscillation fit +α: {}Wb, β₀={:.2f} +ω: {}Hz, β₁={:.2f} +γ: {}Hz, β₂={:.2f} + +χ²: α={:.2f}, α>ε: {} +''', p[0], beta[0] + , p[1], beta[1] + , p[2], beta[2] + , alpha, alpha>epsilon)) + +plt.figure(2) +plt.title('underdamped oscillations (discharge)') +plt.xlabel('time (μs)') +plt.ylabel('tension (V)') +plt.scatter(t*1e6, V, s=12, color='#a17e3e') +plt.plot(x*1e6, f(x), color='#daba8b') +plt.show() + + +## critical dampening +R = ufloat(586, 20) # resistance +L = ufloat(3.82e-3, 0.5e-3) # inductance +C = ufloat(4.45e-8, 0.05e-8) # capacitance +V0 = ufloat(5, 0) # generator peak tension + +# time +t = array(0, 10, 20, 30, 40, 60, 80, 100, 120, 160, 220, 280, + 340, 420, 500, 560, 620, 700, 780, 860, 940, 1040, 1120, 1180, + 1330)*1e-7 + +# tension +V = -array(0.18, 0.91, 1.54, 2.03, 2.52, 3.11, 3.58, 3.70, 3.77, 3.59, 3.04, 2.37, + 1.76, 1.10, 0.72, 0.49, 0.30, 0.19, 0.07, 0.08, 0.05, 0.02, 0.01, 0.03, + 0.01) + + +## fit as critical + +y = R/(2*L) +a = R*C*V0 + +def critic(t, a, y): + return -a * y**2 * t*np.exp(-y*t) + +x = np.linspace(0, t.max(), 500) +f, p = curve(t, V, critic, 0.01, guess=[a.n, y.n]) + +# χ² and compatibility tests +alpha = chi_squared_fit(t, V, f, sigma=0.01, s=2) +beta = [check_measures(*i) for i in zip([a, y], p)] + +print(mformat(''' +critically damped oscillation fit +α: {}Wb, β₀={:.2f} +γ: {}Hz, β₂={:.2f} + +χ²: α={:.2f}, α>ε: {} +''', p[0], beta[0] + , p[1], beta[1] + , alpha, alpha>epsilon)) + +# plot V(t) +plt.figure(3) +plt.subplot(3,1,1) +plt.title('critically damped oscillations (discharge)') +plt.xlabel('time (μs)') +plt.ylabel('tension (V)') +plt.scatter(t*1e6, V, color='#5c6652') +plt.plot(x*1e6, f(x), color='#718062') + + +## fit as overdamped + +w0 = 1/um.sqrt(L*C) +y = R/(2*L) +a = V0 +b = um.sqrt(y**2 - w0**2) + +def over(t, a, b, y): + return a * y/b * (np.exp(-(y-b)*t) - np.exp(-(y+b)*t)) + +x = np.linspace(t.min(), t.max(), 500) +f, p = curve(t, V, over, 0.01, guess=[V0.n, b.n, y.n]) + +# χ² and compatibility tests +alpha = chi_squared_fit(t, V, f, sigma=0.01, s=3) +beta = [check_measures(*i) for i in zip([a, b, y], p)] + +print(mformat(''' +critically damped oscillation fit (as overdamped) +α: {}V, β₀={:.2f} +β: {}Hz, β₁={:.2f} +γ: {}Hz, β₂={:.2f} + +χ²: α={:.2f}, α>ε: {} +''', p[0], beta[0] + , p[1], beta[1] + , p[2], beta[2] + , alpha, alpha>epsilon)) + +# plot V(t) +plt.subplot(3,1,2) +plt.xlabel('time (ms)') +plt.ylabel('tension (V)') +plt.scatter(t*1e3, V, color='#8c4f4a') +plt.plot(x*1e3, f(x), color='#845336') + +## fit as underdamped + +a = R*C*V0 +y = R/(2*L) +w0 = 1/um.sqrt(L*C) +w = um.sqrt(abs(w0**2 - y**2)) + +def under(t, a, w, y): + return a * np.exp(-y*t) * (-y*np.cos(w*t) - w*np.sin(w*t)) + +x = np.linspace(0, t.max(), 500) +f, p = curve(t, V, under, 0.01, guess=[a.n, w.n, y.n]) + +# χ² and compatibility tests +alpha = chi_squared_fit(t, V, f, sigma=0.01, s=3) +beta = [check_measures(*i) for i in zip([a, w, y], p)] + +print(mformat(''' +critically damped oscillation fit (as underdamped) +α: {}Wb, β₀={:.2f} +ω: {}Hz, β₁={:.2f} +γ: {}Hz, β₂={:.2f} + +χ²: α={:.2f}, α>ε: {} +''', p[0], beta[0] + , p[1], beta[1] + , p[2], beta[2] + , alpha, alpha>epsilon)) + +plt.subplot(3,1,3) +plt.xlabel('time (μs)') +plt.ylabel('tension (V)') +plt.scatter(t*1e6, V, s=12, color='#a17e3e') +plt.plot(x*1e6, f(x), color='#daba8b') +plt.show() diff --git a/circuits/C4-2.py b/circuits/C4-2.py new file mode 100644 index 0000000..543a08b --- /dev/null +++ b/circuits/C4-2.py @@ -0,0 +1,314 @@ +# 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() diff --git a/circuits/__pycache__/lab.cpython-36.pyc b/circuits/__pycache__/lab.cpython-36.pyc new file mode 100644 index 0000000000000000000000000000000000000000..97a82900dd2a15edb6ede7c46a5c1743e0575d5f GIT binary patch literal 10552 zcmb_iOKcoRdhXXe`4COhdfRT>vc$0^Qk48!L6&97ma~CPJ1culb-bQ!PW6yO^*mJF zLyGh;NK7Hp23QBd9u_ecxo=LfPYIHf1LSC*a_LI|2RZndLoRzszVGjz9+HymT_8hr zb#+yBb=CiW|Fd|0VxnpqhVkW-hVh@c=RXIv%ebN#-O5mAXtYhwY+IgX^4SXQw&OYN zf>)4dJ1n+KUPP2VxpbZ}(t$sOBs)uU<#Jx{3@)nn@59n*VS9a4wU@{Bs7j^choy`+w-$IMRn| zSxH&7)>~1l+exPSomQi(f_m7B1MPKd#wHvEO7ZaW=HyOCD)-wGKi*4pe>J5Bov1h|#59H{5Kt;s#g1pOszj;vQ8| z#KxLAFm_C=)7&~$b=~W$T^+Cb9pyH=+FdxKTJ502Mf;)aYwd4h3|4D(dPzZ^@!n&UnN6H4Ug%xwhl1;pJ-zNHcQS5hOH}8+L&DnkIuGk-L4@}UB4XVn<(qbB))-xzha3QkNf;{Id z<%-#Ul>_c0^$S?QVz=EQS%qJ2x9N8FP&;4NF1qVSF`GBpN#D#k4Mv-{_pL1{Wn)hB z^BE-%qyBrY?;O_m!H{n8<{z@Uauln&3f7~VsSn|c`&U+oS}Se;{H}E!;kt0;Kd$#O zuIM6)!0=3Ec$P9fTUnl?Y~?@{3#y=sxEEDPm2od+YNVX0k&2oC9ww7wZI;75HOMCbyi(TTFv?zdBpv(>-NR{!?i*1*T$ z?9-T|!;yjF$UwFHPX7A3uhEzd!~^;o0S1Esqjq8(+l2B9)>~M4m(yWH8=ttZi>r=T z1FoFIz9S=Y#8B>diR_(4B{l}8GLIOsxn}KH`gm-o&pYO}89Qr*fw^O@6{%jftpTOe zTGSMaS}CbQ5sLUVf}}uZMBc<(K?lk(?&|mPoH!wftS2x$DfyA`$#d3Q!3K1n-3TLn zlwD@DA6}(;y%TKI>pw7=UtALRbrn~pQwE|Ixhv|FMSOA+ z>agv{^^B*9LlPwp2?_agq9ANOU?tc=y)MdDL~!>Sq^lv#_%81>)v|`aLR^b{KRLvy<3>H46W>rZV)3DM~#XcdAVd z>fF+CVuM!tS&XPz`lP(yU0&Nyv*~(WSY2;YHMX{bdi^ID8y!Spm}S#3tN0tw{-(@H zbL&XX*8Ok({v7VVhv@y^X_f^n|<>eSA2}o@%_o+Eff6vb2 z-CJMz*Z-OQ=&!%}uh;Y}mQgE=2$7hvehzKlLy*LYvGo9hGQ+wxmQjT5i?wP}h6wv@ zU?Zv4>+P=Uhpdm+>$m!TnDsypz~`vfy&O_K-RkOAyjid7r@4YtEEZUh1DY0T7i;Lm zXfg9tG#7A*ZYYmCrel?hWxG64E{YOrPcFY;$ z7*ygRT*^8I@dmR#!(h;Rkvb-HufA~&4~=dcMkh$Y2t%{%FSpuW7+|(T-$G?A~hY7Ah$di<%^CP{;AddjXe_ik}+TupcN% ztQ}~#0np6aVbOoXCu-yVKVtYEV!6@KX+?8NmIqVhK6*f6=?#j#14;1uLx}=ec{oG- zi%S=mDfRG~-g5k`+D?=}>d@k8l&ZC8nR4%loNk_O>t}BFT-C zV;(m3tEla<7dici<``C(`klI)GrDbN8z+Gms6bB~3YYAnGYy zRIs80KoD+06BS?yuYV66ZeLvb_#EkwBFEOR;ES5f)?6(#KRv?4RQtDUH&$DJg{~o3 zUPFD)W;us|g-$=D3flxcp!NU>tG|ZE;l>uur1JAaTZ2RNhiwj9yL4MNCrZuKTx7N- zeF^pbn=)_x3a!H}UFHZ1R4xRx4t91$Ww8|2Oy?enOdee+#L><9=E7xC+x*_x!0d%V$ zbQ*A$G5*S)@yl=q!H=-jWrPuEuor51P31+(B${N2ZY4(rEl(j`{|#5^<~bk z+>+AMVmQV~Dxw}+0hE4ra@)KIP=Ns^F6$n^b;Mw)j6=o|V+VHiZ>(Dm^-pTjY>~&* zw$Ec0#>8`H6Z0l~Me|BxE+*zuiU%sg-&HTto#mo0{XIX-RTPb_Cjt9<1goV3HFq)W z!W>7JP#KO*#i*28%o^Si^$TCD)odmOX!n%>epSD%FQczk(Qo3GnD1&zihfES-Vtab zd|FD(Hxv8nYD<;{?Ja}e;7GnEwk{dlwV>aiH5x}@ROvxit)u2~OTUBX28EeFTDx2_ z-^atiRwmrmC*VxNx&facA+ZSr7CMt?wFYo7*UTLv8sP#J`wlT)I3{Q%zJX}CuyAP| zIyE49p`T!CFkQ4%j#nf$knWH`4Sl=t+UlsF7ob3JKQw(9apM9Y_JH>01$S9wpanUA zQ#Xu}z6XJfkt9Hz?e;j_4+%60Y%h%bLdMM<_r~cQ$vT&jJ?|~<0Vg%WbzBjFBjBuV zL#o#@0H)kCbqYHDItC?X%F_|v^ffg9GbhP__RKE6!ncke4z^-y;ZNXYar5;bv33~+ z{CG%ot(2kcMq)No!!JN9$?oen^xGUS-c#;EB?tC>{vR+#k|IX6Y#y{GY2Uc#T;WM$@HtaAIZ&A!esO0B&BHp&k+Hsv28P@QFROw_mV5FC7 zIl%!Hra1H1wP+0nPWmS-#0R{BnwLT0=Dw{b6H$w)7-T)Uf`7-*ln4&_f6&w)pd}{) zb#zXIukdgm5vWvP{Gcu}BFq@9!>tBM#AN0@XvgSbKS!kmoWVG%p6y~f&#*O}XZTh& zrX&O|urDP#-{h>@<+O){&RQxkw=lTy!is~^pt55%smQFYlRvUs&S zx$i>FPXr4_mr=@~?!febvz#{~O)W6bi#Bn?(rzlr@uz|Ig8u8`TOPh`RyGXmQDtwJIKlo6rN*g7Pi-EHU65&4H{8A)cbwM>nRt(`*Z5`;>zo4S1I z_p@YV}B2e zs7fFmM1aAC@yNt*#-UD#|8Z-|dIDL!P=O}5Vkj3V7K>vPv_3U5-vO5F{y$+>!1d~C(6|}pF*@Bj7{h+Z_yQU5jjkI7 zF{@d+mk|kkbcs5UL(fwwf!|nlyNvVaI01siErZLOW%I_7m|=b}ht9vxZQP4v)$MiOWMf&Lk<*{M7XbHbkm&%lo z?1vP`EE$n0juY56Ns)7Dw?4a1hwvqej%YMW-ov8@I;%glklJV=Dd-9*=%}RAs6aqR zb5$!Om0_e6CFY0XlXwops?;yst00r9{rh^F-=fQB6lcsCn8j&JxA2s!1mcg^E|owX z%@~q60s?WwgxJ&2Ps6u>lIU2)xgXwH&Vydy5EnxI5}%Il0wwSX{m1Gy`;3BQ6oJJ z(g)CpgS1M$E+?_b(dX7(-$LJI5`q9QU^%vq-oM7VY{J|_9U_KHjqadeI`pQcufT8^ za=3OBNWsE8)6!Bz-|5&?2p~V@p^Aa18&ttq`(>QZVcZ2z_zgfuivNmCk1d>>2Kh@}%2U2- zdRavr+w#0mTCA<(58(RfW2Wxl?b?`r0XO|33nmVbAKQ$Q@pYsXLBV>f1Kh2l-E+cj zXGO1b908`eiew_gUn2;^AqKh7;i9QkWt*>Zy_9|CFUB2dYv3{e5qA7VB-BtSMdWBW z%bqPF$Bd$#Oudu-r$8=nc@)f9hL>b3Q|A;AJWOARDvy2-&WNe%2sWGqBUJ!n%9CT9 zxF_a&h}FgEW4xwA7Q(FoYRPyihH%e;HY5dwVFhttys@m07#tE{#I--e>mLrOV1P@L z*4E?uH!`PC zx_Z|u1OyJid(d~Qn)=b3&QFQ?Um|n(H?2UE58lL{j`V<(22kQnZ)Q*LTpv>2+w_X* zJ^e4BW70M^3NLq6wM4F)oK8G_{y$YxkQY+(ROF#9;4ij(hT|9kr8cRQCGzEPQ;w_B z1Q(Ax^ko+Bu%MTbI+o1aNOnN~oHe@g;y{W87jxxRa%hXB2qOiGNvWMeBq;^Gs0EoE zL1CmKlJdp$nD-J{7+pg#Wiu%~ZjaL}+Kec&zf!qWF5@aqzg(_Vr^*$SGu4CT1LbOY Tx_Y2`s62-H`1FyfQrY@Hxxk%` literal 0 HcmV?d00001 diff --git a/circuits/lab.py b/circuits/lab.py new file mode 120000 index 0000000..3e8499d --- /dev/null +++ b/circuits/lab.py @@ -0,0 +1 @@ +../lab.py \ No newline at end of file diff --git a/circuits/lab.pyc b/circuits/lab.pyc new file mode 100644 index 0000000000000000000000000000000000000000..4fbfc71297b7c538c8aaa2860fce2b815ec6b4da GIT binary patch literal 12721 zcmd5?No*WfdVbw(HrWzo(h|i*^ir}dtL4#R%G*GZZS9T|5A77~G1VAPr`T0Z7FAu< zsznsXA~3NuW6xwV10<&yIp&hgC7x3P-v4UkKlkl1_-irY-w2*_PT6GtUNXkV|B{J}>6A>PWIAQjD4P;X z%O)Bzor-BxOlQE~xnW4524 zGmZIv{*Y-LGJeh6F{ZI#{5{G%V*I_zJZk(2WgavBJ`>EFjmJ&9Wc-I%c-xq(4afMC zYH`^3W5%B{cR-aVjsLI;pECY_WsVsCfbpkQ`?Ol(^`N-};vF^qjM~(VKda0!|G{6qj9fmB#1C(TS(wE66kdurY$blsjKVxf zy(qKI2%ad}4AM2P)A90Ayed}*dQ3ioBuLx#rHYsJt?ude;_&^QD^cR*b+*e{U$pt2 zSCHLIy>9+6GOZ*{gUHJRzj?!piXdCE{l`$(O`@%Kn75EMTsEIKiQaw4nJixvX_uMwVrf>Q)#^3*+@1{`AO^4D#+Jrd8gj+R-MzQT3*yjqVQBJOj|{m zXQv`>^<;O;oj~75@noFXR9aj$v!u{#+}HtRdK=T;)=K7fxxc+vQTdha6>%#_b1#hT z)||~SUq_=>5@)#==gy!$A8wn|71+0UvSZ^c57JbO1v7SO(}dSkc2z_$*xOc+W9oQ7hPFUZ)!cl1$VX z?(_nAc^;Vn8W^e-X_*G5F+ziK2d@?5k0>*0{E9MFU079S+_=ki(zJLMgA(qV z$;r36fm4rt!PUu%)mG$X*)cSDExZu`{f_VDo@3YLBt_mW^33UY@szVw=7Lj)quy!r{=OgDJjJ#Z~9)0awmu zAz35|qb&20TDgo2)XFVzrDW#wlG!LzucXJQ#CxSsquV7DqhZBtjItCZ8`YkM-7cG3 zR74zqDMVt#q?{H4Rz%9C6-bkkQ&~3fQttfNn*o?s1CGD#zo6R`bY&Ku#u#_BAnkgM9&6NgPiL*;9^;9toj5RZsWY*Fw z*_Begz{J2ZHhaj?R>%uw$TX{PXi~BuYVVR*Qa+l^EKi%wH#xEdv((}bAI_#d{2Cdu z*O+VLCY@qt*>c~q7Yw%m0IRXlW`zVFFrOu(Dgg=v2bNS-GPa{+OeNzxN@^kXV!c6#?XP#!eUeRR_pD^v^1#M{Bxf|TcO@P(jPrp-Pk2Q?#GbL zri7r`rj-kS5Hy=EO0yXxrZiEilqT_SfA8O1>EJhnYVj7vfo8CJGV3b~syWBYhZ^O@ zZy#|kzTtS81AR3#s`A*J`pR)<`P3I5eH>Ha!h52WxJ01{YKcUyh)BCb99P7_m|wQ2 z}Yp zQ0{@eW*UMscruE7b)s6TmTT48NNr!O+9F^1aH#^t31k2Y*g4P?2!P_6gXID|EHI#1 zrWcUR;(=mXSQeIxpQ9VNyF_gsITtTo!AmRYz=;T~%)oJ2^;ScikG7n*>Nn4N$DB1e z9xF}=nQ>zEjJBKsI07GIq_E`){Ri0k)hMutqJeNe`&&qwKNP@$wRb+heG)wd6WQ~6 zg?l1kB#ka#t-^a2IV+3w`LIXmwZffcDV_&Eut*Db*Td@=%!*;P#o76AJ)-!F$lMw+ z>C+NM_!A=f?Ge373TZq>%-k}Cbw@$()b`ite;4VUtbc$*u_f}?)AuXXg_`2|N8S(! zjPCa&G35N-?DfvBoL!}=!};;8`U;u7xgNmJN}W}*TZ@Duv%f8hGvBj5r|8 zD>p?P?%V?EJaZxrN#dU4E1q2tXwGr{PmuKJVGj(HXL{xI<|WgAjts~F0$3mcEamNr z*~Ym5XE&-;$n!Q_L8(;;H5zP>nKiVpunAr-y|-w#;S|EBL@$h`qr?UrN5F90%-bJ~ zH9LPBR{Hy8bN!szhUch00P6IuF+4Rh2Q_^FzeMadJjj~41<$f-10&>c@o!vu55UmZ zx#B*avko(ztfd~&DYoHv<&1L=T{2M7D`@uS8K-~LWW+(b9ln|Gc&k|w!6*>;*NY%- z;UtOn7w&4mit`X~5QQB^(^#LOesIcEY$ip>1_xA#6z98-@K6R{`TWB#KKS9MzZ$^A z3KH36eFjRsJy(&He%NlyUcmaS1uo-8)=yhR&U&^*$Gw1;y3S$*#fo&*n@TUZ4jZy~ zN$Gb4LD?7IHoT_t%Np$BdbmsI`1Wb6!d7cZ$VGURWh%UW1*+VVNKQT5hs8&sUQ0;)$5ti_DXABXoqpvkn zD4zq+768%%1RoAM0I`lf#Bs&Y%j`K2|Aa+=kKVZ1Ys%YaoK@jE1h(Ov_k#lWK^21niip&cE}MH1!J}BtnW2BrXjAN5Jw8BZ zI>7y1{e-urI?3*m*Vg@bk~v~ZColIhb3~OrML}L{hbt^1n64^q%Rmkx7Rt}m%9>Z* zCA8ZmG%*CY15G$#<$r~ujACQ}0!7ez z1Ol9#5kdxP4X@uT8><^34FCm7paEwQ;evBXuog#=fcU*b#8EhucO4-dpgqpRh>;*0 zqbIT4=aGEb38a8zt+2awl0mIZZREm3nFJ~PdkG|zRGmArlK)J9FtLT!T!Nf;>`!-h zsI8^^9hXOtn9-@y6wj(8Dwi_wPu{}|Y%OCdvKl?XgSXiC9H&7*`_hxU1yC3CW-k~M zHS5)&v>7i|^5AtgaWq=GoL+O&N`jzpa!Id@j`MO}yOY&|VByUzja3s9Vcar+g6afd_L zwosjxn_rrW193!p4gCg|ueKPz0j9zKaoYZ#Dy}i#9`OBXkl;7){rtBs3^_lB_QT+j zwURQXUt+&Z?bosYb<&*jdP(Vd_qSQVVXa0xzaxKPm#{&&T2TBm=i$-ZR1TIWN>7wt zEFCW&F3lpZgBsVtrAx2T?E_&Eai^Wc5P+(0T#;vg!{#P|Z5=3#aD1~aw8QH^-mKt24%JOzC7ng)P3gH-as$)iSX@NWQk{;1HBR$vCSM|na7jV) zmq|;$&dW_4!-pJKEHH)5^>SByC$%)_+i79r*;jBK8n`u-E9C7CQE(1r(gfo1W)QBe z=b0R@RuZ(^uFpPrcENAwvWE2d54>el9Jw-4sv$NwQ+@HnQuy5?8bklA7 zOUd}KQ>fL?(;vGcsueYVAF zKD1m1`qz_ocbM#FJTOLbkyGdyRV6S|xo_ZN7rlJ5}yFUym^vLn673z1Yed zk-S6RNLUH1RO?kq+8RC?$XrfsC%eWR1s46K4)VeT0w^1~eF6gQ&MSKL!B?8>73IV~U5 zvK$S}12LM2YLS%vtG#TN9Oukt?AxUlf2+MmcdUPm7t3J)9IEird$YVyrM1-Y5-dEh zKvlMnydzuXsr4JI4U9r?iZ~vUULvb2{!Gu zwhak-PKLmFgaGN_?nob;@H-qK0ViI3Z{!zfwB_X9HC)-_-Y#>Rh!1%$Ywn=J z0^m#T7%30W!kGTyIg_}{H+WYXLpGd;W-0ew0tHPj2u30fVDSWDK}BV%CrU@khf4G1 z78MSW!bjpCAOmh;eUN${VI;s{n1uqk3U)$urTDB1?1pBD70|#*6B|KmT49)4@Eote zmwkn4$2nBtHI7JrbIRg5Bp~b+-A(g$q9^MP#gvR4en)`MC~!p&4~dLGH_M2E1b%<} z(@);t*2za|8D~t%rSxu;tSvt4QbP5bCPnN^H`@zw8xDUG)uO++(M{otQGqHjwmU3w zpTKyF0~S;kN{{fuw5(W*W~(Og%Nvj2Oq2i_KIjkd^l`^ ziU~&Pb$}DdH{|%0AF&WPnX@tj-@udrz;G5+;oF+^omi`Qv z@jX|9M3I@aH*j8&C-|r!wA9K9e6ix*MiaT2H+lq?Q9HtpG`_n;%^zc^bEGD>Wnk3? z*F7u{l>Z(XfQ@*Cm-zccKvA^3OcKB?1mZu?cP6$YwH^e8 zLyVf2%Ix^NJqZD5K&0g>8IM zM9FdT#2K)x*T6=5p(>@^?ltB;dvVBnPiVtxcd*VOXuivjCoiXIlDaznjB`_Q-sIOf z?vL112%^qh3v&B4oE1`HNb-Fx*b*LIj{cSnRuKj+7w!I#i7@bU)%wZFLZ7?H?kulahKbI;n4oUG;((} zI9Wgsb+5|EUxhC#d;3^+Q6(lnwGl+OkjZD!28R3zz_iwB=RvX7N mHvdwsHZfJh-^BFfOl`8Zr?!7`a&oq|7x{_#`I)ha^8W&T_1I1T literal 0 HcmV?d00001 diff --git a/circuits/transfer.py b/circuits/transfer.py new file mode 100644 index 0000000..722b55f --- /dev/null +++ b/circuits/transfer.py @@ -0,0 +1,59 @@ +# coding: utf-8 +from __future__ import print_function + +from sympy import * +from sympy.utilities.lambdify import lambdify +import matplotlib.pyplot as plt +import numpy as np + +## +## Transfer function calculator +## + +I # imaginary unit +i = Symbol('i', positive=True, real=True) # current +R = Symbol('R', positive=True, real=True) # resistance +Rl = Symbol('Rl', positive=True, real=True) # resistance +L = Symbol('L', positive=True, real=True) # inductance +C = Symbol('C', positive=True, real=True) # capacitance +w = Symbol('w', positive=True, real=True) # angular frequency + +# impedances +Zl = Rl+I*w*L +Zc = 1/(I*w*C) +Zr = R + +# input/output tensions +Va = (Zc+Zr+Zl)*i +Vb = (Zc)*i + +# transfer function +H = Vb/Va +phase = simplify(atan2(re(H), im(H))) +magn = simplify(sqrt(re(H)**2 + im(H)**2)) + +print('arg H =', phase) +print(' |H| =', magn) + + +## +## Plot function +## + +p = (9.94e3, 20, 3.8e-3, 4.45e-8) # R,Rl,L,C +x = np.linspace(100, 150e3, 1500) +f = np.vectorize(lambda nu: lambdify((w, R, Rl, L, C), magn)(2*pi*nu, *p)) +g = np.vectorize(lambda nu: lambdify((w, R, Rl, L, C), phase)(2*pi*nu, *p)) + +plt.subplot(2, 1, 1) +plt.title('transfer function') +plt.xscale('log') +plt.ylabel('magnitude') +plt.plot(x, f(x)) + +plt.subplot(2, 1, 2) +plt.xscale('log') +plt.plot(x, g(x)) +plt.xlabel('frequency (Hz)') +plt.ylabel('phase (rad)') +plt.show() diff --git a/lab.py b/lab.py new file mode 100644 index 0000000..a900d1f --- /dev/null +++ b/lab.py @@ -0,0 +1,351 @@ +# coding: utf-8 + +from __future__ import division, unicode_literals + +# special functions +from sympy.functions.special.gamma_functions import lowergamma, gamma +from sympy.functions.special.error_functions import erf +from mpmath import hyp2f1 + +# uncertainty propagation +from uncertainties.core import UFloat +from uncertainties import ufloat, wrap, correlated_values +import uncertainties.unumpy as unp + +# fit +from numpy.polynomial.polynomial import polyfit + +import sys +import collections +import types +import numpy as np +import string + + +## +## Variables +## + +def array(*args, **kwargs): + """ + Shorthand for 1-dimensional array + """ + return np.array(args) + + +def uarray(err, *val): + """ + Shorthand for uncertain array with + constant uncertainty + """ + return unp.uarray(val, [err]) + + +def nominal(x): + """ + Nominal value of n-dimensional uncertain array + """ + return unp.nominal_values(x) + + +def sigma(x): + """ + Uncertainty of n-dimensional uncertain array + """ + return unp.std_devs(x) + + +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 val(self): + """ + Gives the measure at standard 68% CL as an uncertain float (X, S/√n) + """ + return ufloat(self.mean, self.stdm) + + def tval(self): + """ + Gives the measure at 68% CL, calculated with the Student + t-distribution, as an uncertain float (X, tS/√n) + """ + return ufloat(self.mean, t(0.68, self.n-1) * self.stdm) + + +## +## Normal distribution +## + +def phi(x, mu, sigma): + """ + Normal CDF + computes the probability P(xx) + """ + 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 + dX = np.diff(B) + + E = [N*p(X.mean, X.std, B[k], B[k+1]) for k in range(M)] + Chi = sum((N*O*dX - E)**2/E) + + 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) + + +## +## Student t-distribution +## + +def t(p, n): + """ + Student's t-distribution quantile + """ + + from scipy.stats import t + dist = t(n-1) + tm, tp = dist.interval(p) + return (abs(tm) + tp)/2 + + +## +## 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 + """ + if sigma_Y is None: + weights = None + elif 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) + return correlated_values(coeff, cov) + + +def curve(X, Y, f, sigmaY=None, guess=None, **kwargs): + """ + Any function fit + """ + from scipy.optimize import curve_fit + + if sigmaY is None: + weights = None + elif isinstance(sigmaY, collections.Iterable): + weights = np.asarray(sigmaY) + else: + weights = np.repeat(sigmaY, len(X)) + + coeff, cov = curve_fit(f, X, Y, sigma=weights, p0=guess, **kwargs) + if sigmaY is None: + return lambda x: f(x, *coeff), [ufloat(i, np.nan) for i in coeff] + return lambda x: f(x, *coeff), 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) + 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) + + +class WithEmptyFormatter(string.Formatter): + int_type = (int, long) if sys.version_info < (3,) else (int) + + def vformat(self, *args): + self._automatic = None + return super(WithEmptyFormatter, self).vformat(*args) + + def get_value(self, key, args, kwargs): + if key == '': + if self._automatic is None: + self._automatic = 0 + elif self._automatic == -1: + raise ValueError("cannot switch from manual field specification " + "to automatic field numbering") + key = self._automatic + self._automatic += 1 + elif isinstance(key, int_type): + if self._automatic is None: + self._automatic = -1 + elif self._automatic != -1: + raise ValueError("cannot switch from automatic field numbering " + "to manual field specification") + return super(WithEmptyFormatter, self).get_value(key, args, kwargs) + + +class MeasureFormatter(WithEmptyFormatter): + def format_field(self, value, format_spec): + if isinstance(value, UFloat): + return value.format(format_spec+'P') + else: + return super(MeasureFormatter, self).format_field(value, format_spec) + +mformat = MeasureFormatter().format + +epsilon = 0.05 diff --git a/optics/__pycache__/lab.cpython-36.pyc b/optics/__pycache__/lab.cpython-36.pyc new file mode 100644 index 0000000000000000000000000000000000000000..3ac8d04331e094457835322fb6b3472d748fd38a GIT binary patch literal 10550 zcmb_iOKcoRdhXXe`4COhdfRT>vc$0^Qk48!L6&97ma~CPJ1culb-bQ!PW6yO^*mJF zLyGh;NK9d+4X_S^9D*2|+&8D#rv%B#0dllYx%4G~gB*O!A(y=*-}iS<4@pV(E|4L* zy1J^my6XSG|5-dgF;V^7#o)%5Pa4L58oPcD>X&gvv1u5}P-bYfP0wswo@Mgc3hlPz zIqia1kY_tAwo6_~>P}d0SG>6myovUtH`zYm9gwzCIJNeOH;oy~ z;iK&tZzgMf%zF%NrSPD4Q0|AkLvlat9g+J{@2IM%>K(&7rpDAb?#I=Hn#BEabwEww z{={2`dPGg%G1PQ$WbMfv+jG^UY6d+|sTb8_>fjyIds-b*htcwkI--u^enP#Zj;qJf z^Q?MVJ)xdN?WA(mQ|LRTo>tG`KC4csXK}BoSJX*$3L~CVKT@-*hT8M$IrTjHPOH=E z4DM&toSMgd4r@897O<9u;ON>s)^Zjt=hS($oDa_8c_CQXG1ks`=hX||3#xS8s9pGu z?5J5uS+&+%QLEcYruv;$qpN~?*op(~hY@_#?LO%{vecccHL9~d`rCvAOY_{SCSHPc@UGL%^ zRZ+yonmI6bOsvz~I#qSu>#JQIulgP3HoMwgIHOwapu6p=X#rZ67N3F)p4nYBhL*$SAJ}>alM5r;z(IOHn3vuSh9uJ?%P7Y(+G6z zw>s(a-Hlefib0KTCyM<}?B@M(wl%wty%qbz&4C3vu|ZSWR$5Hs(|QKQ2`)rdT9D^_ zrCc%FuX4bBgnj`FSnRf2B&u-B?Ka)c9$M$i+C_E!C}#5}JL#JlpTTDH_P(_xrEJV; zemZI;75HOMCbyi(TTFv?zdBpv(>-NR{!?i*1*T$ z?9-T|!;yjF$UwFHPX7A3uhEzd!~^;o0S1Esqjq8(+l1;1)>~M4m(yWH8=ttZi>r=T z1FoFIz9S=Y#8B>diR_(4B{l}8GLIOsxn}KH`gm-o&pYO}89Qr*fw^O@6{%ddtpTOe zTGSMaS}CbQ4T|_Rf}}uZMBc<(K?kZZ?&|mPoH!wftS2x$DfyA`$#d3Q!3H#+-3TLn zlwD@DA6}(;y%TKI>pw7=UtALRbrn~pQwE|Gxhv|FMSOA+ z%CPOn^^B*9LlPwp2?_agq9ANOU?tc=y)LR&L~!>Sq^lv#_%81>)Ut+OpaM!7H_Lil z8P_u%Csega(}7X9q4%7u?w~HTMtcGEV%At@wtgbhy8e3aIz5Wx7jSiwP%(Cg2Rs4--e^ch9 zxpgFG>;5-?cMkX8!}I`rD*n+>1}MF*ACCm(rB|MHufFa2kqa>#DQL07?(EWO_xk*| z_df4XkHARUx`JtiD(S{e+4#OaYwoR%D^t^3)oeqs3ww4@#2CTH6R1G;fYBIQ5qIdt z!^RHyy5`K_4$Rv9O$>VYj?#^sc4>n9;xvd^l#O4H`zq+f@^TE51f;f>`_vq#zh~$1 z?yax<>;KGt^cP?K*K2wf%cvDbgh$bFL)ORZ^;>;E%zB^);BnOJUJj|AZgq7l-mKU4(_FzR77HxM0Zogvi#2p% zw3vA+nhUr@H-IGR&0?`JLZgW z3@Y&uE@d5qc!OD=VKC^uNF5WpSKqjXheo#zqZ6cHgrQmXms{`_H{JJUZ(sCl?ur9J44Q#`57^FO==5oLO7eH8zkeyox_nVzLzDdUEV<>7) zQV@V^0x*^TIkwwk-V>|c*R z)(*7W0BC0Iu;{y64Ou2<^9w%ds`R~k>p0n zF%O&iRn+#_i=2K$a||m?{Z8G@atClrpQ+Eqy|?;J48$V7x%<_Z8OVmzlBOIx5cQNU zDp=70APBdhi3+fU*T07jw=XVze2#QTkz?yu@I_5#YpxcWpB`ais{Pxw8>_89N7oQ6 zuc5wYvz$Y~LZ=^6g>3>JPWw8|2Oy?enOdee+#L><9=E7xC+x*_x!0d%V$ zbQ*A$G5*S)@yl=q!H=-jW;o3L2>rPuEuor51P31+(B${N2ZY4(rEl(j`{|#4^<~bk z+>+AMVmQV~Dxw}+0hE4ra@)KIP=Ns^F6$n^b;Mw)j6=o|V+VHiFRfb+^-pTjY>~&* zw$Ec0M#OVx6Z0l~Me|BxE+*zuiU%sg-&QZuo#mo0{XIX-RTPb_Cjt9<1goV3HFq)W z!W>7JP#KO*#i*28%o^Si^$TCD)odmOX!n%>epSD%FQczk(Qo3GnD1&zihfES-Vtab zd|FD(Hxv8nYD<;{?Ja}e;7GnEwk{dlwV+?4H5x}@ROvxit)u2~OTUBX28EeFTDx2_ z-^atiRwmrmC*VxNx&facA+ZSr7CMt?wFYo7*UTLv8sP#J`wlT)I3{Q%zJX}CuyAP| zIyE49p`T!CFkQ4%j#s2LknE5^4Sl=t+UlsF7ob3JKQw(9apM9Y_JH>01$S9wpanUA zQ#Xu}z6XJfkt9HzW$uEh5duvD+Y2MVka2Uzy>U86vd(2>&wGn|z)6j89alu)2so?T zkm|JzfGPJ(oq|rkjzNi;@^pkZeGSe3%tJ{|Ahbq=->1n+L5)+BfbwSNXJi zBm^-Uy@Z013lo>yNLm6krh-xgSnPofNTK$Twe1L`pbQ6vLjaSwR)sXBm3p+N;7jfC zXI)5MCvJrdn$yVC)aUTA-Oi6U7ByPE%{hksk@R|HD55~a7fozvw^}h(1<&zCzD@Cr zV$X`(xO)3YAtU|`qoOhjqrgyRGW}}qM{@I&FxfmgNog9ub3{ODY}=3PT*@-)A~is>(RQ0XmoG9dXuU(RXQ0C80n>2 zPH;elDb750En0(tll}<{@d2-(=4DX0xo_*qMATv`202fz;NLMcC4xi#A2js`Xvv8{ z9i0>5D?Hps1S%C6Kd6h02r~xjaH~NQF`0P}+A(_A&rvA>XE2VcXS`RHxH#zHeIqe~#vz7|XEetNau;QRJsO(rxDl%*9fkc5fzW&Oj>30xOv(<0lhzsx?VzS z?z>R)6TyPfWmGcFKdHc!&RO_pc9s@;PFb;yWQ8Oa&^Lppf`T15W`H;#;0{C=Z6z`V zqAdwNRrpmAB3*zrfY&l8iXA9!xgei91BX5J-CSn3ndn5&1f_<~MG@Wiuh9x`AWzg? z_nc-w3}G9OW0A0L-7pM9#Yai`$IQJWE*jA}DQcJqKh{EO8b3XTB2#)=&@QeXA)k~s zg4W7v93_QDH)uAs$MO4AUv;#`|3-H-g~D*kW)XrvWzD4G|8unEMzw+mBA*|QXoPd! zSA_1fl??-XRN335&MbUg~w1GxPIls0B zg}=5wgVy||pjYO+L&vgdNnW*_W>xPQG~UE+)>lE!wH#OA83EgZ>FU-&(EW8JMu)aD zV@0M=lS-qn*8>`zML5gU5}wsNo0~kQm#FE*AER^6l{U?L#@|_=(YZGlg`!z988g%| zM*k6uk$ov&4E^@9@-N{i%81ZsY#ox%?zVI3i2TE}j3l$zTBgRu)=nXH2|}gVOUP@DQC%L_%Soj+_o5JBzqzUA7#*bOcvA>5! zR3(rOBEaCncw}NY<4`BW|F|_}J%Oy=aZ?9q&#gc4KqI)<0MdDc38dYJEr97;LwJM$ z$4u``s6dljF_a4wi^VYtTAv!3?*L17|DP}`;Cgj6XxxnQ7@h7MjA1`ye1Q!3M%Rsk znAI%Z%ZLO%xxya&V9)Yrg*C44i{Ly^62Dcy-|nG6E7#;{+=H zW=B#3m@8&tfL~6P2t8SN^?^o+#5eKrrEfp|Gf7CG@+q3cFpGkNJ`%psvQb$YPXGm5NU@Y zL1K57p;JovTFw;V6$pszASsqQb8xzmEK~hPtFz)F4HRso=Ot-E!*?TO*r4%ud+B(W z2cE(V#11&4RmdYEnHKvuaTe6&QZOc$)jgWyDCi-YV|p9GBK>rN^4Kt3w1i)&OJzz( z_CtzemW;?0#|dnkq{z9nTc6#hL--O!M>HBG@8Qt{oz)*&NNu!`6m*3YbX3x5R3MXk2cu~5KnezsHABJwpfME0Ia}w!+G}#S!$fK( z@GfPQ&>9_);e-&UQ}CBDUTm8%wNFX=Wc#ZxKH07n6C06%5Xmcq-Iay2ngWVIz1#1o zq?jHdA@_Fzt*Iq+dOB!eFcTm91XMqRQ42X;oN050haZ--H_(^=?5q~iGZ`M225<{h z5%f436oO*Y>R3U+E8T@@P4h-xS?Z1|+$npNyCAMtg^DaC(;p%u`ewTqZ>Cg%kMteC zum}EZpbIDmGk`q~Fk+lWz~4ZcYG7uT8>!=J-S5X;oSd~9Ng?Xvv{bk8UQ$}msF5B9 z=>urQL0Y9=my=lJ=yU6?Z=vrp2|)lDupC=Q?_c9wHev3e4iUqpMt4v!9eUH!S70~{ zIb1smq+sElX=y2?(JPPJ4!%MQ2sdEYl9;Poul}En^XT(%>1RG9*kt%>O<;gKl z+!OOX#OmVoF<#Rl3*puPwPZXML%8Qa8S3dt6{(kcs59ZC(vf>{Lm|2HK z;|+1f;$rPZ9u7A9G0t$oNlu_iN$*s&*+xER7_CHeB4U0xC;xJVH_{AXG&;6;_$$z6 zUA^lS0s;r%J?OhtP5tOi=cmN{FOfO?lUAU~2XA6eM|wa?11RxkGP9?5t`8~iZQ_5D zwUIy7J;ciUk*Q5?L5sLosDDDLrnF(=6JID6+p&xl}IWDowv!u2iSW6_hj8gXII| VYI(YPpn9l0hWhyQk*QMI`aikJpI-m~ literal 0 HcmV?d00001 diff --git a/optics/interf.py b/optics/interf.py new file mode 100644 index 0000000..7fb70b6 --- /dev/null +++ b/optics/interf.py @@ -0,0 +1,148 @@ +# coding: utf-8 + +from __future__ import division, print_function, unicode_literals +from lab import * +import uncertainties.umath as um +import uncertainties.unumpy as unp +import numpy as np +import matplotlib.pyplot as plt + +# convert micrometer readings to SI unit +length = lambda x,y,z: x*1e-4 + y*1/4*1e-4 + z*1e-6 + +## +## Part I +## + +k = 60 # number of wavefronts +l0 = ufloat(633, 1)*1e-9 # laser wavelength + +# micrometer calibration +d0 = k*l0/2 # expected value +d = length(0,3,10) # measured value +df = d - d0 # offset + +# laser wavelength measure +d = sample(length(*i)-df.n for i in [(0,3,9.5), (0,3,10), + (0,3,10.5), (0,3,10)]).tval() +l1 = 2*d/k + +print(mformat(''' +# Part I (Fabry-Perot interferometer) +offset: {} μm +λ: {} nm +''', df*1e6, l1*1e9)) + + +## +## Part II +## + +df = length(2,0,1) # micrometer offset +k = 70 # number of wavefronts + +# laser wavelength measure +d = sample(length(*i)-df for i in [(2,0,22.5), (2,0,23.5), (2,0,23.5), + (2,0,23.5), (2,0,23), (2,0,23)]).tval() +l2 = 2*d/k + +print(mformat(''' +# Part II (Michelson interferometer) +λ: {} nm''', l2*1e9)) + + +## compare wavelengths +alpha = check_measures(l1, l2) +l = combine_measures(l1, l2) + +print(mformat(''' +comparibility test λ₁ - λ₂ + α={:.2f} α>ε: {} +λbest: {} nm +''', alpha, alpha>epsilon, l*1e9)) + + + +## +## Part III +## + +d = 3e-2 # air cavity thickness +Patm = 101.3e3 # atmosferic pressure + +k = array(17, 20, 17, 18, 19) # number of wavefronts +P = array(82e3, 83e3, 84e3, 83e3, 84e3) # gauge pressure + +m = sample(*(k*l.n/(2*d*P))) + +print(mformat(''' +# Part III (air refraction index) +m: {} Pa⁻¹ +''', m.tval())) + + +## +## Part IV +## + +d = ufloat(6, 1)*1e-3 # glass pane thickness +t = ufloat(12, 1)*np.pi/180 # angle +k = sample(162, 112, 124, 122).tval() # number of wavefronts +n = (2*d-k*l) / (2*d-k*l/(1-um.cos(t))) # refraction index + +print(mformat(''' +# Part IV (glass refraction index) +n: {} +''', n)) + + + +## +## Part V +## + +# peak position +S = array(8.9, 10.7, 11.7, 12.7, 13.4, + 14.1, 14.7, 15.6, 16.3, 16.7, + 17.3, 18.0, 18.3, 18.9)*1e-2 + +# uncertainty +sigmaS = 2e-3 + +# peak number +n = np.arange(1, len(S)) + +d = ufloat(1e-3, 0) # slit length +L = ufloat(125.1, 0.5)*1e-2 # screen distance +ti = um.atan(S[0]/L) # incidence angle + +# drop first peak +S = S[1:] + +x = np.linspace(0, 13, 200) +peak = lambda n, t, l: L.n*np.sqrt(1/(np.cos(t)-n*l/d.n)**2 - 1) +f, (tio, lo) = curve(n, S, peak, sigmaS, guess=[ti.n, l.n]) + +plt.xlabel('peak order') +plt.ylabel('distance from center (m)') +plt.errorbar(n, S, sigmaS, color='#604848', linestyle='', linewidth=1.1) +plt.scatter(n, S, color='#a17e3e') +plt.plot(x, f(x), color='#65788f') +plt.show() + +alpha = check_measures(l, lo) +beta = chi_squared_fit(n, S, f, sigmaS) + +print(mformat(''' +# Part V (ruler diffraction) +λₒ: {} nm +θₒ: {} rad + +compatibility test λ - λₒ: + α={:.2f}, α>ε: {} + +χ² test: + α={:.2f}, α>ε: {} +''', lo*1e9, tio + , alpha, alpha>epsilon + , beta, beta>epsilon)) diff --git a/optics/lab.py b/optics/lab.py new file mode 120000 index 0000000..077fd87 --- /dev/null +++ b/optics/lab.py @@ -0,0 +1 @@ +/home/rnhmjoj/doc/bicocca/Lab 2/calcoli/optics/../lab.py \ No newline at end of file diff --git a/optics/lab.pyc b/optics/lab.pyc new file mode 100644 index 0000000000000000000000000000000000000000..928d1a61ff18875b148361eba88741e70897df2e GIT binary patch literal 12657 zcmc&)No*WfdVbYxHrWzo(h|i*^ir}dtL4#R%G*GZ7p)#A9@;6|W2!NpPO+<+EULSz zRf{N&MPOoS#-7P!21rgZa-2&tmv~MIkdp^Ukc)lFC6@pO0wl;O2y)0N`M&S%^MFBylW zCDR@;owBKyO=r~9M{VhdX;(~V%+$y1d)c(dO{Z$=Rr@|_+IviAuc_~~?-kRYFr9s- zzRz?XF!cvaXVTOsZQYn@PnnGeO?|(XF>cxiOlR8Erw7#sP5q$i&~(PsXZrbBQ=jeU z=S+RRpFd>khfGj4w~eVUm|%}G51C-EG7p95BJOYM)d~ydE^SLA;|Tm{FUW31*c!W`a3o zo>DIiJFk}CP)EET((7>(EU4KD6Fj8MNfSJ*%qbH*V!}ga<7qATQ5I^sizaYQxTw;@ zgVNI`Jk9!zXH5NB6Fg?>&vB~d+T&n+^gsBkm5|HFkp$88D2tLfpTbKKN6jP%8|^3$ zQ@@?rW&}?=*$mS)zti#aQM@Ww1A0t8f+S2^_N9WC^{wvd*5dH}?aS@N&ueU#v%YBY zJuf4>nfl%QL1daqnucvZ4}->azg>jclI=f+x^B|mYDIYyS;J-Xd4uTPhfI!LI_4(G z;b4DUP;nH=9@kx7Ptttdj{~=rr0(L0AnJrMcj~uYKTZ8Dw6Xi+fN@v*QiUz$`uub) z)0(etGFwBDQIb7}+^O}X6P`-r_0C4JaVkifr&d9}X460Qy1(k4KGpQw&7>WjO1gQ} z%ucoa)sx*VZvs6Z#*=YcQ*Lt2%#uE@a$5(G>1|6pTXW2Wki@os4_Ovioc_jcY%V86dovtIb|+i^ zFMB+VKG@^ieHP|hUvr1!oHP#_L3lmW-n?m)?8??C=~=Ya@tFEFrv(;3ChY;OfaepbYMl9apNu5NYUbH3`(GD zCO6;ehHfnm1XMRGR-0`<%Z{PJtI_okxOW3T_g%X#H!1RNk!NnlkGEWiKptg~Kp?1_ zCS6EM6lQ18H9q6wFD^Q(VM?NRLO(0gFhdmphf1VTr|4MVvDCl8Qvb%CrJ;|;y{R!q z%$8xmmSGV2@u2&4KSg1$Azo0|4AB_e8;c#ft!4zMr;;QrYkSk}i zh%AzHQHuFUpr`p-Ae!4L*2#Q6-`AGjUTj*Q0RL zS_rPX+0MK{a!Eq&&PLn@;IG1+tYh=h; zW3G;ybc&f}%X!OMFxxU8Y{oVlWfFV9c$SQ+1PBlmSW;2R*p8Ahm5lEwsj6g;=En2@ zn&YtV9+mIqTpN3cEf!Abbz(JlU;l>pIGfPgS*f(;VW3|Kqh;E)G*ZHvAiHX0sEvl%4DnQ+R^B>wI1{hMGc z5%T-HL>i3vMLSB) z3E4-Hp*E9j=CU!&%goF{t?{L7rsXnmxy2TeyA#6h|F}zMt&j&Q_eu0x&WiweQs*2# z2gEK~EvS=^--BTJPKvGl`oC)5{rM;V^0I|4E$WbdLTc_k#X_NNS4k;E=Wv=}k$dDE ztyPhi-w$QzcK|x+@J6GP1Vx*#R>8oTMSdqtyJ?i9QN9J39_D(VWJ13{9P{W85Baa$ zefY>!{%7!H6!pqP#i^94mFh@!U$xRCPx)}60>lYq00`JN&=cH;!kUA%0yHcrpi!n5 zkj&zNLRwfBii@A28#uc}9o}{?T)d2zX3~M<5L#)0v#{!~MmP&?xo_5Pob`{nYjQSL z+z8U*#_HMLatFW&1B{Wvekb%9VBc5Up+ytTg8S*;LdJrj=nX8q``N9N=qZTEhS$p6 z6Tu=m^!REO-P6cfDWsQ&Jwl%q-YiS;-0$H;y05z)UBggT04q(-&4;TIy`M+s=7>q3 zlo-N;5Xo4YFGK>Gn3y)v3qJl1_HTIhnEeD=X7 zzv`jBPjVFvzzi*Ci5rA#X$YBsuh~Kd;70d*5c9^_mG7MVIG`lc1E=RMY5oN=AO^N?K?I!ErNhe=yeU z{4H4K@0HB8SIjm%MfCwlr*DqoshT;c=L7g9TDRdbR?SU#k`){NAb*R0oe33PLztkqzKI5FbZ+veD@F>%HK<$een7FKltQV1CUrD zBD<`gK#8|!D6%q$S}oZLSe>=dV|2)RW{b#KPqyfI=kZe0xr?A!k(PQx>3Pp(Ll!S7 z{k8xoyAs%jS5m8K4m$;hE-HE*l3kAQ2F{`0iv6voQKZ*t~#RG4VmQ16D^^!jqga#Ip{$);ytn z4lr8)L=O;pIMx8dI_MC^WkbKR=QjKk76CR&R}K(}@FYT)ib0@m2Lhuw2)iK?MB)0d zLyIeC0Gh8TZ=G>hh3W{%hBMwj{@(*rJpL>8QDc~q(M>`}QJXVEf1XjL*u8RmfX#G( z`aAjwZA&$h-6N;1)A2ZSM3l~1-X-RUDSKvuwAc=pSwtvZQQDGW90DwoZ>g0vu6Rpm zw<}m;m~RJ`aGuKl3Pl-3#_XAJYNdnDBi3}%l2VNbO|K#Yh67FZ?gKFnhj2gRLMtvI&O7#}yBoCD zLjI16BS6gPlrzP1DM`wO47`$e@d6vm*oo{#Pw3zw_Pxbvkk7vKWNrbp1-01=W<jq^jkkD>iANMx;} zjOZ7bFBAJ!%zuqEr>tI7dd~Y63plLRW#@L}&+iH{2vrM-f95PahMV%i(uDJv^SpDs zbl90iUI#O-f(w^Mqq7IHBBoCJh-=ePYOMGoS5I=wg99zo7v%)oK(DMKr7wPnuz>Fd zkf8Mfc4ZYGdZ3zh;F^prq#9}yjDx0xlLj4#D8j9g5voYc3FC@DALO;?xp6y%iv(;7 zPT;s^-DZ~;d%QWp@f%8;f=Vij$eYq|&Ez_!#UZ$8x0h;k5Uf#}XEXT_NgLM?G<%u! z<7>RO#G!l0Z$uc zyVXp>R?7?62hV-`>Ri;28vlW}Y>MNSC!8uGfHS3s5T=Hvr5<2gE_N?3T>cdq@M~zq zDllSzz?dL#gc*1<@mZ-${m1aA)Uh}3KbQ}U`M{yj0)s!ob+IuYLVn(zliLrUUV^@{my%6q?Jp(Y z!%h`;$|u;#ssf4sI5t}_J_@Lu!_7|LL+YbP7u$jm>C-II>mj`Xdp_q0tNV15)qH5P z4)m`c?d~wyPkBma^->Klzm*WxSKyi?V#^Ys{$F(i`c*pIEG5yd+s zjYO5GN~KnabH{~g4S&HBO(p>2% zLh{d9FS76MNatQb!=6EgS(f{VlZoY~%W6c3{~TI(_=_i;E`^1yl3~-xTDZ-BLIVW91Bnab~LvdE~(XG)!?i zrdx~2S3ZCL4;6$%=6w{&_p-5Xg1{=~73Q8|BHznpP7xw1dqp?pw5}}LpVRU?EyK~k z+y|S9rX~s5zqHF{$#2eU#=hw^`CI8dx;y<7yjcDM+)#lR-<{<}DQ%<?p9!k%sk@PTl4&5t17|K*`GOaKM<_OkrQYisJmrec`tSVGy4;UX| zZ6NW5q=PU70(5;|6-MI^@iknp8E! z2_K1mhzxjz)j{HUWRdWKT^0!7DVPbRmEt2UFdJGSR`3GnN^AqIX?tN{!Ew9-UG^0Q z9j8tO(l{cy%_)m#kbtb0b?40chn}3f6j8Eu_;mn2j=;q`ydtsz-7FgllJ~vsPdK#P;FLk>4mL&t+J&7!O==xh9p)e*7-rysY&|U$Jhzo~JST^A( zJq_>y`MO-)+h_@g#(tur&ri6mXaq@<2qxS4B5m*}GJ*6!16;FQhogjP#L%DOlD+3X zkR&pX_I}MT@&um>M3z=rfv-@!TWBK3@_LWFGG<5Ekw$k{X!#=y^a^RoO&M6T!F3J` z1mnL)20$ZL;TQe^k+`8<@(^)U(8Xo7G4Uh)D$o-xFOdLn1cCPV_2r1|NPPz(;Si(d zg%UgdUQapzbC4?iWq%y|ZNg88cXUC^YlA9Gfg7e@Y6_epGzDtmbO3Brg_i2@O);G` zDK4DsmqSYhw`yTekFfY%yf*zfPH^ju+ag%MR*DNcd>jO|1>Z4j<4&#{Hlr545u((% zdEyQj)@xv6x=@l*Yj+#-uDv+q-6yo+v)fqb5HjE4XONfDG)X-jcgDG?IB)R#8t;c} zD*RAeu7$b%w#^DC5hU5Z8g2giMP#KVzQmYYu<;*c*L&98HY<(?zcDpBH4>v zKUHDHIb5oi_1){Y?U!F)Rw~p0QYsH36H2KB#(-1!xD8*bv0C9TY6w1(s27oeVE{B? z8;E$9c{+zO33UIKK?$q%sXY3jXylhMyuFcHDmeO^HYi01cwDsi1194C&yW)#t@^-G z2S~Pe!-Bvl>gChvsN8ErpbzPJAEf*_Mt`SIN&~h(Sh8~unC-F2aIX919k=@(PDnLb zan?TiGyeU>%Qa5e`--7ee*2`pVr~q|hZmIh z?1ismlXOQfqE=S2(}jv(_eggctD(v7l8x>}C;k!>-}Up26*BFE3Gav&%?kMrIJxo*J2=(=oyWVDH~pb*ws$XKen3 jYIS0&ioc2J$(ib8bx(Ew {}, mode ? "shell" }: + +with pkgs; + +let + matplotlib = pythonPackages.matplotlib.override { + enableGtk3 = true; + }; + +in stdenv.mkDerivation rec { + name = "lab-data"; + source = "."; + + buildInputs = with pythonPackages; [ + matplotlib numpy scipy + sympy uncertainties + ]; + + shellHook = '' + exec fish + ''; +}