function [x,v] = oscillator( x0, v0, omega, tau, t )

% Returns the position of a damped harmonic oscillator of the form
%   x(t) = ( A*cos(omega*t) + B*sin(omega*t) ) * exp(-t/tau)
% Input:
%   x0    : initial position (m)
%   v0    : initial velocity (m/s)
%   omega : angular frequency (rad/s)
%   tau   : damping time (s)
%   t     : time(s) at which x is to be found (s)
% Output:
%   x     : position(s) of the oscillator at time t (m)
%   v     : velocity(s) of the oscillator at time t (m/s)
% Usage:
%   [x,v] = oscillator( 1.5, 25., 3.5, 5., 10. )
% Author:
%   J.M.Soler. Sept.2010

% Find A and B such that
%   x0 = x(0) = A
%   v0 = dx/dt(0) = B*omega - A/tau
A = x0;
B = (v0 + A/tau) / omega;

% Find x and v
x = ( A*cos(omega*t) + B*sin(omega*t) ) .* exp(-t/tau);
v =((-A*sin(omega*t) + B*cos(omega*t) ) * omega ...
   -( A*cos(omega*t) + B*sin(omega*t) ) / tau ) .* exp(-t/tau);

end
