gistfile1.txt

################################################################################
# undriven RLC circuit 
# solved for analytical solution in Maxima
# see these MIT OCW notes for a good intro:
#  http://ocw.mit.edu/courses/physics/8-02sc-physics-ii-electricity-and-magnetism-fall-2010/undriven-rlc-circuits/
# Author: Josh Stults
# Date: 28 Aug 2011
# Website: www.variousconsequences.com  
#
# governing equations as defined in Maxima:
# depends(I,t);
# gov_eqn : expand((diff(I,t,2) * L + R*diff(I,t) + I/C = 0)/L);
# sol_pos : ode2(gov_eqn, I, t); /* C(CR^2-4L) > 0 */
# sol_neg : ode2(gov_eqn, I, t); /* C(CR^2-4L) < 0 */
# sol_zero : ode2(gov_eqn, I, t); /* C(CR^2-4L) == 0 */
# sol_pos_ic : ic2(sol_pos, t=0, I=0, diff(I,t)=v0/L);
# sol_neg_ic : ic2(sol_neg, t=0, I=0, diff(I,t)=v0/L);
# sol_zero_ic : ic2(sol_zero, t=0, I=0, diff(I,t)=v0/L);

import scipy as sp
from scipy.integrate import ode 
from scipy.optimize import fmin 
bootindex = sp.random.random_integers 
import matplotlib.pylab as plt

# dump these expressions from Maxima using f90; then add 'sp.' and remove '&'  
def undriven_rlc(t,v0,R,L,C):
    check = C*(C*R**2 - 4*L)
    if(check>0):
        I = v0*sp.exp(t*(sp.sqrt(R**2/L**2-4/(C*L))-R/L)/2.0)/(L*sp.sqrt(R**2/L**2-4/(C*L)))-v0*sp.exp(t*(-sp.sqrt(R**2/L**2-4/(C*L))-R/L)/2.0)/(L*sp.sqrt(R**2/L**2-4/(C*L)))
    elif(check==0):
        I = t*v0*sp.exp(-t*R/L/2.0)/L
    else:
        I = 2*v0*sp.exp(-t*R/L/2.0)*sp.sin(t*sp.sqrt(4/(C*L)-R**2/L**2)/2.0)/(L*sp.sqrt(4/(C*L)-R**2/L**2))
    return(I)

def undriven_rlc_ls(X,v0,t1,i1,t2,i2):
    R = X[0]
    L = X[1]
    C = X[2] 
    err_1 = undriven_rlc(t1,v0,R,L,C) - i1
    err_2 = undriven_rlc(t2,v0,R,L,C) - i2 
    return(sp.sqrt(sp.sum(sp.append(err_1,err_2)**2)))

################################################################################

v0 = 800.0 # V initial potential on the capacitor 
R = 39e-3 + 1e-3 + 33e-3 + 227e-3 # Ohms
C = 3e-4 # F
L = (C*R**2)/4.0 #6.75e-6 # H
X0 = sp.array([R,L,C])
print (C*R**2)/4.0 
print C*(C*R**2 - 4*L)

t = sp.linspace(0,1e-3,4000)
I0 = undriven_rlc(t,v0,R,L,C)
d = 0.25 
Ip = undriven_rlc(t,v0,R,L/d,C)
Im = undriven_rlc(t,v0,R,L*d,C)

# read in g3data scraped current traces 
t1,i1,t1err,i1err = sp.loadtxt('IMG_1663.JPG.dat',unpack=True)
t2,i2,t2err,i2err = sp.loadtxt('IMG_1663.JPG_2.dat',unpack=True)
toff = t1.min()
t1 = t1 - toff
t2 = t2 - toff

# fit least squares estimate to the data
xopt = fmin(undriven_rlc_ls, X0,(v0,t1,i1,t2,i2)) 
tfit = sp.linspace(0,0.0025,4000)
ifit = undriven_rlc(tfit,v0,xopt[0],xopt[1],xopt[2]) 
# take a look at the autocorelation of the residuals
resid1 = i1 - undriven_rlc(t1,v0,xopt[0],xopt[1],xopt[2])
fresid1 = sp.fft(resid1)
autocor1 = sp.ifft(fresid1 * fresid1.conj()).real
autocor1 = autocor1/autocor1.max()
resid2 = i2 - undriven_rlc(t2,v0,xopt[0],xopt[1],xopt[2])
fresid2 = sp.fft(resid2)
autocor2 = sp.ifft(fresid2 * fresid2.conj()).real
autocor2 = autocor2/autocor2.max()
# bootstrap the residuals to give confidence intervals 
nboots = 100
xboot = sp.zeros((nboots,3))
for i in xrange(nboots):
    bi1 = bootindex(0,i1.shape[0]-1,i1.shape[0])
    bi2 = bootindex(0,i2.shape[0]-1,i2.shape[0])
    xboot[i] = fmin(undriven_rlc_ls, xopt, (v0,t1,i1+resid1[bi1],t2,i2+resid2[bi2]))

plt.figure() 
plt.plot(t, I0, label='L=%g H'%L)
plt.plot(t, Ip, label='L=%g H'%(L/d))
plt.plot(t, Im, label='L=%g H'%(L*d))
plt.xlabel('t (s)')
plt.ylabel('I (A)')
plt.legend(loc=0) 
ax = plt.gca()
plt.text(0.5,0.5,'initial voltage: %g V\nresistance: %g Ohms\ncapacitance: %g F'%(v0,R,C),fontsize='large',transform=ax.transAxes)
plt.savefig("undriven_rlc_pulse.png")

plt.figure()
plt.plot(t1,i1,'o') 
plt.plot(t2,i2,'o')
plt.plot(tfit,ifit,'k-')  
ax = plt.gca()
plt.text(0.2,0.6,'Given\n  initial voltage: %g V\n  IMG_1663.JPG scraped w/g3data\nFit\n  resistance: %1.2e [%1.2e,%1.2e] Ohms\n  inductance: %1.3g [%1.2e,%1.2e] H\n  capacitance: %1.2e [%1.2e,%1.2e] F'%(v0,xopt[0],xboot[:,0].min(),xboot[:,0].max(),xopt[1],xboot[:,1].min(),xboot[:,1].max(),xopt[2],xboot[:,2].min(),xboot[:,2].max()),fontsize='large',transform=ax.transAxes)
plt.xlabel('t (s)')
plt.ylabel('I (A)')
plt.legend(loc=0) 
plt.savefig("undriven_rlc_fit.png")

plt.figure()
# significant autocorrelation? 
plt.plot(autocor1[0:autocor1.shape[0]/2],'o')
plt.plot(autocor2[0:autocor2.shape[0]/2],'o')
plt.ylabel('autocorrelation')
plt.xlabel('lag')
plt.savefig("fit_resid_autocor.png")

plt.figure()
# look about normal?
nbins = 24
plt.hist(resid1,bins=nbins,normed=True,alpha=0.3)
plt.hist(resid2,bins=nbins,normed=True,alpha=0.3)
plt.savefig("fit_resid_hist.png")
plt.show()