close all

dt = 0.01;
tau = 0:dt:10;
i = 3;
j = 3;
L = ComputeIRF(w,F,tau);


plot(tau,L{i,j});


T = 2;  %seconds
f = 1/T;


t = 0:dt:30;
tc = 5;  %Ramp up time constant (seconds)
x = (1-exp(-t/tc)).*cos(2*pi*f*t);
xdot = -(1-exp(-t/tc)).*sin(2*pi*f*t)*2*pi*f+cos(2*pi*f*t).*exp(-t/tc)/tc;
xddot = -(1-exp(-t/tc)).*cos(2*pi*f*t)*(2*pi*f)^2 -2*pi*f*sin(2*pi*f*t).*exp(-t/tc)/tc -cos(2*pi*f*t).*exp(-t/tc)/(tc^2) -2*pi*f*exp(-t/tc).*sin(2*pi*f*t)/tc;


figure
plot(t,x,t,xdot,t,xddot);
xlabel('t');
legend('x','xdot','xddot');

Force = -real(F{i,j}(end))*xddot;

for k = 2:length(Force)
NN = min(k,length(tau));
Integrand = zeros(1,NN);
  for n = 1:NN
  Integrand(n) = xddot(k-n+1)*L{i,j}(n);
  end
I = trapz(tau(1:NN),Integrand);
    
Force(k) = Force(k)-I;
end



%Compute Frequency Domain Forces
a = spline(w,real(F{i,j}),2*pi*f);
b = spline(w,imag(F{i,j}),2*pi*f);

ForceFreq = cos(2*pi*f*t)*(2*pi*f)^2*a+sin(2*pi*f*t)*2*pi*f*b;


figure
plot(t,Force,t,ForceFreq)