%% prog diferencias centrales

%este programa resuelve la EDO de segundo orden con condiciones inciales
%usando diferencias finitas

clearvars;
close all;
clc;

%% datos inciales

m = 1;        %masa
k = 10;       %resorte
c = 0.5;      %amortiguador viscoso

%% tiempo

dt = 0.01;          %intervalo de muestreo
tini = 0;           %tiempo incial 
tfin = 200;         %tiempo final
t = tini:dt:tfin;   %vector de tiempo
N = length(t);      %numero de muestras total

%% condiciones inciales

x0 = 0;     %desplazamiento incial
v0 = 0;     %velocidad incial

%% fuerza externa

f(1:N) = 0;             %incializacion de fuerza
tf1 = 10;               %tiempo1
tf2 = 30;               %tiempo2
tf3 = 50;               %tiempo2
nf1 = round(tf1/dt);    %muestra1
nf2 = round(tf2/dt);    %muestra2
nf3 = round(tf3/dt);    %muestra3
f(nf1:nf2) = 50;
f(nf2+1:nf3) = -50;

%% desplazamiento

%paso 1 determinar la aceleracion

f0 = f(1);
a0 = (f0 - x0 - v0)/m;

%paso 2 determinar el desplazamiento x1

x1 = x0 + v0*dt +(1/2)*a0*dt^2;

%paso 3 determinamos el resultado mediante recurrencia

x(1:N) = 0;   %incializacion de la solucion
x(1) = x0;
x(2) = x1;

C1 = m/(dt^2) + c/(2*dt);
C2 = k - (2*m)/(dt^2);
C3 = m/(dt^2) - c/(2*dt);

for n = 2:N-1
    x(n+1) = ( f(n) - C2*x(n) - C3*x(n-1) )/C1;
end;


%% velocidad

v(1:N) = 0;   %incializacion de la velocidad

for n = 2:N-1
    v(n) = ( x(n+1) - x(n-1) )/(2*dt);
end;

%% resultados

figure(1)
plot(t,f,'b');
title('fuerza');
xlabel('tiempo t(s)');
ylabel('fuerza f(t) (N)');
legend('f(t)');
axis([tini, tfin, -1.1*max(abs(f)), 1.1*max(abs(f))]);
grid on;
box on;


figure(2)
plot(t,x,'b');
title('desplazamiento diferencias centrales');
xlabel('tiempo t(s)');
ylabel('desplazamiento x(t) (m)');
legend('x(t)');
axis([tini, tfin, -1.1*max(abs(x)), 1.1*max(abs(x))]);
grid on;
box on;


figure(3)
plot(t,v,'b');
title('velocidad diferencias centrales');
xlabel('tiempo t(s)');
ylabel('velocidad v(t) (m/s)');
legend('v(t)');
axis([tini, tfin, -1.1*max(abs(v)), 1.1*max(abs(v))]);
grid on;
box on;


figure(4)
plot(x,v,'b');
title('diagrama de fase diferencias centrales');
xlabel('desplazamiento x(t) (m)');
ylabel('velocidad v(t) (m/s)');
legend('x(t) vs v(t)');
axis([-1.1*max(abs(x)), 1.1*max(abs(x)), -1.1*max(abs(v)), 1.1*max(abs(v))]);
grid on;
box on;