In [1]:
import numpy as np
from scipy import integrate

import matplotlib as mpl       # As of July 2017 Bucknell computers use v. 2.x 
import matplotlib.pyplot as plt
mpl.style.use('classic')
plt.rc('figure', figsize = (6, 4.5))            # Reduces overall size of figures
plt.rc('axes', labelsize=16, titlesize=14)
plt.rc('figure', autolayout = True)
In [2]:
%matplotlib notebook
In [3]:
def yLin(t):
    tau = m/bA
    return -(vtA*t - vtA*tau*(1 - np.exp(-t/tau))) + y0

def vLin(t):
    tau = m/bA
    return vtA*(1 - np.exp(-t/tau))

def yQuad(t):
    return -vtB**2*np.log(np.cosh(g*t/vtB))/g + y0

def vQuad(t):
    return vtB*np.tanh(g*t/vtB)
In [4]:
data = np.loadtxt('sim.dat')
t, y = data.T
In [5]:
m = 0.0038
bA = 0.0175
bB = 0.01287
y0 = 1.153
vy0 = 0
g = 9.81
vtA = vt = m*g/bA 
vtB = np.sqrt(m*g/bB)
print(vtA, vtB)

2.1301714285714284 1.701911718187373
In [6]:
plt.figure()
plt.scatter(t,y)
t2 = np.linspace(0,0.6,601)
plt.plot(t, yLin(t))
plt.grid()
plt.xlabel('$t$')
plt.ylabel('$y$')
plt.title('linear drag; $B_A = %s$'%bA)
plt.ylim(0,1.2);
In [7]:
plt.figure()
plt.scatter(t, y - yLin(t))
t2 = np.linspace(0,0.6,601)
plt.grid()
plt.title('residuals with linear drag; $B_A = %s$'%bA);
#plt.ylim(0,1.2)

Sum of the squares of deviations:

In [8]:
sum_sq = np.sum((y-yLin(t))**2)
print('sqrt of sum of squared deviations = ', np.sqrt(sum_sq))

sqrt of sum of squared deviations =  0.06523985444355067
In [9]:
plt.figure()
plt.scatter(t,y)
t2 = np.linspace(0,0.6,601)
plt.plot(t, yQuad(t))
plt.grid()
plt.xlabel('$t$')
plt.ylabel('$y$')
plt.title('quadratic drag; $B_B = %s$'%bB)
plt.ylim(0,1.2);
In [10]:
plt.figure()
plt.scatter(t, y - yQuad(t))
t2 = np.linspace(0,0.6,601)
plt.grid()
#plt.ylim(0,1.2)
plt.title('residuals with quadratic drag; $B_B = %s$'%bB);
In [11]:
sum_sq = np.sum((y-yQuad(t))**2)
print('square root of sum of squared deviations = ', np.sqrt(sum_sq))

square root of sum of squared deviations =  0.0024367953182602574

Estimating $v$ from $y$ data

Use numpy features to make life easy

  • calculate $\Delta y$
  • calculate $\Delta t$
  • calculate $v$

NOTE: I plotted values of $v$ at the midpoint of the time intervals.

This improves things quite a bit.

In [12]:
y
Out[12]:
array([1.153, 1.148, 1.132, 1.107, 1.075, 1.036, 0.992, 0.944, 0.894,
       0.842, 0.789, 0.735, 0.68 , 0.625, 0.569, 0.514, 0.458, 0.402,
       0.345])
In [13]:
np.roll(y,-1)
Out[13]:
array([1.148, 1.132, 1.107, 1.075, 1.036, 0.992, 0.944, 0.894, 0.842,
       0.789, 0.735, 0.68 , 0.625, 0.569, 0.514, 0.458, 0.402, 0.345,
       1.153])
In [14]:
deltaY = y - np.roll(y,-1)
In [15]:
t
Out[15]:
array([0.   , 0.033, 0.066, 0.099, 0.132, 0.165, 0.198, 0.231, 0.264,
       0.297, 0.33 , 0.363, 0.396, 0.429, 0.462, 0.495, 0.528, 0.561,
       0.594])
In [16]:
np.roll(t,-1)
Out[16]:
array([0.033, 0.066, 0.099, 0.132, 0.165, 0.198, 0.231, 0.264, 0.297,
       0.33 , 0.363, 0.396, 0.429, 0.462, 0.495, 0.528, 0.561, 0.594,
       0.   ])
In [17]:
deltaT =  np.roll(t,-1) - t
In [18]:
v = deltaY[:len(deltaY)-1]/deltaT[:len(deltaY)-1]
In [19]:
plt.figure()
plt.scatter(t[:len(v)]+0.033/2,-v)
t2 = np.linspace(0,0.6,601)
plt.plot(t, -vLin(t))
plt.grid()
plt.xlabel('$t$')
plt.ylabel('$v$')
plt.title('linear drag; $B_A = %s$'%bA);
#plt.ylim(0,1.2)
In [20]:
plt.figure()
plt.xlabel('$t$')
plt.ylabel('$v$')
plt.title('quadratic drag; $B_B = %s$'%bB);
plt.scatter(t[:len(v)]+0.033/2,-v)
t2 = np.linspace(0,0.6,601)
plt.plot(t, -vQuad(t))
plt.grid()
#plt.ylim(0,1.2)
In [ ]: