function [u ud udd] = newmark_sdof(m, k, xi, t, F, u0, ud0, plotflag)
% This function computes the response of a linear damped SDOF system 
% subject to an arbitrary excitation. The input parameters are:
% m         - scalar, mass, kg
% k         - scalar, stiffness, N/m
% xi        - scalar, damping ratio
% t         - vector of length N, in equal time steps, s
% F         - vector of length N, force at each time step, N
% u0        - scalar, initial displacement, m
% v0        - scalar, initial velocity, m/s
% plotflag  - 1 or 0: whether or not to plot the response
% The output is:
% u         - vector of length N, displacement response, m
% ud        - vector of length N, velocity response, m/s
% udd       - vector of length N, acceleration response, m/s2
% Written by Dr Colin Caprani - www.colincaprani.com

% Set the Newmark Integration parameters
% gamma = 1/2 always
% beta = 1/6 linear acceleration
% beta = 1/4 average acceleration
gamma = 1/2;
beta = 1/6;

N = length(t);  % the number of integration steps
dt = t(2)-t(1); % the time step
w = sqrt(k/m);  % rad/s - circular natural frequency
c = 2*xi*k/w;   % the damping coefficient

% Calulate the effective stiffness
keff = k + (gamma/(beta*dt))*c+(1/(beta*dt^2))*m;
% Calulate the coefficients A and B
Acoeff = (1/(beta*dt))*m+(gamma/beta)*c;
Bcoeff = (1/(2*beta))*m + dt*(gamma/(2*beta)-1)*c;

% calulate the change in force at each time step
dF = diff(F);

% Set initial state
u(1) = u0;
ud(1) = ud0;
udd(1) = (F(1)-c*ud0-k*u0)/m; % the initial acceleration

for i = 1:(N-1) % N-1 since we already know solution at i = 1
    dFeff = dF(i) + Acoeff*ud(i) + Bcoeff*udd(i);
    dui = dFeff/keff;
    dudi = (gamma/(beta*dt))*dui-(gamma/beta)*ud(i)+dt*(1-gamma/(2*beta))*udd(i);
    duddi = (1/(beta*dt^2))*dui-(1/(beta*dt))*ud(i)-(1/(2*beta))*udd(i);
    u(i+1) = u(i) + dui;
    ud(i+1) = ud(i) + dudi;
    udd(i+1) = udd(i) + duddi;
end

if(plotflag == 1)
    subplot(4,1,1)
    plot(t,F,'k'); 
    xlabel('Time (s)');
    ylabel('Force (N)');
    subplot(4,1,2)
    plot(t,u,'k'); 
    xlabel('Time (s)');
    ylabel('Displacement (m)'); 
    subplot(4,1,3)
    plot(t,ud,'k'); 
    xlabel('Time (s)');
    ylabel('Velocity (m/s)');
    subplot(4,1,4)
    plot(t,udd,'k'); 
    xlabel('Time (s)');
    ylabel('Acceleration (m/s2)'); 
end