Solution to dy/dan's Fan Slowing Problem

FanProblem.m

% Kenneth Finnegan, 2012
% kennethfinnegan.blogspot.com
% 
% Predict end of spinning object
% Solution to http://blog.mrmeyer.com/?p=13245

clear all;
close all;

% Single step through video and pull blade hits on bottom of frame
rawdata = [ 0.36, 0.53, 0.70, 0.88, 1.06 1.23 1.42 ...
    1.61 1.80 1.98 2.16 2.37 2.56 2.76 2.96 3.16 ...
    3.37 3.58 3.79 4.00 4.21 4.43 4.65 4.88 5.10 ...
    5.33 5.57 5.83 6.07 6.31 6.56 6.80 7.06 7.30 ...
    7.57 7.83 8.09 8.35 8.62 8.90 9.17 9.46 9.74 ... 
    10.03 10.32 10.61 10.91 11.21 11.51 11.83 12.14 ...
    12.47 12.81 13.13 13.46 13.81 14.14 14.51 14.83 ...
    15.19 15.54 15.93 16.31 16.71 17.08 17.49 17.90 ...
    18.31 18.74 19.16 19.60 20.05 20.48 20.95 21.42 ...
    21.88 22.36 22.85 23.35 23.86 24.39 24.92 25.45];

% Convert blade hit times to rotational velocity / 4
%       The /4 is because I counted all four blades hitting bottom
%       of the frame, not just one, to get more data
for index = 1:(length(rawdata) - 1)
    rotvel(index) = 1 / (rawdata(index + 1) - rawdata(index));
end

% Clean up the data a little bit due to frame rate aliasing and
%   me being human
for index = 1:(length(rotvel) - 2)
    rotvel(index) = mean(rotvel(index:index+2));
end

figure(1);
plot(rawdata(1:length(rawdata)-1), rotvel);
xlabel('Time [seconds]');
ylabel('Rotational Velocity [blades per second]');

% Take derivative to find acceleration
for index = 1:(length(rotvel) - 1)
    rotaccel(index) = (rotvel(index + 1) - rotvel(index)) / ...
        (rawdata(index + 1) - rawdata(index));
end

figure(2);
plot(rawdata(1:length(rawdata)-2), rotaccel);
xlabel('Time [seconds]');
ylabel('Acceleration [blades / sec / sec]');

figure(3);
hold on;
scatter(rotvel(1:(length(rotvel) - 1)), rotaccel);
xlabel('Rotational Velocity [blades / sec]');
ylabel('Rotational Acceleration [blades /sec / sec');

% Considering that the acceleration (a = F / m) is somehow 
% related to the velocity of the fan, we'll say that it's
% probably a second order polynomial, because air drag is N^2,
% and we might be able to get away with calling friction N + C.
%  (Friction is really probably a mess, and we're ignoring the 
%   motor's induction implications)

accelfromvel = polyfit(rotvel(1:(length(rotvel) - 1)), rotaccel, 2);
plot(1.5:0.1:6, polyval(accelfromvel, 1.5:0.1:6), 'r');

figure(4); hold on;
% Now that we have this acceleration(velocity) = N^2 + N + C model,
% we will take each velocity point we got from the original video
% and extrapolate them to the end
stepsize = 0.01;
for index = 1:(length(rotvel) - 1)
    clear model_v model_t model_index
    model_index = 1;
    model_t(1) = rawdata(index);
    model_v(1) = rotvel(index);
    
    % Given those intial conditions, walk out to the future
    while model_v(model_index) > 0
        model_t(model_index + 1) = model_t(model_index) + stepsize;
        model_v(model_index + 1) = model_v(model_index) + ...
            polyval(accelfromvel, model_v(model_index)) * stepsize;
        model_index = model_index + 1;
    end
    
    endtimes(index) = model_t(length(model_t));
    plot(model_t, model_v, 'r');
end

plot(rawdata(1:length(rawdata)-1), rotvel, 'bd-');
axis([0 (max(endtimes) + 5) 0 6]);
xlabel('Time [seconds]');
ylabel('Rotational Velocity [blades per second]');

figure(5);
hist(endtimes, linspace(min(endtimes), max(endtimes), 10));
xlabel('Modeled End Time [seconds]');
ylabel('Number of simulations [#]');

PredictedFinished = mean(endtimes)
% Result for 0.01 stepsize with 3 point smoothing = 46.7 seconds
%
% Due to aliasing of the data, The model oscillates between ~46 sec
% for odd smoothing and ~62 seconds for even smoothing