function x = dho_position( m, c, k, x0, t )

% Find the position of a damped harmonic oscillator spring with zero
% initial velocity:
%   x = x0*w0/w * exp(-g*t) * cos(w*t-phi)
% with
%   w0 = sqrt(k/m)
%   w = w0*sqrt(1-z^2)
%   z = c/2/sqrt(m*k)
%   g = c/2/m
%   phi = acos(sqrt(1-z^2))
% Input:
%   m   : mass (kg)
%   c   : damping constant (kg/s)
%   k   : spring constant (N/m)
%   x0  : initial position (m)
%   t   : time (s)
% J.M.Soler, Oct.2015

w0 = sqrt(k/m);
z = c/2/sqrt(m*k);
w = w0*sqrt(1-z^2);
g = c/2/m;
phi = acos(sqrt(1-z^2));
x = x0*w0/w * exp(-g*t) .* cos(w*t-phi);

end
