clear all
close all
clc

% define solution parameters:
f0=1;
T=5;
dt=0.1;

M=1;
C=0.2;
K=1;

beta=0.25;
gamma=0.5;

% i.c. on disp., velocity & acc. respectively:
d(1)=0; 
v(1)=0;
a(1)=0;

% time axis:
t=[0:dt:4*T];

% define force vector:
for i=1:length(t)
    if t(i)<=T/2
        f(i)=2*f0/T*t(i);
    elseif t(i)<=T
        f(i)=-2*f0/T*(t(i)-T);
    else
        f(i)=0;
    end
end

% solve with Newmark:
for i=2:length(t)
    
    % define predictors:
    dpred(i) = d(i-1) + dt*v(i-1) + dt^2/2*(1-2*beta)*a(i-1);
    vpred(i) = v(i-1) + (1-gamma)*dt*a(i-1);
    
    % solve for acc.: 
    a(i) = ( K*f(i) - C*vpred(i) - K*dpred(i) ) / ( M + gamma*dt*C + beta*dt^2*K );
 
    % evaluate disp. & velocity:
    d(i) = dpred(i) + beta*dt^2*a(i); 
    v(i) = vpred(i) + gamma*dt*a(i);
    
end

plot(t,d,'b')

hold on

% evaluate analytic solution vector:

t2 = t - T/2;
t3 = t - T;

for i=1:length(t)
  
    r1(i)=t(i)/K - 2*xi/K/wn + exp(-xi*wn*t(i)) * ( (2*xi^2-1)*sin(wd*t(i))/wd/K + 2*xi*cos(wd*t(i))/wn/K );
    
    r2(i)=t2(i)/K - 2*xi/K/wn + exp(-xi*wn*t2(i)) * ( (2*xi^2-1)*sin(wd*t2(i))/wd/K + 2*xi*cos(wd*t2(i))/wn/K );
    
    r3(i)=t3(i)/K - 2*xi/K/wn + exp(-xi*wn*t3(i)) * ( (2*xi^2-1)*sin(wd*t3(i))/wd/K + 2*xi*cos(wd*t3(i))/wn/K );
    
    x_an(i) = 2*K*f0/T * ( r1(i) - 2*r2(i)*(t(i)>=T/2) + r3(i)*(t(i)>=T) );
end

plot(t,x_an,'r')