% Corso di Meccanica applciata alle macchine 
% Sapienza Università di Roma 
% Esercitazione n. 12
% Ulteriori riferimenti
% http://dma.ing.uniroma1.it/users/m_appl1_c1/index.html
% 
clear all
close all
%
% Simulazione di un elementare oscillatore libero smorzato
%
%Dati del problema
%
omega = 0.2 % rad/s
zita = 2
m = 10 % kg
x0 = 0  % m
v0 = 0.01 % m/s
xp0 = v0;
%
%  Scelta dell'intervallo di integrazione per il 
% metodo numerico
%
DeltaT =input("Digita il valore di Delta t (passo di integrazione)");
%
% Calcololo degli esponenti
%
alpha = zita* omega
betA = omega*sqrt(zita^2 - 1)
alpha1 =   -alpha + betA
alpha2 =    -alpha - betA 
%
% Calcolo delle ampiezze
%
X1 = v0/(2*betA)
X2=-X1
%
% Tempo di osservazione
%
interT =input("Digita il valore del tempo max (esempio 10, 20 o 30 s)");
t = 0:DeltaT:interT;
%
%  Soluzione analitica
%
x = X1*exp(alpha1 .* t)+X2*exp(alpha2 .* t);
%
%
figure(1)
set(1,'Color','w');
plot(t, x,'-r','LineWidth',1);
%
title(" Soluzione Analitica ");
ylabel("  Spostamento della massa ");
xlabel("  tempo  ");
grid "minor" on;
hold off
% 
% Soluzione numerica
% Metodo elementare
%
N = length(t)
XN = zeros(1,N);  % inizializzo il vettore degli spostamenti numerici
VN = XN;  % inizializzo il vettore delle velocità nimeriche
%
% 
%
XN(1) = x0;
VN(1) = v0;
for i = 1:(N-1)
	xpp = - 2 * zita * omega * VN(i) - XN(i)*omega^2;
	VN(i+1) = VN(i) + DeltaT * xpp;
	XN(i+1) = XN(i) + DeltaT*VN(i);
endfor;
%
% Confronto tra i due metodi
hold on
plot(t,XN,'g','LineWidth',1);
legend ("Soluzione Analitica","Soluzione Numerica");
hold off
%
% confronta tu le velocità analitica con quelle otttenute numericamente
% e raprpesentale in un grafico
%
%  Fai le prove con diversi valori del DeltaT !!!!
%  e .... buon divertimento