import numpy as np
import pylab as plt
_volts, _A, _B, _C, _AB, _AC, _BC, _ABC = np.loadtxt('spi_indistinguishable.txt',\
unpack=True, skiprows=9, usecols=(2,3,4,5,6,7,8,9), \
delimiter=',')
print '_volts =', _volts
print '_ABC =', _ABC
print '_AC =', _AC
volts = np.arange(40,71)
print 'volts =', volts
N = len(volts)
AC = np.zeros(N)
AC_err = np.zeros(N)
nsamp = np.zeros(N)
print 'range(N) =', range(N)
for i in range(N):
print "i =", i
print "voltage =", volts[i]
indices = np.where(_volts == volts[i])
print "indices =", indices
print "_volts at indices =", _volts[indices]
print "_AC at indices =", _AC[indices]
AC[i] = np.mean(_AC[indices])
AC_err[i] = np.std(_AC[indices])
nsamp[i] = len(_AC[indices])
plt.clf()
plt.errorbar(volts, AC, AC_err/np.sqrt(nsamp), fmt='r.', label='AC')
plt.xlabel('Piezoelectric Voltage (V)')
plt.ylabel('Count (per 0.01 s)')
plt.title('Indistiguishable Paths')
plt.legend(loc='lower left')
plt.savefig('spi_data.png')
# Try linear fit a*cos(2*pi*V/period)+b*sin(2*pi*V/period)+c
period=9 #Volts
cos = np.cos(2*np.pi*volts/period)
sin = np.sin(2*np.pi*volts/period)
one = np.ones(len(volts))
MT = np.array([cos, sin, one]) #Transpose of projection matrix
M = np.transpose(MT) #Projection matrix
Ninv = np.diag(1/(AC_err**2/nsamp)) #Inverse noise matrix
# np.dot function produces
D = np.dot(MT,np.dot(Ninv,M))
b = np.dot(MT,np.dot(Ninv,AC))
cov = np.linalg.inv(D) #Parameter covariance matrix
param = np.dot(cov,b) #Best fit parameters: a=(D^{-1}M^TN^-1)(data)
print "a =", param[0], '+/-', np.sqrt(cov[0][0])
print "b =", param[1], '+/-', np.sqrt(cov[1][1])
print "c =", param[2], '+/-', np.sqrt(cov[2][2])
plt.plot(volts, param[0]*cos+param[1]*sin+param[2]*one, label='Model')
plt.legend(loc='lower left')
plt.savefig("spi_data_linear_model.png")