Analysis of Preliminary Trail Data Pulled from Cavendish Balance



In [1]:

    
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
import math as m
from scipy.signal import argrelextrema as argex
plt.style.use('ggplot')



In [2]:

    
data_dir = '../data/'
trial_data = np.loadtxt(data_dir+'20171010_cavendish_trial.txt', delimiter='\t')



In [3]:

    
plt.plot(trial_data[:,0], trial_data[:,1])
plt.title("Trial Data from Cavendish Balance")
plt.ylabel("Anglular Positon (mrads)")
plt.xlabel("Time (s)")









    Out[3]:





<matplotlib.text.Text at 0x1042161d0>

The weird behavior at the beginning occured when we were making an alteration to the experimental setup itself (doing one of our many adjustments to try and zero out our $\theta_e$. We can ignore this as it is not indicative of our data and look at where it moves into a damped harmonic oscialltion.



In [30]:

    
plt.plot(trial_data[3900:,0], trial_data[3900:,1])
plt.title("Trial Data from Cavendish Balance, Adjusted")
plt.ylabel("Anglular Positon (mrads)")
plt.xlabel("Time (s)")









    Out[30]:





<matplotlib.text.Text at 0x111e6bb70>



In [5]:

    
np.savetxt(data_dir+'20171010_cavendish_trial_useable.txt', trial_data[3200:,], delimiter=',')



In [31]:

    
x_data = trial_data[3900:,0]
y_data = trial_data[3900:,1]



In [32]:

    
#I want to try to extract the peaks of the corresponding sine waves to show the exponential decay and use that 
#to fit my curve to.

angles = []
time = []

for i in range(1, len(y_data)-1):
    if y_data[i] >= np.average(y_data):
        if (y_data[i] > y_data[i-1]) and (y_data[i] > y_data[i+1]):
            angles.append(float(y_data[i]))
            time.append(float(x_data[i]))
            
inds = argex(np.array(angles), np.greater)
times = np.array(time)[inds]
turning_points = np.array(angles)[inds]
plt.plot(times, turning_points, 'g--')
plt.plot(x_data, y_data, 'b-')









    Out[32]:





[<matplotlib.lines.Line2D at 0x111ea2a20>]

Trying out the scipy.optimzie library to fit this to a decaying sinuisodal curve.



In [33]:

    
from scipy.optimize import curve_fit



In [41]:

    
def decay(t, a, b, w, phi, theta_0):
    return a*np.exp(-b*t)*np.cos(w*t + phi) + theta_0



In [48]:

    
popt, pcov = curve_fit(decay, x_data, y_data , p0 = (-32, 1.3e-3, 3e-2, -6e-1, 0))



In [49]:

    
popt









    Out[49]:





array([ -3.17805302e+01,   1.33728675e-03,   3.08477941e-02,
        -6.52543320e-01,  -3.01036426e+01])



In [69]:

    
plt.plot(times, turning_points, 'b--', label = 'Decay from Data')
plt.plot(x_data, y_data, 'r-', linewidth = 5, label = "Raw Data")
plt.plot(x_data, decay(x_data, *popt), 'g-', label = 'Fit of a*np.exp(-b*t)*np.cos(w*t + phi) + theta_0')
plt.title("Free Oscillation Data, No Large Masses")
plt.ylabel("Angle, Not Calibrated (mrad)")
plt.xlabel("Time (s)")
plt.legend()









    Out[69]:





<matplotlib.legend.Legend at 0x1164712e8>



In [68]:

    
round(popt[1],5)









    Out[68]:





0.00134

$$b = 1.34 \times 10 ^{-3} \frac{1}{s}$$



In [ ]: