Substituting the more accurate expression for $\eta$ into the expression $-m'g - 6\pi\eta a v_g$ gives:
$$
-\frac{4}{3}\pi a^3\rho'g = \frac{6\pi\eta_0av_g}{1+\frac{b}{ap}},
$$
which rearranges to a quadratic equation in the droplet radius, $a$:
$$
a^2 + \frac{b}{p}a + \frac{9\eta_0v_g}{2\rho'g} = 0.
$$
Taking the positive root gives the following expression for $a$ (recall that $v_g < 0$):
$$
a = -\frac{b}{2p} + \frac{1}{2}\sqrt{\frac{b^2}{p^2} - \frac{18\eta_0v_g}{\rho'g}},
$$
which is used in the function get_a below.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
df = pd.read_csv('eg10-millikan-data.txt', delim_whitespace=True, index_col=[0, 1])
# Air viscosity parameters, eta0 (kg.m-1.s-1) and b (Pa.m).
eta0 = 1.869e-5
b = 7.88e-3
# Air density, oil density, density difference (kg.m-3).
rho_air = 1.17
rho_oil = 917.3
rhop = rho_oil - rho_air
# Air pressure (Pa), acceleration due to gravity (m.s-2).
p = 100.82e3
g = 9.803
# Fall/rise distance (m), electric field (V.m-1).
d = 11.09e-3
E = -322.1e3
def get_eta(a):
"""Return the viscosity adjusted for droplet radius."""
return eta0 / (1 + b / a / p)
def get_a(tg):
"""Get the droplet radius from fall time."""
# Droplet fall terminal velocity.
vg = - d / tg
c = b / p
a = (-c + np.sqrt(c**2 - 18 * eta0 * vg / rhop / g)) / 2
return a
for drop in df.index.levels[0]:
drop_df = df.loc[drop].T
# Rename the index to simple 'tg' and 'te'.
drop_df.index = drop_df.index.str.slice(0,2)
# Average over all tg values.
tg = drop_df.loc['tg'].values.mean()
# Average te for each experiment.
te = drop_df.loc['te'].mean()
a = get_a(tg)
eta = get_eta(a)
q = 6 * np.pi * eta * a * d / E * (1/tg + 1/te)
df.loc[drop, 'q'] = q.values
sorted_q = sorted(-df.loc[:, 'q'])
dq = np.diff(sorted_q)
e_estimate = dq[(dq>1.e-19) & (dq<2.e-19)].mean()
df['N'] = (df['q'] / e_estimate).astype(int)
e = (df['q']/df['N']).mean()
print('Estimated value for e = {} C'.format(e))
The output is a slightly better estimate of $e$:
Estimated value for e = 1.5976041889076869e-19 C