commit f0fa70042718d4085764486f9f135a9f30ae9afe Author: rnhmjoj Date: Sun Mar 18 18:02:21 2018 +0100 initial commit 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 0000000..97a8290 Binary files /dev/null and b/circuits/__pycache__/lab.cpython-36.pyc differ 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 0000000..4fbfc71 Binary files /dev/null and b/circuits/lab.pyc differ 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 0000000..3ac8d04 Binary files /dev/null and b/optics/__pycache__/lab.cpython-36.pyc differ 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 0000000..928d1a6 Binary files /dev/null and b/optics/lab.pyc differ diff --git a/optics/light.py b/optics/light.py new file mode 100644 index 0000000..6d59547 --- /dev/null +++ b/optics/light.py @@ -0,0 +1,173 @@ +# coding: utf-8 +from __future__ import division, print_function, unicode_literals +from lab import * +from math import pi, sqrt +import matplotlib.pyplot as plt +import matplotlib.mlab as mlab + +## +## measure focal length +## + +# convert micrometer readings (mm) +length = lambda x,y,z: x+0.5*y+0.01*z + +# object - lens distance (mm) +p = ufloat(4608, 10, tag='p') + +# object size (mm) +h0 = array(10.00, 15.00, 18.60, 20.00, 25.00, 30.00, 35.00, 40.00, 45.00, 50.00) + +# image size (mm) +h1 = array(length(11, 0, 11.5) - length(10, 1, 1.5), + length(11, 0, 39.5) - length(10, 1, 1.5), + length(11, 0, 27.5) - length(10, 0, 20.5), + length(11, 0, 41.0) - length(10, 0, 27.0), + length(11, 0, 34.5) - length( 9, 1, 42.0), + length(11, 0, 46.5) - length( 9, 1, 28.0), + length(11, 1, 3.0) - length( 9, 1, 3.0), + length(11, 1, 32.0) - length( 9, 1, 1.5), + length(11, 0, 23.5) - length( 8, 1, 14.5), + length(12, 0, 5.0) - length( 9, 0, 22.0)) + + +# linear fit +sigma_h = 0.03 +x = np.linspace(5, 55, 10) +k = simple_linear(h0, h1, sigma_h) +k.tag = 'k' +f = ufloat(252, 1) +f0 = p*k/(k+1) + +# tests +alpha = chi_squared_fit(h0, h1, lambda x: k.n*x, sigma_h, s=1) +beta = check_measures(f0, f) + +print(mformat(''' +# focal lenth measure +k: {} +fₒ: {} +χ² test: α={:.2f} +t-test fₒ-f: α={:.2f} +''', k, f0, alpha, beta)) + +plt.figure(1) +plt.title('focal length (f=252mm)') +plt.xlabel('object size (mm)') +plt.ylabel('image size (mm)') +plt.scatter(h0, h1, color='#206a99') +plt.errorbar(h0, h1, sigma_h, color='#1a567c',linestyle='') +plt.plot(x, k.n*x, color='#1a567c') +plt.show() + +f = f0*1e-3 + + +## +## measure Δs +## + +ds = sample( + length(12,0,5.5) - length(11,1,21.5), + length(12,0,4.5) - length(11,1,20.0), + length(12,0,5.5) - length(11,1,21.0), + length(12,0,4.5) - length(11,1,20.5), + length(12,0,6.5) - length(11,1,22.5), + length(12,0,6.5) - length(11,1,21.0), + length(12,0,6.5) - length(11,1,21.0), + length(12,0,3.5) - length(11,1,20.0), + length(12,0,5.5) - length(11,1,21.0), + length(12,0,5.5) - length(11,1,20.5), + length(12,0,6.5) - length(11,1,22.0), + length(12,0,4.5) - length(11,1,21.0), + length(12,0,3.5) - length(11,1,20.0), + length(12,0,5.5) - length(11,1,21.5), + length(12,0,7.5) - length(11,1,22.5), + length(12,0,4.5) - length(11,1,21.0), + length(12,0,4.0) - length(11,1,21.0), + length(12,0,6.5) - length(11,1,22.0), + length(12,0,6.0) - length(11,1,21.0), + length(12,0,6.5) - length(11,1,22.0), + length(12,0,4.5) - length(11,1,20.0), + length(12,0,6.5) - length(11,1,21.5), + length(12,0,6.5) - length(11,1,21.0), + length(12,0,5.5) - length(11,1,22.0), + length(12,0,5.5) - length(11,1,20.5), + length(12,0,5.5) - length(11,1,22.0), + length(12,0,7.5) - length(11,1,22.5), + length(12,0,4.5) - length(11,1,20.0), + length(12,0,7.5) - length(11,1,22.5), + length(12,0,5.5) - length(11,1,21.5), + length(12,0,4.5) - length(11,1,21.0), + length(12,0,5.5) - length(11,1,21.0), + length(12,0,5.5) - length(11,1,20.0), + length(12,0,5.5) - length(11,1,20.0), + length(12,0,8.5) - length(11,1,22.5), + length(12,0,6.5) - length(11,1,20.0), + length(12,0,6.5) - length(11,1,22.0), + length(12,0,5.5) - length(11,1,20.5), + length(12,0,4.5) - length(11,1,20.5), + length(12,0,5.5) - length(11,1,20.0), + length(12,0,6.5) - length(11,1,22.5), + length(12,0,4.5) - length(11,1,20.5), + length(12,0,4.5) - length(11,1,22.0), + length(12,0,7.5) - length(11,1,22.5), + length(12,0,4.0) - length(11,1,20.0), + length(12,0,4.5) - length(11,1,20.0), + length(12,0,5.5) - length(11,1,22.0), + length(12,0,5.5) - length(11,1,21.0), + length(12,0,4.5) - length(11,1,20.5), + length(12,0,7.5) - length(11,1,22.5)) + +## Δs histogram +plt.figure(2) +plt.xlabel('displacement (mm)') +plt.ylabel('frequency') + +# gaussian best fit +x = np.linspace(ds.min-1e-3, ds.max+1e-3, 300) +gauss = mlab.normpdf(x, ds.mean, ds.std) + +# χ² test +O, B, _ = plt.hist(ds, bins=6, normed=1, facecolor='#506d72') +alpha = chi_squared(ds, O, B) +print(''' +# Δs histogram +χ² test: α={:.2f} +'''.format(alpha)) + +# plot +plt.plot(x, gauss, 'black', linewidth=1) +plt.show() + +ds = sample(*1e-3*ds).val() +ds.tag = 'Δs' + + +## +## calculate c +## + +B = ufloat(0.292, 0.001, tag='B') # lens-rotating mirror distance +D = ufloat(5.491, 0.005, tag='D') # rotating-fixed mirrors distance +w = ufloat(3000, 2, tag='ω') # motor angular velocity + +A = (B+D)*f/(B+D-f) +c0 = 8*np.pi * A*D**2 * w / ((D+B)*ds) # speed of light + +print(mformat(''' +# speed of light +c: {} +σc/c: {:.2f}% + +error components:''',c0, c0.s/c0.n*100)) + +for p, sigma in c0.error_components().items(): + print(' {}: {:.2f}%'.format(p.tag, abs(sigma/c0.n)*100)) + +# t-test cₒ-c +c = ufloat(299792458,0) +alpha = check_measures(c0, c) + +print(mformat('t-test cₒ-c: α={:.2f}', alpha)) +exit() diff --git a/optics/micro.py b/optics/micro.py new file mode 100644 index 0000000..3aa2225 --- /dev/null +++ b/optics/micro.py @@ -0,0 +1,339 @@ +# 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 + +from scipy import signal +from scipy import interpolate + + +## +## Reflection +## + +## reflection +i = array(20, 30, 40, 45, 50, 60, 70, 80) # incidency angle (deg) +r = array(321, 300, 282, 270, 261, 241, 223, 203)*(-1) + 360-i # reflection angle (deg) + +# fit and plot θi - θr +x = np.linspace(10, 80, 10) +k = simple_linear(i, r, 1) +f = lambda x: k.n*x +alpha = chi_squared_fit(i, r, f, 1) + +print(mformat(''' +# reflection +k: {} +χ² test: + α={:.3f} +''', k, alpha)) + +plt.figure(1) +plt.clf() +plt.title('reflection') +plt.xlabel('incidency angle (deg)') +plt.ylabel('reflection angle (deg)') +plt.scatter(i, r, color='#003cad') +plt.plot(x, f(x), color='#0069ad') +plt.show() + + +## +## Refraction +## + +# incidency angle +y = ufloat(5.3, 0.1) +x = ufloat(13.4, 0.1) +ti = 180/np.pi * um.atan(y/x) + +# refraction angle +to = ufloat(180,1) - array(170, 168) + ti +n = to/ti + +print(mformat(''' +# refractive index +n₁: {} (styrene?) +n₂: {} (styrene) + +t-test: α={:.2f} +''', n[0], n[1], check_measures(*n))) + + +## +## Brewster angle +## + +# current, amplification factor, angle +data = array( + (0.08, 3, 5), + (0.10, 10, 10), + (0.12, 3, 15), + (0.38, 10, 20), + (0.64, 10, 25), + (0.62, 10, 30), + (0.64, 10, 35), + (0.44, 10, 40), + (0.34, 10, 45), + (0.38, 3, 50), + (0.36, 3, 55), + (0.32, 1, 60), + (0.08, 1, 65)) + +t = data[:,2]*np.pi/180 # angle +I = data[:,0]*data[:,1] # normalized current +sI = 0.02*data[:,1] # normalized error + +# plot Ι - θ +plt.title('Brewster\'s angle') +plt.xlabel('incidency angle (rad)') +plt.ylabel('vertical polarization intensity (mA)') +plt.ylim(0, 8) +plt.scatter(t, I, color='#109d1d') +plt.errorbar(t, I, yerr=sI, linestyle='', color='#0c7616') +plt.show() + + +## +## Polarization - Malus Law +## + +## distance 1 + +d = (104.8 - 20 + 13.5)*1e-2 # transmitter - receiver distance (m) +y = array(0, 15, 30, 45, 60, 90, 105, 120, 135, 150, 180) # receiver angle (deg) +I = array(0.120, 0.100, 0.080, 0.050, 0.020, 0.000, 0.010, 0.020, 0.045, 0.100, 0.120) # current (mA) + +# fit, plot γ - I, with I = I₀cos(γ²) +x = np.arange(0,180,1) +k = simple_linear(np.cos(np.pi/180*y)**2, I, 0.01) +f = lambda x: k.n*np.cos(np.pi/180*x)**2 +alpha = chi_squared_fit(y, I, f, 0.01, s=1) + +print(mformat(''' +# malus law +## at d={}m +I₀: {} +χ² test: + α={:.3f}''', d, k, alpha)) + +plt.figure(2) +plt.title('malus law') +plt.xlabel('horizontal angle (deg)') +plt.ylabel('intensity') +plt.xlim(0,180) +plt.ylim(0,0.2) +plt.scatter(y, I, color='#d6a800') +plt.plot(x, f(x), color='#e9b700') + + +## distance 2 + +d = (78.3 - 34.7 + 13.5)*1e-2 # transmitter - receiver angle (m) +y = array(0, 15, 30, 45, 60, 90, 105, 120, 135, 150, 180) # receiver angle (deg) +I = array(0.175, 0.160, 0.120, 0.065, 0.035, 0.000, 0.010, 0.045, 0.085, 0.130, 0.175) # current (mA) + +k = simple_linear(np.cos(np.pi/180*y)**2, I, 0.01) +f = lambda x: k.n*np.cos(np.pi/180*x)**2 +alpha = chi_squared_fit(y, I, f, 0.01, s=1) + +print(mformat(''' +## at d={}m +I₀: {} +χ² test: + α={:.3f}''', d, k, alpha)) + +plt.scatter(y, I, color='#178200') +plt.plot(x, f(x), color='#1ead00') + + +## distance 3 + +d = (69.2 - 34.7 + 13.5)*1e-2 # transmitter - receiver angle (m) +y = array(0, 15, 30, 45, 60, 90, 105, 120, 135, 150, 180) # receiver angle (deg) +I = array(0.200, 0.180, 0.160, 0.085, 0.045, 0.000, 0.010, 0.050, 0.090, 0.150, 0.200) # current + +k = simple_linear(np.cos(np.pi/180*y)**2, I, 0.01) +f = lambda x: k.n*np.cos(np.pi/180*x)**2 +alpha = chi_squared_fit(y, I, f, 0.01, s=1) + +print(mformat(''' +## at d={}m +I₀: {} +χ² test: + α={:.3f} +''', d, k, alpha)) +plt.scatter(y, I, color='#003cad') +plt.plot(x, f(x), color='#0069ad') +plt.show() + + +## polarizer + +t = array(0, 45, 90) +V = array(0.185, -0.080, -0.04) + + +## +## Standing waves +## + +# receiver position from transmitter (m) +d = np.arange(30, 86.5, 0.5)*1e-2 + +# tension (V) +V = array(0.851, 1.019, 1.971, 0.810, 1.116, 1.591, 0.820, 1.064, 1.448, 0.804, + 1.196, 0.953, 0.772, 1.319, 0.846, 0.773, 1.211, 0.761, 0.806, 1.143, + 0.810, 1.049, 0.629, 0.785, 0.985, 0.571, 0.837, 0.816, 0.563, 0.867, + 0.739, 0.551, 0.839, 0.638, 0.542, 0.796, 0.586, 0.541, 0.828, 0.535, + 0.544, 0.678, 0.556, 0.833, 0.546, 0.537, 0.790, 0.480, 0.511, 0.705, + 0.451, 0.504, 0.637, 0.437, 0.524, 0.574, 0.422, 0.519, 0.525, 0.413, + 0.517, 0.476, 0.396, 0.499, 0.416, 0.379, 0.474, 0.369, 0.371, 0.433, + 0.319, 0.358, 0.383, 0.275, 0.345, 0.335, 0.245, 0.327, 0.297, 0.227, + 0.326, 0.264, 0.224, 0.322, 0.258, 0.235, 0.328, 0.242, 0.253, 0.341, + 0.232, 0.249, 0.315, 0.199, 0.239, 0.264, 0.173, 0.236, 0.227, 0.143, + 0.221, 0.188, 0.125, 0.206, 0.162, 0.112, 0.193, 0.130, 0.109, 0.193, + 0.119, 0.111, 0.172) + +# plot +plt.figure(5) +plt.viridis() +plt.title('standing waves') +plt.xlabel('receiver position (m)') +plt.ylabel('tension (V)') +plt.xlim(0.28,0.87) + +# mark peaks +mins = signal.argrelextrema(V, np.less)[0] +maxs = signal.argrelextrema(V, np.greater)[0] +plt.scatter(d[mins], V[mins], marker='X', color='#d8a300') +plt.scatter(d[maxs], V[maxs], marker='X', color='#1ead00') + +# spline interpolation +x = np.arange(d.min(), d.max(), 5e-4) +plt.plot(x, interpolate.spline(d, V, x), color='#0069ad') +plt.show() + +# calculate wavelength +l1 = 2*np.diff(d[maxs]).mean() +l2 = 2*np.diff(d[mins]).mean() +l = ufloat((l1 + l2)/2, 0.005) + +print(mformat(''' +# standing waves +λ: {} m +''', l)) + + +## +## Double slit +## + + +t = 102+0.05*7 # total width (mm) +s = 20+0.05*9 # slits width (mm) +d = (t-s)*1e-3 # slits spacing (m) + +# maxima order +n = array(1, 2, 3) + +# peaks (right side) +r = array(200, 223, 262)*np.pi/180 - np.pi # 194, 204, 222 + +# peaks (left side) +l = np.pi - array(160, 136, 98)*np.pi/180 #166, 155, 137 + +# calculate wavelength +l = sample(*np.append(np.sin(r)/n, np.sin(l)/n)*d) + +print(mformat(''' +# double slit interference +λ: {} m +''', l.tval())) + + +## +## lloyd mirror +## + +a = 37.5e-2 # transmitter position (from axis) +b = 32.5e-2 # receiver position (from axis) +m = 11.0e-2 # reflector position (first minimum) +h = array(14.4, 17.8, 20.8, 23.5)*1e-2 # reflector position (maxima) + +# calculate wavelength (first distance) +n = 1 + np.arange(1, len(h)+1) +l1 = sample(*(np.sqrt(a**2 + h**2) + np.sqrt(b**2 + h**2) - (a+b))/n).tval() + +a = 27.5e-2 # transmitter position (from axis) +b = 27.5e-2 # receiver position (from axis) +m = 9.5e-2 # reflector position (first minimum) +h = array(12.6, 15.8, 18.4, 20.1)*1e-2 # reflector position (maxima) + +# calculate wavelength (second distance) +n = 1 + np.arange(1, len(h)+1) +l2 = sample(*(np.sqrt(a**2 + h**2) + np.sqrt(b**2 + h**2) - (a+b))/n).tval() + + +print(mformat(''' +# lloyd mirror interference +λ₁: {} m +λ₂: {} m + +compatibility test: α={:.3f} +λbest: {} m +''', l1, l2, check_measures(l1, l2), combine_measures(l1, l2))) + + +## +## Fabry-Perot +## + +# reflection relative position (maxima) +m = array(19.4, 20.9, 22.3, 23.8, 25.3, 26.7, 28.0, 29.3, 30.7, 32.2)*1e-2 + +# calculate wavelength +l = sample(*np.diff(m)*2) + +print(mformat(''' +# fabry-perot interferometer +λ: {} m +''', l.tval())) + + +## +## Michelson +## + +# interferometer arms length (m) +l1 = 0.307 +l2 = 0.423 +l3 = 0.246 + +# reflector position (m) +d = array(0.092, 0.106, 0.120, 0.133, 0.150, 0.165, 0.180, 0.194, + 0.220, 0.235, 0.245, 0.260, 0.276, 0.290, 0.305, 0.322, + 0.338, 0.347, 0.360, 0.376, 0.392, 0.410) + +# fit and plot n - d +n = np.arange(len(d)) +a,b = linear(n, d, 2e-3) +f = lambda x: a.n + b.n*x +alpha = chi_squared_fit(n, d, f, 4e-3) + +print(mformat(''' +# michelson interferometer +λ: {} m +χ² test: α={:.3f} +''', 2*b, alpha)) + +plt.figure(3) +plt.title('michelson interferometer') +plt.xlabel('n-th maximum') +plt.ylabel('reflector position (m)') +plt.scatter(n, d, color='#178200') +plt.plot(n, f(n), color='#1ead00') +plt.show() diff --git a/optics/spectra.py b/optics/spectra.py new file mode 100644 index 0000000..1c68851 --- /dev/null +++ b/optics/spectra.py @@ -0,0 +1,192 @@ +# 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 degrees to radians +angle = lambda x, y, z: (x + y/60 + z/3600) * np.pi/180 + +t0 = angle(180,0,0) # angle offset +te = angle(0,0,1) # angle uncertainty + + + +## +## Part I +## + +l = 5893e-10 # sodium doublet wavelength +n = array(1, 2) # line order + +# line angles +t = array(angle(200,0,16) - angle(159,0,8), + angle(225,0,24) - angle(134,1,5))/2 + +# spacing (m⁻¹) +d = ufloat((n*l/np.sin(t)).mean(), 0.01e-6) + +print(mformat(''' +# Part I: diffraction grating spacing + +## Na lamp +d: {} lines/mm +''', 1/d*1e-3)) + + +## +## Part II +## + + +## lamp B + +# line angles +t = uarray(te, angle(195,0,2 ), # violet + angle(195,0,15), + angle(195,1,1 ), + angle(195,1,20), + angle(195,1,30), + angle(199,1,9 ), # green + angle(200,1,21)) # orange + +# wavelengths +l = d*unp.sin(t-t0) + +print(''' +# Part II + +## lamp B - Xe? +λ:''', '\n '.join(mformat('{}', i) for i in l*1e9), 'nm') + + +## lamp C + +# line angles +t = uarray(te, angle(193,1,11), # violet + angle(195,1,14), + angle(196,1,11), # green + angle(197,0,25), + angle(197,1,23), + angle(200,1,27), + angle(204,0,1 )) + +# wavelengths +l = d*unp.sin(t-t0) + +print(''' +## lamp C - Kr? +λ:''', '\n '.join(mformat('{}', i) for i in l*1e9), 'nm') + + + +## lamp D + +# line angles +t = uarray(te, angle(194,1,12), # violet + angle(194,1,23), + angle(195,0,8 ), + angle(195,0,15), + angle(195,0,25), + angle(196,0,2 ), # blue + angle(196,0,10), + angle(196,1,13), # green + angle(198,0,15), + angle(198,0,22), + angle(198,1,5 ), + angle(199,1,2 ), # yellow + angle(199,1,12), + angle(200,0,0 ), + angle(200,0,11), # orange + angle(200,0,23), + angle(201,0,2 ), + angle(201,0,8 ), + angle(201,1,11)) # red + +# wavelengths +l = d*unp.sin(t-t0) + +print(''' +## lamp D +λ:''', '\n '.join(mformat('{}', i) for i in l*1e9), 'nm') + + +## +## Part III +## + +## Hg lamp + +# prism angle +a = ufloat(angle(60,0,0), te) + +# spectral line angles +t = array( + (angle(230,50,11), angle(230,44, 1), angle(230,1,21), angle(229,1,19), angle(228,16,4), angle(228,1,25), angle(228,1,19)), + (angle(230,50,25), angle(230,44,10), angle(230,1,28), angle(229,1,21), angle(228,16,4), angle(228,1,22), angle(228,1,21)), + (angle(230,50,11), angle(230,44, 5), angle(230,1,15), angle(229,1,12), angle(228,16,3), angle(228,1,13), angle(228,1,12)), + (angle(230,50,15), angle(230,44, 7), angle(230,1,20), angle(229,1,14), angle(228,16,2), angle(228,1,11), angle(228,1,10))) + +# calculate mean values +t = uarray(1e-3, *map(np.mean, t.T)) + +# refractive index +n = unp.sin((a+t-t0)/2) / unp.sin(a/2) + +# spectral lines wavelengths +l = array(4047e-10, # violet I + 4077e-10, # violet II + 4358e-10, # blue + 4916e-10, # green I + 5460e-10, # green II + 5769e-10, # yellow I + 5790e-10) # yellow II + +# fit and plot n(λ) +f, (a,b) = curve(l, nominal(n), lambda x, a, b: a+b/x**2, sigma(n), guess=[1, 0.001]) +alpha = chi_squared_fit(l, nominal(n), f, sigma(n)) +x = np.linspace(400e-9, 580e-9, 100) + +print(mformat(''' +# Part III + +## Hg lamp +a: {} +b: {} + +χ² test: α={} +''', a, b, alpha)) + +plt.figure(1) +plt.xlabel('wavelength (nm)') +plt.ylabel('refractive index') +plt.scatter(l*1e9, nominal(n), color='#4f28fc') +plt.errorbar(l*1e9, nominal(n), yerr=sigma(n), linestyle='', color='#286bfc') +plt.plot(x*1e9, f(x), color='#5671c9') +plt.show() + + +## lamp B +t = array(angle(124,1,15), # violet + angle(125,0,0 ), + angle(125,0,5 ), + angle(125,0,16), # blue + angle(125,0,24), + angle(127,0,0 ), # green + angle(127,0,19), # yellow + angle(127,0,26), # orange + angle(127,1,12), # red + angle(128,0,25)) # red + +## lamp C +t = array(angle(232,1,24), # red + angle(233,0,18), + angle(234,1,1 ), # yellow + angle(236,1,3 ), # green + angle(237,0,7 ), + angle(237,1,28), # blue + angle(239,0,10)) # violet diff --git a/optics/spectra/ar b/optics/spectra/ar new file mode 100644 index 0000000..1246ae4 --- /dev/null +++ b/optics/spectra/ar @@ -0,0 +1,24 @@ + I λ + 100 4131.724 + 200 4277.528 + 250 4348.064 + 130 4426.001 + 130 4545.052 + 130 4579.350 + 130 4589.898 + 200 4609.567 + 130 4657.901 + 200 4726.868 + 100 4735.906 + 250 4764.865 + 200 4806.020 + 250 4879.864 + 300 6965.431 + 300 7067.218 + 300 7383.980 + 600 7503.869 + 400 7514.652 + 700 7635.106 + 400 7723.761 + 300 7724.207 + 600 7948.176 diff --git a/optics/spectra/he b/optics/spectra/he new file mode 100644 index 0000000..f48b398 --- /dev/null +++ b/optics/spectra/he @@ -0,0 +1,10 @@ + I λ + 200 3888.6456 + 300 3888.6489 + 200 4471.479 + 100 5015.678 + 500 5875.6148 + 250 5875.6404 + 120 5875.9663 + 200 6678.1517 + 100 7065.1771 diff --git a/optics/spectra/hg b/optics/spectra/hg new file mode 100644 index 0000000..f72071c --- /dev/null +++ b/optics/spectra/hg @@ -0,0 +1,10 @@ + I λ + 600 3650.153 + 1000 3983.931 + 400 4046.563 + 100 4347.494 + 1000 4358.328 + 500 5460.735 + 200 5677.105 + 250 6149.475 + 250 7944.555 diff --git a/optics/spectra/kr b/optics/spectra/kr new file mode 100644 index 0000000..e348e78 --- /dev/null +++ b/optics/spectra/kr @@ -0,0 +1,48 @@ + I λ + 100 3718.02 + 150 3778.046 + 150 3783.095 + 100 4057.037 + 100 4065.128 + 150 4088.337 + 150 4273.9694 + 200 4292.923 + 150 4317.81 + 150 4319.5794 + 1000 4355.477 + 130 4376.1216 + 100 4386.54 + 150 4431.685 + 200 4436.812 + 100 4453.9175 + 130 4463.6900 + 250 4475.014 + 130 4489.88 + 100 4502.3543 + 130 4523.14 + 250 4577.209 + 100 4582.978 + 150 4615.292 + 300 4619.166 + 250 4633.885 + 700 4658.876 + 150 4680.406 + 1000 4739.002 + 100 4762.435 + 300 4765.744 + 100 4811.76 + 100 4825.18 + 250 4832.077 + 250 4846.612 + 100 4945.59 + 130 5125.73 + 150 5208.32 + 150 5333.41 + 300 5570.2894 + 130 5681.89 + 500 5870.9160 + 100 6420.18 + 130 7289.78 + 130 7407.02 + 100 7524.46 + 150 7587.4136 diff --git a/optics/spectra/n b/optics/spectra/n new file mode 100644 index 0000000..c3304e7 --- /dev/null +++ b/optics/spectra/n @@ -0,0 +1,5 @@ + I λ + 150 5752.50 + 150 7423.64 + 200 7442.29 + 200 7468.31 diff --git a/optics/spectra/na b/optics/spectra/na new file mode 100644 index 0000000..8a79b52 --- /dev/null +++ b/optics/spectra/na @@ -0,0 +1,9 @@ + I λ + 500 3711.07 + 150 4113.70 + 600 4392.81 + 400 4405.12 + 500 4455.23 + 400 4490.87 + 1000 5889.950 + 500 5895.924 diff --git a/optics/spectra/ne b/optics/spectra/ne new file mode 100644 index 0000000..fe226f4 --- /dev/null +++ b/optics/spectra/ne @@ -0,0 +1,63 @@ + I λ + 250 3713.079 + 250 3727.107 + 800 3766.259 + 1000 3777.133 + 100 3818.427 + 120 3829.749 + 150 4219.745 + 100 4233.850 + 120 4250.649 + 120 4369.862 + 150 4379.550 + 100 4385.059 + 200 4391.991 + 150 4397.990 + 150 4409.299 + 100 4413.215 + 100 4421.389 + 100 4428.516 + 100 4428.634 + 150 4430.904 + 150 4430.942 + 120 4457.049 + 100 4522.720 + 100 4537.7545 + 100 4569.057 + 150 4704.3949 + 120 4708.8594 + 100 4710.0650 + 150 4712.0633 + 150 4715.344 + 100 4788.9258 + 100 4827.338 + 100 4884.9170 + 100 5341.0938 + 200 5400.5618 + 200 5852.4879 + 100 5881.8952 + 100 6029.9969 + 100 6074.3377 + 100 6143.0626 + 100 6163.5939 + 100 6217.2812 + 100 6266.4950 + 100 6334.4278 + 100 6382.9917 + 200 6402.248 + 150 6506.5281 + 100 6598.9529 + 1000 6929.4673 + 300 7024.0504 + 800 7032.4131 + 100 7059.1074 + 800 7173.9381 + 150 7213.200 + 150 7235.188 + 800 7245.1666 + 150 7343.945 + 300 7488.8712 + 100 7492.102 + 150 7522.818 + 300 7535.7741 + 130 7544.0443 diff --git a/optics/spectra/xe b/optics/spectra/xe new file mode 100644 index 0000000..3ba3d53 --- /dev/null +++ b/optics/spectra/xe @@ -0,0 +1,66 @@ + I λ + 300 4180.10 + 150 4193.15 + 100 4208.48 + 100 4213.72 + 100 4223.00 + 130 4238.25 + 150 4245.38 + 150 4296.40 + 150 4310.51 + 300 4330.52 + 150 4393.20 + 150 4395.77 + 150 4448.13 + 300 4462.19 + 150 4480.86 + 700 4844.33 + 130 4972.71 + 100 4988.77 + 300 5080.62 + 100 5122.42 + 100 5188.04 + 130 5191.37 + 150 5260.44 + 150 5261.95 + 700 5292.22 + 100 5309.27 + 300 5313.87 + 700 5339.33 + 150 5372.39 + 1000 5419.15 + 250 5438.96 + 100 5445.45 + 130 5460.39 + 300 5472.61 + 200 5531.07 + 100 5616.67 + 100 5659.38 + 200 5667.56 + 150 5726.91 + 150 5751.03 + 100 5758.65 + 100 5776.39 + 100 5893.29 + 150 5945.53 + 100 5971.13 + 700 5976.46 + 300 6036.20 + 700 6051.15 + 200 6093.50 + 500 6097.59 + 130 6101.43 + 150 6194.07 + 150 6270.82 + 130 6277.54 + 130 6343.96 + 200 6356.35 + 100 6512.83 + 300 6595.01 + 130 6597.25 + 100 6694.32 + 300 6805.74 + 250 6942.11 + 700 6990.88 + 150 7164.83 + 100 7548.45 diff --git a/shell.nix b/shell.nix new file mode 100644 index 0000000..0d1da3f --- /dev/null +++ b/shell.nix @@ -0,0 +1,22 @@ +{ pkgs ? import {}, 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 + ''; +}