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