In [1]:

    

import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
import scipy.optimize as opt


    

In [2]:

    

#2.6mm wide
#1.9mm verticle
#Taken at 640x480
mmperpix = .004
fps = 1000
n = 4 #Number of epx files
#epx files renamed test0.txt through testn.txt
dt = 1/fps


    

In [3]:

    

def filterindex(maj):
    z = 0
    for entry in maj:
        if entry<2*maj[0] and entry>.5*maj[0] and entry<400:
            """When the particle tracker loses the particle it creates a very large ellipse where the minor axis is
            about the width of the screen. Also if tracking a small particle that overlaps with a dust spec, sometimes
            the tracker will shrink the size of the ellipse and fit it to the dust. This if statement along with the 
            list slicing below discard any entries in the position arrays after there has been an error in tracking"""
            z = z + 1
    return z


    

In [4]:

    

alldata = np.array([np.delete(np.loadtxt('new' + str(j) +'.txt', unpack=True),0,0) for j in range(n)])
alldata = alldata*mmperpix
alldata = np.array([entry for box in alldata for entry in box])


    

In [5]:

    

N = int(len(alldata)/5)
allx = np.array([alldata[5*i] for i in range(N)])
ally = np.array([alldata[1+5*i] for i in range(N)])
allmaj = np.array([alldata[2+5*i] for i in range(N)])
allmin = np.array([alldata[3+5*i] for i in range(N)])


    

In [6]:

    

x = np.array([allx[i][0:filterindex(allmaj[i])] for i in range(len(allx))])
x = np.array([entry for entry in x if len(entry)>=20])
y = np.array([ally[i][0:filterindex(allmaj[i])] for i in range(len(ally))])
y = np.array([entry for entry in y if len(entry)>=20])


    

In [7]:

    

xavg = np.array([np.mean(entry) for entry in x])
alldy = np.array([np.gradient(entry) for entry in y])
dyavg = -np.array([np.mean(entry) for entry in alldy])
vyavg = dyavg/dt


    

In [8]:

    

particles = len(x)
particles


    

    
Out[8]:





32

In [9]:

    

def u(x,dp):
    b = 1.25*.001
    mu = .045
    return -(x/mu)*dp*(b-x/2)


    

In [10]:

    

theta_best, theta_cov = opt.curve_fit(u, xavg[vyavg>0]*.001, vyavg[vyavg>0]*.001)
dp = theta_best[0]
dp


    

    
Out[10]:





-264.1062560520308

In [11]:

    

print('dp/dx = {0:.3f} +/- {1:.3f}'.format(dp, np.sqrt(theta_cov[0,0])))


    

    




dp/dx = -264.106 +/- 3.803

In [16]:

    

plt.title('Average Y Velocity vs Average X Position');
plt.xlabel('X Position (mm)');
plt.ylabel('Y Velocity (mm/s)');
plt.scatter(xavg[vyavg>.2], vyavg[vyavg>.2]);
plt.plot(np.linspace(0,2.5,1000),u(np.linspace(0,2.5*.001,1000),dp)*1000)
plt.tick_params(axis='x',top='off',direction='out')
plt.tick_params(axis='y',right='off',direction='out')
plt.text(.7, 1, 'dp/dx = {0:.3f} +/- {1:.3f} Pa/m'.format(dp, np.sqrt(theta_cov[0,0])));
plt.text(.7, .5, 'Number of Particles = {} '.format(particles));
plt.plot(np.zeros(5),np.linspace(0,5,5),color = 'black',linewidth = '6')
plt.plot(np.array([2.5]*5),np.linspace(0,5,5),color = 'black',linewidth = '6')
plt.xlim(0,2.5)
plt.ylim(0,5)
plt.savefig('Firstrun.png')


    

In [ ]: