% Studio del sistema massa-molla-smorzatore:
% evoluzione libera, forzata (ingresso a gradino) e completa per i casi: 
%
% 1. attrito elevato (autovalori reali distinti)
% 2. attrito medio (autovalori reali coincidenti)
% 3. attrito moderato (autovalori complessi coniugati)
% 4. attrito nullo (autovalori immaginari puri)

close all
clear all

% parametri del sistema
m=3;                % massa
k=10;               % costante elastica
d=[3*sqrt(k*m) 2*sqrt(k*m) sqrt(k*m) 0];     %coefficiente di attrito
dtext=["elevato","medio","moderato","nullo"];

% EVOLUZIONE LIBERA

% condizione iniziale
x0=[0.1; 0];

% calcolo e plot evoluzione libera per i vari casi
for i=1:4,
 A=[0 1;-k/m -d(i)/m];    % matrice dinamica
 B=[0;1/m];               % matrice ingressi
 C=[1 0];                 % matrice uscite
 D=0;                     % non c'e' legame diretto ingresso/uscita 
 disp('autovalori');disp(eig(A))
 [y,x,t]=initial(A,B,C,D,x0);
 figure(i)
 plot(t,y,'r',t,x(:,2),'b','LineWidth',2)
 title({['y = x_1 (rosso) e x_2 (blu) in evoluzione libera'];['attrito ',dtext{i}]})
 grid
 pause
end

close all

% RISPOSTA FORZATA (ingresso a gradino)

% calcolo e plot risposta forzata (ingresso a gradino) per i vari casi
close all
for i=1:4,
 A=[0 1;-k/m -d(i)/m];    % matrice dinamica
 B=[0;1/m];               % matrice ingressi
 C=[1 0];                 % matrice uscite
 D=0;                     % non c'e' legame diretto ingresso/uscita
 [y,x,t]=step(A,B,C,D);
 figure(i)
 plot(t,y,'r','LineWidth',2)
 title({['y = x_1 (rosso) in risposta forzata (ingresso a gradino)'];['attrito ',dtext{i}]})
 grid
 pause
end

close all

% RISPOSTA COMPLETA (ingresso a gradino)

% calcolo e plot risposta completa (ingresso a gradino) per il
% solo caso 3. (attrito moderato)

close all 
A=[0 1;-k/m -d(3)/m];    % matrice dinamica
B=[0;1/m];               % matrice ingressi
C=[1 0];                 % matrice uscite
D=0;                     % non c'e' legame diretto ingresso/uscita
t=0:0.01:15;             % istanti di tempo per l'integrazione
u=1+0*t;                 % ingresso a gradino unitario
x01=[0.1; -0.1];         % tre diverse condizioni iniziali
x02=[0.0; 0.1];
x03=[-0.05; -0.2];
[y1,x1]=lsim(A,B,C,D,u,t,x01);
[y2,x2]=lsim(A,B,C,D,u,t,x02);
[y3,x3]=lsim(A,B,C,D,u,t,x03);
plot(t,y1,t,y2,t,y3,'LineWidth',2), grid
title('risposta completa per tre diverse condizioni iniziali, attrito moderato')
pause
 
% si noti come il valore di regime sia:
%  - indipendente dalla condizione iniziale (per definizione)
%  - pari al guadagno della funzione di trasferimento

guadagno=dcgain(A,B,C,0)

% RISPOSTA FORZATA (ingresso sinusoidale)

% calcolo e plot risposta forzata (ingresso sinusoidale con pulsazione
% pari a quella naturale) 

close all 
for i=1:4,
 A=[0 1;-k/m -d(i)/m];    % matrice dinamica
 B=[0;1/m];               % matrice ingressi
 C=[1 0];                 % matrice uscite
 D=0;                     % non c'e' legame diretto ingresso/uscita
 t=0:0.01:15;             % istanti di tempo per l'integrazione
 omega_n=sqrt(k/m);
 u=sin(omega_n*t);        % ingresso sinusoidale
 x0=[0;0];                % condizioni iniziali nulle
 [y,x]=lsim(A,B,C,D,u,t,x0);
 figure(i)
 plot(t,y,'r','LineWidth',2)
 title({['risposta forzata (ingresso sinusoidale)'];['attrito ',dtext{i}]})
 grid
 pause
end

% diag di Bode

close all
for i=1:4,
 A=[0 1;-k/m -d(i)/m];    % matrice dinamica
 B=[0;1/m];               % matrice ingressi
 C=[1 0];                 % matrice uscite
 D=0;                     % non c'e' legame diretto ingresso/uscita
 omega_n=sqrt(k/m);
 [numW,denW]=ss2tf(A,B,C,D);
 modW_n=bode(numW,denW,omega_n)
 figure(i)
 bode(numW,denW)
 %plot(t,y,'r','LineWidth',2)
 title({['risposta armonica'];['attrito ',dtext{i}]})
 pause
end



