df = pd.read_excel('Data.xlsx', sheetname=None)

df['Run 1']

keys = ['Run 1', 'Run 2', 'Run 3', 'Run 4']

# Fall time in seconds when V = 0
tf = np.array([df[key]['t_f (s)'] for key in keys])
tf

# Average Fall times for each run (seconds)
avg_tf = np.array([np.mean(tf[i]) for i in np.arange(4)]) 
avg_tf

# Standard Deviation Fall times for each run (seconds)
std_tf = np.array([np.std(tf[i]) for i in np.arange(4)])
std_tf

"""
Explanation: Read-In
End of explanation
"""


tf = np.array([ma.masked_outside(entry, np.mean(entry) - 2.5*np.std(entry), np.mean(entry) + 2.5*np.std(entry)) 
                     for entry in tf])

"""
Explanation: Mask points that lie outside 2.5$\sigma$
End of explanation
"""


# Standard deviation of the mean (seconds)
tf_u = np.array([std_tf[i]/np.sqrt(len(tf[i])) for i in np.arange(4)])

tf_u

"""
Explanation: Uncertainty
End of explanation
"""


# Units of m/s
speeds = np.array([(.5*1e-3)/entry for entry in avg_tf ])
speeds

"""
Explanation: Fall Speed
All measurements were of the time to travel .5mm
End of explanation
"""


#Units of m/s
speeds_u = np.array([((.5*1e-3)/avg_tf[i]**2)*tf_u[i] for i in np.arange(4)])
speeds_u

"""
Explanation: \begin{equation}
    v = \frac{h}{t} 
\end{equation}
\begin{equation}
    \delta v = \frac{h}{t^2} \delta t 
\end{equation}
End of explanation
"""


d = 7.63*1e-3 #Electrode spacer thickness in meters
V = 301.8 #Volts. Constant through all measurements

E = V/d #V/m
E

"""
Explanation: E Field
End of explanation
"""


Resistances = np.array([2.065, 2.011, 1.957, 1.959])*1e6 #From thermistor. Used to find temperature. Ohms

temperatures = np.array([24, 25, 26, 26]) #From Thermistor table. Celcius

rho = 885 #kg/m^3. Density of oil drop.
g = 9.80 #gravity

def viscosity(T):
    #Returns the viscosity of air for a given temp in celsius
    # Viscosity units [Ns/m^2]
    return (1.8 + (T-15)/209)*1e-5

viscosities = np.array([viscosity(entry) for entry in temperatures])
viscosities

"""
Explanation: Air Viscosity
End of explanation
"""


def radius(speed,visc):
    return np.sqrt(9*visc*speed/(2*g*rho))

radii = np.array([radius(speeds[i],viscosities[i]) for i in np.arange(4)]) #Blob radius in meters
radii

"""
Explanation: Drop Radius
End of explanation
"""


def mass(r):
    return (4/3)*np.pi*rho*r**3

masses = np.array([mass(entry) for entry in radii])
masses

"""
Explanation: Droplet Mass
End of explanation
"""


masses_u = np.array([np.pi*(radii[i]**2)*np.sqrt(18*viscosities[i]*rho/g)*np.sqrt(1/speeds[i])*speeds_u[i] for i in np.arange(4)])
masses_u

"""
Explanation: \begin{equation}
    m = \frac{4}{3} \pi r^3 \rho 
\end{equation}
\begin{equation}
    \delta m = 4 \pi r^2 \rho \delta r
\end{equation}
\begin{equation}
    r = \sqrt{\frac{9 \eta v_f}{2 g \rho}}
\end{equation}
\begin{equation}
    \delta r = \frac{1}{2}\sqrt{\frac{9 \eta}{2 g \rho}} \frac{\delta v_f}{\sqrt{v_f}}
\end{equation}
\begin{equation}
    \delta m =  \pi r^2 \sqrt{\frac{18 \eta \rho}{ g}} \frac{\delta v_f}{\sqrt{v_f}}
\end{equation}
End of explanation
"""


def gamma(r):
    b = .082*1e-7
    return (1+b/r)**(-3/2)

gammas = np.array([gamma(entry) for entry in radii])
gammas

"""
Explanation: Correction Factor
End of explanation
"""


# Rise time in seconds when V = +
tu = np.array([df[key]['t_u (s)'] for key in keys])
tu

#### tu0: Time Up 0'th drop.
#### Average Time
#### Std of Times
#### Resulting Speed
#### Speed Uncertainty

tu0 = ma.masked_array(tu[0],mask=[0,0,0,0,0,0,0,0,0,1])
avg_tu0 = np.mean(tu0)
std_tu0 = np.std(tu0)
tu0_speed = .5*1e-3/avg_tu0
tu0_speed_u = .5*1e-3/(avg_tu0**2)*std_tu0/np.sqrt(len(tu0))

tu1 = ma.masked_array(tu[1],mask=[0,0,1,1,0,0,0,0,0,0])
avg_tu1 = np.mean(tu1)
std_tu1 = np.std(tu1)
tu1_speed = .5*1e-3/avg_tu1
tu1_speed_u = .5*1e-3/(avg_tu1**2)*std_tu1/np.sqrt(len(tu1))

tu2_1 = ma.masked_array(tu[2],mask=[0,1,1,0,1,0,0,0,0,0])
avg_tu2_1 = np.mean(tu2_1)
std_tu2_1 = np.std(tu2_1)
tu2_1_speed = .5*1e-3/avg_tu2_1
tu2_1_speed_u = .5*1e-3/(avg_tu2_1**2)*std_tu2_1/np.sqrt(len(tu2_1))

tu2_2 = ma.masked_array(tu[2],mask=[1,0,0,1,0,1,1,1,1,1])
avg_tu2_2 = np.mean(tu2_2)
std_tu2_2 = np.std(tu2_2)
tu2_2_speed = .5*1e-3/avg_tu2_2
tu2_2_speed_u = .5*1e-3/(avg_tu2_2**2)*std_tu2_2/np.sqrt(len(tu2_2))

tu3 = tu[3]
avg_tu3 = np.mean(tu3)
std_tu3 = np.std(tu3)
tu3_speed = .5*1e-3/avg_tu3
tu3_speed_u = .5*1e-3/(avg_tu3**2)*std_tu3/np.sqrt(len(tu3))



#This data can't be used unless the free fall time data is analyzed
tu4 = np.array([4.20, 3.51, 3.44, 3.21, 3.56, 3.27, 3.25]) #This data from an extra column not entered into excel
avg_tu4 = np.mean(tu4)
std_tu4 = np.std(tu4)
tu4_speed = .5*1e-3/avg_tu4
tu4_speed_u = .5*1e-3/(avg_tu4**2)*std_tu4/np.sqrt(len(tu4))

avg_tu0

std_tu0/np.sqrt(len(tu0))

def up_charge(vu,vf,gamma,mass):
    return gamma*mass*g*((vu/vf) + 1)/E  

qu0 = up_charge(tu0_speed, speeds[0], gammas[0], masses[0])
qu1 = up_charge(tu1_speed, speeds[1], gammas[1], masses[1])
qu2_1 = up_charge(tu2_1_speed, speeds[2], gammas[2], masses[2])
qu2_2 = up_charge(tu2_2_speed, speeds[2], gammas[2], masses[2])
qu3 = up_charge(tu3_speed, speeds[3], gammas[3], masses[3])
qu = np.array([qu0,qu1,qu2_1,qu2_2,qu3])
qu*1e19

"""
Explanation: Upward Movement
End of explanation
"""


def delta_qu(qu,dm, m, dvu, vu, dvf, vf):
    return np.sqrt((dm/m)**2 + (dvu/(vu+vf))**2 + ((vu/vf)*dvf/(vf+vu))**2)*qu

qu0_u = delta_qu(qu0,masses_u[0],masses[0],tu0_speed_u, tu0_speed, speeds_u[0], speeds[0])
qu1_u = delta_qu(qu1,masses_u[1],masses[1],tu1_speed_u, tu1_speed, speeds_u[1], speeds[1])
qu2_1_u = delta_qu(qu2_1,masses_u[2],masses[2],tu2_1_speed_u, tu2_1_speed, speeds_u[2], speeds[2])
qu2_2_u = delta_qu(qu2_2,masses_u[2],masses[2],tu2_2_speed_u, tu2_2_speed, speeds_u[2], speeds[2])
qu3_u = delta_qu(qu3,masses_u[3],masses[3],tu3_speed_u, tu3_speed, speeds_u[3], speeds[3])
qu_u = np.array([qu0_u,qu1_u,qu2_1_u,qu2_2_u,qu3_u])
qu_u*1e19

"""
Explanation: \begin{equation}
    \frac{\delta q_u}{q_u} =  \sqrt{(\frac{\delta m}{m})^2 + (\frac{\delta v_u}{v_u + v_f})^2 + (\frac{v_u}{v_f}\frac{\delta v_f}{v_f+v_u})^2 }
\end{equation}
End of explanation
"""


df['Run 1']

# Fall time in seconds when V = -
td = np.array([df[key]['t_d (s)'] for key in keys])
td

#### tu0: Time Up 0'th drop.
#### Average Time
#### Std of Times
#### Resulting Speed
#### Speed Uncertainty



td0 = td[0] #All values here are about the same
avg_td0 = np.mean(td0)
std_td0 = np.std(td0)
td0_speed = .5*1e-3/avg_td0
td0_speed_u = .5*1e-3/(avg_td0**2)*std_td0/np.sqrt(len(td0))

td1 = np.array([td[1]]) #All values here are about the same
avg_td1 = np.mean(td1)
std_td1 = np.std(td1)
td1_speed = .5*1e-3/avg_td1
td1_speed_u = .5*1e-3/(avg_td1**2)*std_td1/np.sqrt(len(td1))

td2 = np.array([td[2]]) #All values here are about the same
avg_td2 = np.mean(td2)
std_td2 = np.std(td2)
td2_speed = .5*1e-3/avg_td2
td2_speed_u = .5*1e-3/(avg_td2**2)*std_td2/np.sqrt(len(td2))


td3 = td[3] #All values here are about the same
avg_td3 = np.mean(td3)
std_td3 = np.std(td3)
td3_speed = .5*1e-3/avg_td3
td3_speed_u = .5*1e-3/(avg_td3**2)*std_td3/np.sqrt(len(td3))

avg_td0

std_td0/np.sqrt(len(td0))

"""
Explanation: Drop Charge
Downward Movement
End of explanation
"""


def down_charge(vd,vf,gamma,mass):
    return gamma*mass*g*((vd/vf) - 1)/E  

def delta_qd(qu,dm, m, dvu, vu, dvf, vf):
    return np.sqrt((dm/m)**2 + (dvu/(vu-vf))**2 + ((vu/vf)*dvf/(vf-vu))**2)*qu

qd0 = down_charge(td0_speed, speeds[0], gammas[0], masses[0])
qd1 = down_charge(td1_speed, speeds[1], gammas[1], masses[1])
qd2 = down_charge(td2_speed, speeds[2], gammas[2], masses[2])
qd3 = down_charge(td3_speed, speeds[3], gammas[3], masses[3])
qd = np.array([qd0,qd1,qd2,qd3])
qd*1e19

qd0_u = delta_qd(qd0,masses_u[0],masses[0],td0_speed_u, td0_speed, speeds_u[0], speeds[0])
qd1_u = delta_qd(qd1,masses_u[1],masses[1],td1_speed_u, td1_speed, speeds_u[1], speeds[1])
qd2_u = delta_qd(qd2,masses_u[2],masses[2],td2_speed_u, td2_speed, speeds_u[2], speeds[2])
qd3_u = delta_qd(qd3,masses_u[3],masses[3],td3_speed_u, td3_speed, speeds_u[3], speeds[3])
qd_u = np.array([qd0_u,qd1_u,qd2_u,qd3_u])
qd_u*1e19

"""
Explanation: Drop Charge
End of explanation
"""


qu*1e19

#Check with known e
np.array([qu[0]/1.6, qu[1]/1.6, qu[2]/1.6, qu[3]/1.6,qu[4]/1.6])*1e19

up_e = np.array([qu[0]/5, qu[1]/7, qu[2]/4, qu[3]/3,qu[4]/4])
up_e*1e19

"""
Explanation: Search for Charge Unit
Upwards Charge
End of explanation
"""


qu_u*1e19

up_e_u = np.array([qu_u[0]/5, qu_u[1]/7, qu_u[2]/4, qu_u[3]/3,qu_u[4]/4])
up_e_u*1e19

"""
Explanation: Upwards Uncer
End of explanation
"""


qd*1e19

#Check with known e
np.array([qd[0]/1.6, qd[1]/1.6, qd[2]/1.6, qd[3]/1.6])*1e19

down_e = np.array([qd[0]/5, qd[1]/7, qd[2]/5, qd[3]/5])
down_e*1e19

"""
Explanation: Down Charge
End of explanation
"""


qd_u*1e19

down_e_u = np.array([qd_u[0]/5, qd_u[1]/7, qd_u[2]/5, qd_u[3]/5])
down_e_u*1e19

"""
Explanation: Down Uncer
End of explanation
"""


Charges = np.hstack((up_e,down_e))
Charges*1e19

Charge_Best = np.mean(Charges)
Charge_Best*1e19

Charges_u = np.hstack((up_e_u,down_e_u))
Charges_u*1e19

stdom_e = np.std(Charges)/np.sqrt(len(Charges))

#delta_e_meas = np.sqrt(np.sum(np.array([(Charges_u[i]/Charges[i])**2 for i in np.arange(len(Charges))])))*Charge_Best

delta_e = stdom_e

delta_e*1e19

print('e = (%.2f +/- %.2f)x10^(-19) C' % (Charge_Best*1e19,delta_e*1e19))

"""
Explanation: Determine Electron Charge
End of explanation
"""


# From a previous lab
eperm = 1.648*1e11 #C/kg
eperm_u = .287*1e11

m = Charge_Best/eperm
m

dm = np.sqrt( (delta_e/Charge_Best)**2 + (eperm_u/eperm)**2 )*m
dm

print('m = (%.3f +/- %.3f)x10^(-31) Kg' % (m*1e31,dm*1e31))

"""
Explanation: Charge to Mass Ratio
e/m = (1.648 +/- 0.287)x10^11 C/kg
End of explanation
"""


eo = 8.854*1e-12

E_r = 13.77 #eV
E_r = E_r*1.602*1e-19
dE_r = .01 #Ev

h = np.sqrt(m * Charge_Best**4/(8*eo**2 * E_r))
h

dh = np.sqrt((dm/m)**2 + (.01/13.77)**2 + (delta_e/Charge_Best)**2)*h
dh

print('h = (%.3f +/- %.3f)x10^(-34) Js' % (h*1e34,dh*1e34))

"""
Explanation: Plank Constant
End of explanation
"""


