clc
close all
clear all

% Sección Biofísica y Biología de Sistemas, Facultad de Ciencias, UdelaR
% Taller de Modelización Matemática y Computación para Biociencias (2021)

% OSCILADOR ARMÓNICO
% ------------------------------------------------------------------------------

% Introducción al Modelo
disp('Sección Biofísica y Biología de Sistemas, Facultad de Ciencias, UdelaR')
disp('Taller de Modelización Matemática y Computación para Biociencias (2021)')
disp('-------------------------------------------------------------------------')
disp('OSCILADOR ARMÓNICO')
disp('-------------------------------------------------------------------------')
disp(' ')
disp('Ecuaciones del modelo')
disp(' ')
disp('Posición')
disp('x(t)=2*exp(-gamma*t)*(U*cos(fi*t)+W*sin(fi*t))')
disp(' ')
disp('Velocidad')
disp('y(t)=2*exp(-gamma*t)*(fi*W*cos(fi*t)-fi*U*sin(fi*t))-2*gamma*exp(-g*t)*(W*sin(fi*t)+U*cos(fi*t))')
disp(' ')

% Parámetros del modelo
t = linspace(0,30,250);      % tiempo
fi =input('Ingrese un valor para la Frecuencia Angular de la Oscilación:  ')
gamma = input('Ingrese un valor para el Factor de Amortiguamiento:  ') % factor de amortiguamiento 
% (gamma = 0 sin fricción, gamma > 0 para amortiguamiento, gamma < 0 para amplificación)
% Ej.: con fricción gamma = 1/8
% Ej.: con incremento gamma = -1/8
x0 = input('Ingrese un valor para la posición inicial (x0):  ')
y0 = input('Ingrese un valor para la velocidad inicial (y0): ')

U = x0/2;  
W = (y0+gamma*x0)/(2*fi); 

% Posición vs. Tiempo
x = 2*exp(-gamma*t).*(U*cos(fi*t)+W*sin(fi*t));
figure(1)
hold on
plot(t,x)
if gamma ~= 0
  r=U*2*exp(-gamma*t); % Envolventes de la Oscilación 
  plot(t,r,'r')
  plot(t,-r,'r')
end
hold off
grid on
title('Posición vs. Tiempo')
xlabel('Tiempo (s)')
ylabel('Posición (m)')

% Velocidad vs. Tiempo
y = 2*exp(-gamma*t).*(fi*W*cos(fi*t)-fi*U*sin(fi*t))-2*gamma*exp(-gamma*t).*(W*sin(fi*t)+U*cos(fi*t));
figure(2)
hold on
plot(t,y)
hold off
grid on
title('Velocidad vs. Tiempo')
xlabel('Tiempo (s)')
ylabel('Velocidad (m/s)')

% Velocidad vs. Posición
figure(3)
plot(x,y)
grid on
title('Espacio de fases')
xlabel('Posición (m)')
ylabel('Velocidad (m/s)')