% Codigo Matlab: "OSC"

%-----------------------------------------------------------------------
% Parametros libres:
Ntau    = 100    ;% Numero de pasos temporales;
k	= 1. 	  ;% constante elastica del resorte (en Newton/metro);
l0	= 0. 	  ;% longitud natural del resorte (en metros);
m	= 1.	  ;% masa (en kg);

V0      =  0.	;% Velocidad inicial.;
Z0      =  2.   ;% Posicion inicial (en metros);

ga_w0   = 0.	;% gamma/w_0;
Noscil  = 3.	;% Numero de oscilaciones del experimento.;
%-----------------------------------------------------------------------

   w0 = sqrt(k/m)	     ;% frecuencia natural del resorte
   ga = ga_w0 * w0           ;% constante de amortiguamiento gamma
Omega = sqrt(abs(ga^2-w0^2)) ;% frecuencia de oscilacion

Tau_max = Noscil * (2.*pi) / Omega         ;% Tiempo maximo
Tau     = Tau_max*(0:Ntau-1)/(Ntau-1)      ;% Array de tiempos
epsilon = Tau_max/(Ntau-1)                 ;% paso temporal

  Z=zeros(1,Ntau);  % Array de posicion
  V=zeros(1,Ntau);  % Array de velocidad
  A=zeros(1,Ntau);  % Array de aceleracion
Fnc=zeros(1,Ntau);  % Array de fuerza no conservativa
Wnc=zeros(1,Ntau);  % Array de trabajo no conservativo
Wc =zeros(1,Ntau);  % Array de trabajo conservativo
Fc =zeros(1,Ntau);  % Array de fuerza conservativa

% Primer punto de grilla
    Z(1) = Z0;			% posicion  inicial
    V(1) = V0;			% velocidad inicial
  Fnc(1) = -2.*ga*m* V0; 	% Fuerza no conservativa a t=0
   Fc(1) = -k*(Z0-l0);	        % Fuerza conservativa    a t=0
    A(1) = (Fc(1)+Fnc(1))/m;	% Aceleracion inicial

% Evolucion temporal:
for it=2:Ntau,
	Z(it) = Z(it-1) + V(it-1) * epsilon;
	V(it) = V(it-1) + A(it-1) * epsilon;
      Fnc(it) = -2.*ga*m* V(it);
      Wnc(it) = Wnc(it-1) + Fnc(it-1) * ( Z(it) - Z(it-1));
       Fc(it) = -k*(Z(it)-l0);
        A(it) = (Fc(it)+Fnc(it))/m;
end

% Energias:
T=0.5*m*(V.^2);
U=0.5*k*(Z.^2);
E=T+U;

% Solucion Analitica para comparar:
Zeq = l0;

if ga < w0
if abs(Z0-Zeq) > 0
BB   = atan(-(V0/(Z0-Zeq)+ga)/Omega);
AA   = (Z0-Zeq)/cos(BB);
else
BB   =  pi/2.;
AA   = -V0/Omega;
  end
Z_analitic = ( AA * exp(-ga*Tau) ) .* (cos(Omega*Tau+BB))   + Zeq;
V_analitic = (-AA * exp(-ga*Tau) ) .* (ga*cos(Omega*Tau+BB) + Omega*sin(Omega*Tau+BB));
end

ga1 =0.;
ga2 =0.;
if ga > w0
ga1 = ga+Omega;
ga2 = ga-Omega;
DD = (ga1*(Z0-Zeq)+V0)/(ga1-ga2);
CC = Z0-Zeq - DD;
Z_analitic =  CC * exp(-ga1*Tau)       +DD * exp(-ga2*Tau) + Zeq;
V_analitic = -DD * exp(-ga1*Tau) * ga1 -DD * exp(-ga2*Tau) * ga2;
end

E_analitic = 0.5 * m * V_analitic.^2 + 0.5 * k * (Z_analitic-l0).^2;

break
% Graficos:

% Z(t) numerica (roja) versus analitica (azul)
plot(Tau,Z,'r','linewidth',2)
hold on
plot(Tau,Z_analitic,'b','linewidth',2)
hold off

% V(t) numerica (roja) versus analitica (azul)
plot(Tau,V,'r','linewidth',2)
hold on
plot(Tau,V_analitic,'b','linewidth',2)
hold off

% Sol. Numerica: V(t) (roja); Z(t) (azul)
plot(Tau,Z,'linewidth',2) % Z(t)
hold on
plot(Tau,V,'r','linewidth',2) % V(t)
plot(Tau,zeros(1,Ntau),'g')
hold off

% Espacio fases analitico (rojo) y numerico (azul)
plot(Z,V,'r','linewidth',2)
hold on
plot(Z_analitic,V_analitic,'b','linewidth',2)
hold off

% Energias
plot(Tau,E,'b') % E(t)
hold on
plot(Tau,T,'r') % K(t)
plot(Tau,U,'g') % U(t)
% Variacion de energia y trabajo no-conservativo
plot(Tau,(E-E(1))/E(1))  %  (E-E0)/E0
plot(Tau,Wnc/E(1),'r')       %   Wnc  /E0
hold off

% Verificacion balance energetico:  D(E) - Wnc = 0
plot(Tau,(E-E(1)-Wnc)/E(1),'b') % (E(t)-E0-Wnc)/E0