%%
% date: October 12 
% Last modified:

%
%%

close all; clear all; clc;

load elcentro.dat;

g = 9810;                % mm/sec^2
m = 1;                   % kN-sec^2/mm
Tn = 0.2;                  % sec
wn = 2*pi/Tn;            %
k = (2*pi/Tn)^2 * m;     % kN/mm
Dt = 0.02;               % sec
z = 0.02;
c = z * 2 * m *wn;       %
uo = 0;                  % initial displacement (t=0)
vo = 0;                  % initial velocity (t=0)
po = -m * elcentro(1,2)*g % force vector (kN)


%% Initial Conditions

% a: relative acceleration

a = (po - c * vo - k *uo)/m;
u_1 = uo - Dt * vo + (Dt)^2 * a/2;
k_h = m/(Dt^2) + c/(2*Dt);
alpha = m/(Dt^2) - c/(2*Dt);
beta = k - 2*m/(Dt^2);


%% Calculations for time step i
Time = 30;           % [sec] of response
Number_of_steps = Time /Dt; 
to = 0;
for i = 1:Number_of_steps
    po = -m*elcentro(i,2)*g;

    p_hut(i) = po - alpha * u_1 - beta * uo;
    
    u(i) = p_hut(i)/k_h;
    
    u_1 = uo;

    uo = u(i);
   
    a(i) = (uo - 2*u(i) + u_1)/(Dt^2);
    v(i) = (u(i) - uo)/(2*Dt);
    t(i) = Dt * i;
end

figure
plot(elcentro(:,1),elcentro(:,2)*g)
xlabel('Time [sec]')
ylabel('Acc [mm/sec^2]')

figure
plot(t, u)
xlabel('Time [sec]')
ylabel('Rel. Displ. [mm]')