function [f] = drop_rhs(u, p)

% Returns rhs of particle equation given the current state (u)
% and the parameters of the problem (p)

% Unpack the parameters from p
rhog = p(1); % atmospheric density
a    = p(2); % particle radius
mu   = p(3); % atmospheric viscosity
g    = p(4); % gravity
mp   = p(5); % particle mass

% Calculate drag
Re = rhog*u*(2*a)/mu;
if (abs(Re) > 0),
  CD = 24/Re + 6/(1+sqrt(Re)) + 0.4;
  D = 0.5*rhog*u^2*pi*a^2*CD;
else
  D = 0;
end

% Return rhs (total force/particle mass)
f = g - D/mp;