%% prog diferencias finitas

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

clearvars;
close all;
clc;

%% datos inciales

m = 1;        %masa
k = 25;       %resorte
c = 0.0;      %amortiguador viscoso

%% tiempo

dt = 0.01;          %intervalo de muestreo
tini = 0;           %tiempo incial 
tfin = 10;          %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*cos(5*t);

%% 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;   %nincializacion de la velocidad

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

%% resultados

figure(1)
plot(t,f);
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);
title('desplazamiento diferencias finitas');
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);
title('velocidad diferencias finitas');
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);
title('diagrama de fase diferencias finitas');
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;


for n = 1:N

  figure(5)
  plot(x(n),v(n));
  title('diagrama de fase diferencias finitas');
  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;
  hold on;
  
  pause(0.01);
  
end;