function [t u] = sdof_damped(m,k,xi,u0,v0,duration,plotflag)
% This function returns the free vibration displacement of a damped 
% SDOF system with parameters:
% m - mass, kg
% k - stiffness, N/m
% xi - damping ratio
% u0 - initial displacement, m
% v0 - initial velocity, m/s
% duration - length of time of required response
% plotflag - 1 or 0: whether or not to plot the response
% This function returns:
% t - the time vector at which the response was found
% u - the displacement vector of response
% Written by Dr Colin Caprani - www.colincaprani.com

Npts = 1000;    % compute the response at 1000 points
delta_t = duration/(Npts-1);

w = sqrt(k/m);              % rad/s - circular natural frequency
wd = w*sqrt(1-xi^2);        % rad/s - damped circular frequency
ro = sqrt(u0^2+((v0+xi*w*u0)/wd)^2);    % m - amplitude of vibration
theta = atan((v0+u0*xi*w)/(u0*w));      % rad - phase angle

t = 0:delta_t:duration;
u = ro*exp(-xi*w.*t).*cos(w*t-theta);

if(plotflag == 1)
    plot(t,u);
    xlabel('Time (s)');
    ylabel('Displacement (m)');
end