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