function out=waterwheel(t,all_vector,p)

% extract ODE parameter values from input vector
m         = all_vector(1:p.N);  % grab mass vector from input
theta1    = all_vector(p.N+1);  % theta for cup #1
theta     = p.thetadist+theta1; % map theta1 onto all cups
omega     = all_vector(end);    % grab omega vector from input

% mass as a function of time (diff eq)
dmdt      = real(-p.a*sqrt(2*p.g*p.rho*m/pi/p.rc^2) ...
    + p.Q*sin(theta/2).^p.large_n.*(m<p.m_max));

% theta as a function of time (diff eq)
dtheta1dt = omega;

% omega as a function of time (diff eq)
domegadt  = (p.Rw*p.g*sin(-p.phi_0)*sum(m.*sin(theta))-p.alpha*omega ...
    - omega*p.Rw^2*sum(dmdt))/(sum(m)*p.Rw^2);

% combine output for ode output
out       = [dmdt;dtheta1dt;domegadt];
