function [t u] = sdof_forced(m,k,xi,u0,v0,F,Omega,duration,plotflag)
% This function returns the displacement of a damped SDOF system which is 
% subjected to harmonic excitation with parameters:
% m - mass, kg
% k - stiffness, N/m
% xi - damping ratio
% u0 - initial displacement, m
% v0 - initial velocity, m/s
% F - amplitude of forcing function, N
% Omega - frequency of forcing function, rad/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

beta = Omega/w;         % frequency ratio
D = DAF(beta,xi);       % dynamic amplification factor
ro = F/k*D;             % m - amplitude of vibration
theta = phase(beta,xi); % rad - phase angle

% Constants for the transient response
Aconst = u0+ro*sin(theta);
Bconst = (v0+u0*xi*w-ro*(Omega*cos(theta)-xi*w*sin(theta)))/wd;

t = 0:delta_t:duration;
u_transient = exp(-xi*w.*t).*(Aconst*cos(wd*t)+Bconst*sin(wd*t));
u_steady = ro*sin(Omega*t-theta);
u = u_transient + u_steady;

if(plotflag == 1)
    plot(t,u,'k'); 
    hold on;
    plot(t,u_transient,'k:');
    plot(t,u_steady,'k--');
    hold off;
    xlabel('Time (s)');
    ylabel('Displacement (m)');
    legend('Total Response','Transient','Steady-State');
end