## Determining Photogate Timing Uncertainty

In order to do an uncertainty analysis of Front Range’s conservation of momentum experiments, we need to determine the uncertainty of the times reported by the photogates. The key to experiments of this sort is to collect redundant measurements, use least squares techniques to come up with a model of the expected results and then take statistics on the differences between the measurements and the fitted model.

The results, developed below, show that the photogate timing uncertainty is roughly 1.3 ms, and that the resistance force is proportional to the glider’s velocity.

The experiment to determine the photogate timing uncertainties uses four photogates spaced evenly along an air track, connected to a laptop computer via a USB data acquisition system. A glider slides along the air track, and a fin on the glider interrupts a beam of infrared light across the photogate. The data acquisition system records the times that the beam transitions between the broken and clear states. The air track provides the glider with a low-friction environment, but the resistance is not zero, so we will model the glider as having a constant deceleration during each pass along the track.

In the experiment, the glider is given an initial velocity, and allowed to bounce off the ends of the track, reversing direction each time. The data acquisition system collects data for 60 seconds, allowing on the order of 10 round trips along the track for each run.

load run5
n = length(t_meas)
n_trips = fix(n/16)
n = n_trips*16
t_meas = t_meas(1:n);
l_fin = 0.07; % 7 cm fin length
t_est = (t_meas(1:2:end)+t_meas(2:2:end))/2;
v_est = l_fin./(t_meas(2:2:end) - t_meas(1:2:end));
plot(t_est, v_est, '+')
xlabel('Time (s)')
ylabel('Velocity (m/s)')
title('Run 5 velocity history from time measurements')
n =
182
n_trips =
11
n =
176

## Estimating Photogate Spacing

The glider goes through photogates 1, 2, 3 and 4, and then bounces off the far end of the track and goes through photogates 4, 3, 2 and 1. It then bounces off the near end of the track and repeats the process. I should be able to estimate the spacing of the photogates from the measurements, but I need to be sure to use the correct times to work out the three distances.

% Spacing between photogates 1 and 2, travelling right
d1_right = (v_est(2:8:end)+v_est(1:8:end))/2.*(t_est(2:8:end)-t_est(1:8:end));
% Spacing between photogates 1 and 2, travelling left
d1_left = (v_est(7:8:end)+v_est(8:8:end))/2.*(t_est(8:8:end)-t_est(7:8:end));
d1_data = [d1_right; d1_left];
d1 = mean(d1_data)
d1_unc = std(d1_data)
% Spacing between photogates 2 and 3, traveling right
d2_right = (v_est(3:8:end)+v_est(2:8:end))/2.*(t_est(3:8:end)-t_est(2:8:end));
% Spacing between photogates 2 and 3, traveling left
d2_left = (v_est(6:8:end)+v_est(7:8:end))/2.*(t_est(7:8:end)-t_est(6:8:end));
d2_data = [d2_right; d2_left];
d2 = mean(d2_data)
d2_unc = std(d2_data)
% Spacing between photogates 3 and 4, traveling right
d3_right = (v_est(4:8:end)+v_est(3:8:end))/2.*(t_est(4:8:end)-t_est(3:8:end));
% Spacing between photogates 3 and 4, traveling left
d3_left = (v_est(5:8:end)+v_est(6:8:end))/2.*(t_est(6:8:end)-t_est(5:8:end));
d3_data = [d3_right; d3_left];
d3 = mean(d3_data)
d3_unc = std(d3_data)
d1 =
0.2398
d1_unc =
7.4334e-005
d2 =
0.2380
d2_unc =
1.3359e-004
d3 =
0.2422
d3_unc =
1.9750e-004

## Analyzing the First Rightward Pass

This will take a few steps…

## Define the Locations of Block and Clear Events

t_pass = t_meas(1:8) - mean(t_meas(1:8))
x_photogate = d2*[-1 -1 +1 +1]/2 + [-d1 0 0 d3]
A = [1 0 0 0; 1 0 0 0; 0 1 0 0; 0 1 0 0; 0 0 1 0; 0 0 1 0; 0 0 0 1; 0 0 0 1]';
x_event = x_photogate*A + [-1 +1 -1 +1 -1 +1 -1 +1]*l_fin/2
t_pass =
-0.5413
-0.4459
-0.2138
-0.1179
0.1134
0.2101
0.4489
0.5465
x_photogate =
-0.3588 -0.1190 0.1190 0.3612
x_event =
-0.3938 -0.3238 -0.1540 -0.0840 0.0840 0.1540 0.3262 0.3962

## Basic Linear Fit

% Find the midpoints of the occluded intervals
t_mid = (t_pass(1:2:end) + t_pass(2:2:end))/2
% Now find estimated velocities
v_est = l_fin./(t_pass(2:2:end) - t_pass(1:2:end))
% Naive fit
[v_fit, v_unc, a, v_0, a_unc, v_0_unc] = uncertainFit(t_mid, v_est);
plot(t_mid, v_est, '+-', t_mid, v_fit, 'x:')
xlabel('Time (s)')
ylabel('Velocity (m/s)')
title('Naive Velocity Fit')
t_mid =
-0.4936
-0.1659
0.1618
0.4977
v_est =
0.7339
0.7293
0.7242
0.7169

## Find an initial condition to minimize the squared error

I’ll use the fitted acceleration and velocity to predict the event times based on a range of initial positions, and then choose the initial position that minimizes the sum of the squares of the differences between the estimated and measured event times.

t_pred = t_pass; % create array of correct size, the values will be replaced
x_init = (-5:5)*d2/240;
sse_val = x_init*0; % again, creating a correctly sized array
for j = 1:length(x_init)
for i = 1:8
t_pred(i) = 2*(x_event(i)-x_init(j))/(v_0 + sqrt(v_0^2 + 2*a*(x_event(i)-x_init(j))));
end
sse_val(j) = sum((t_pred-t_pass).^2);
end
plot(x_init, sse_val, 'x-')
title('First attempt at finding best initial condition')
xlabel('Initial position (m)')
ylabel('RSS error')
p = polyfit(x_init, sse_val, 2)
x_best = -p(2)/p(1)/2
min(sse_val)
best_sse_est = polyval(p, x_best)
p =
15.1811 -0.0542 0.0000
x_best =
0.0018
ans =
8.1425e-007
best_sse_est =
2.1171e-007

## Now find an improved estimate of the initial position

for i = 1:8
t_pred(i) = 2*(x_event(i)-x_best)/(v_0 + sqrt(v_0^2 + 2*a*(x_event(i)-x_best)));
end
best_sse_calc = sum((t_pred-t_pass).^2)
plot(x_init, sse_val-polyval(p, x_init), 'x-', x_best, best_sse_est, 'o')
x_init = (-5:5)*d2/960+x_best;
sse_val = x_init*0;
for j = 1:length(x_init)
for i = 1:8
t_pred(i) = 2*(x_event(i)-x_init(j))/(v_0 + sqrt(v_0^2 + 2*a*(x_event(i)-x_init(j))));
end
sse_val(j) = sum((t_pred-t_pass).^2);
end
plot(x_init, sse_val, 'x-')
title('Refined attempt at finding best initial condition')
xlabel('Initial position (m)')
ylabel('RSS error')
best_sse_calc =
2.2420e-007

## Use a minimization routine to find best x_0, v_0 and a

The calculations above suggest that finding values of initial position, initial velocity and acceleration that minimize the sum of the squares between modeled and measured times is feasible. Finding an analytical expression for the minimizing values would be challenging, but a minimization algorithm can find these values quickly.

v = version; % for checking if it's running in Matlab or Octave

## Loop through all passes along track

x_best_list = [];
v_best_list = [];
a_best_list = [];
e_best_list = [];
sse_list = [];
x_fwd = x_event;
x_back = -x_event(8:-1:1);
for i = 1:n_trips
i_offset = (i-1)*16;
% forward trip
t_pass = t_meas(i_offset+1:i_offset+8);
t_pass = t_pass-mean(t_pass);
t_mid = (t_pass(1:2:end) + t_pass(2:2:end))/2;
v_est = l_fin./(t_pass(2:2:end)-t_pass(1:2:end));
[v_fit, v_unc, a, v_0, a_unc, v_0_unc] = uncertainFit(t_mid, v_est);
if v(1) == '7' % it's Matlab
[X_best, sse_best] = fminsearch(@(x) sse(x, x_fwd, t_pass), [0 v_0 a]);
% Nelder-Mead unconstrained nonlinear minimization
elseif v(1) == '3' % it's Octave
[X_best, sse_best, info, iter, nf, lambda] = sqp ([0 v_0 a]', @(x) sse(x, x_fwd, t_pass));
% sequential quadratic programming algorithm
else
warning('Unknown environment')
end
x_best_list = [x_best_list X_best(1)];
v_best_list = [v_best_list X_best(2)];
a_best_list = [a_best_list X_best(3)];
e_best_list = [e_best_list sse_best];
% backward trip
t_pass = t_meas(i_offset+9:i_offset+16);
t_pass = t_pass-mean(t_pass);
t_mid = (t_pass(1:2:end) + t_pass(2:2:end))/2;
v_est = l_fin./(t_pass(2:2:end)-t_pass(1:2:end));
[v_fit, v_unc, a, v_0, a_unc, v_0_unc] = uncertainFit(t_mid, v_est);
if v(1) == '7' % it's Matlab
[X_best, sse_best] = fminsearch(@(x) sse(x, x_back, t_pass), [0 v_0 a]);
% Matlab uses Nelder-Mead unconstrained nonlinear minimization
elseif v(1) == '3' % it's Octave
[X_best, sse_best, info, iter, nf, lambda] = sqp ([0 v_0 a]', @(x) sse(x, x_back, t_pass));
% Octave uses a sequential quadratic programming algorithm
else
warning('Unknown environment')
end
x_best_list = [x_best_list X_best(1)];
v_best_list = [v_best_list X_best(2)];
a_best_list = [a_best_list X_best(3)];
e_best_list = [e_best_list sse_best];
end % of the n_trips/2 round trips along the track

## The Results

The main goal was to determine the timing uncertainty of the photogate measurements, but the multiple passes at decreasing speeds allow us to make an estimate of the resistance force as a function of glider velocity. Before doing those calculations, let’s look at the results.

The optimization algorithm provided the resulting sum of the squares of the timing discrepancies, so these provide the information needed to find the timing uncertainty. The number of data points (8 per one-way pass) minus the number of fitted parameters (3: position, velocity and acceleration) gives the number of degrees of freedom to divide the sum by.

t_unc = sqrt(sum(e_best_list)/(5*n_trips*2)) % the uncertainty estimate
m = 0.1478; % kg (mass of glider)
m_unc = 0.00005; % kg (uncertainty of mass measurements)
F = -m*a_best_list;
plot(v_best_list, F, '+-')
title('Analyzing the Resistance Force')
xlabel('Velocity (m/s)')
ylabel('Force (N)')
t_unc =
0.0013

This is a disturbing result. It looks like the forward and backwards trips have different acceleration vs velocity profiles. The likely cause of this is a tilted track; while we made an effort to level the track, it apparently wasn’t as level as we had hoped. To estimate the acceleration bias, take the average of the acceleration differences between the forward and backward passes:

F_bias = mean(F(1:2:end)-F(2:2:end))/2
F_adj = F - F_bias*[1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1];
F_bias =
6.1755e-004

The results are close enough that we can make a plausible force vs. velocity fit. Theoretically, aerodynamic drag results in a quadratic dependence on velocity, while viscosity in the layer of air between the glider and track would give a linear relationship between force and velocity. Let’s check both.

p_quad = polyfit(v_best_list, F_adj, 2);
v_quad = (0:0.05:1)*max(v_best_list);
F_quad = polyval(p_quad, v_quad);
[F_lin, F_unc, C_v, F_0, C_v_unc, F_0_unc] = uncertainFit(v_best_list, F_adj)
F_0/F_0_unc
plot(v_best_list, F_adj, '+-', v_quad, F_quad, '--', v_quad, C_v*v_quad+F_0, ':')
legend('Data', 'Quadratic fit', 'Linear fit', 'Location', 'Northwest')
title('Force after removing bias')
xlabel('Velocity (m/s)')
ylabel('Force (N)')
F_lin =
Columns 1 through 8
0.0018 0.0018 0.0017 0.0016 0.0015 0.0015 0.0014 0.0013
Columns 9 through 16
0.0013 0.0012 0.0011 0.0011 0.0010 0.0010 0.0009 0.0009
Columns 17 through 22
0.0008 0.0008 0.0007 0.0007 0.0006 0.0006
F_unc =
1.9712e-004
C_v =
0.0027
F_0 =
-1.0732e-004
C_v_unc =
2.9286e-004
F_0_unc =
1.4411e-004
ans =
-0.7447

The quadratic fit between velocity and force shows a substantial force at zero velocity, which seems unlikely. The linear fit nearly goes through zero, and in fact the uncertainty in the zero velocity force is greater than the magnitude of the zero velocity force, suggesting that the resistance force is a simple multiple of the velocity.

Published with MATLAB® 7.1