SPI Data Reading and Linear Fit

This is a program that reads in the file media:spi_indistinguishable.py.

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")