function y=risposta_libera(den,y0,T);
% RISPOSTA_LIBERA : traccia la risposta libera di un modello ingresso-uscita
% ******************************************************************************
%                ## SYNTAX ##
%
% RISPOSTA_LIBERA(DEN,Y0): 
%    Traccia l'evoluzione libera di un sistema la cui funzione di trasferimento e'
%                  NUM
%         W(s) =  -----
%                  DEN
%    a partire da condizioni iniziali dell'uscita e delle sue derivate date dal
%    vettore Y0 = [y(0) y'(0) ... y^(n-1)(0)]. La scala dei tempi e' scelta 
%    automaticamente. Si noti che la conoscenza di NUM non e' necessaria.
%    Il valore analitico di questo segnale viene visualizzato nel command window.    
%
% RISPOSTA_LIBERA(DEN,Y0,T): 
%    Usa il vettore dei tempi T.  
%
% Y=RISPOSTA_LIBERA(DEN,Y0): 
% Y=RISPOSTA_LIBERA(DEN,Y0,T):
%     Non traccia il diagramma ma costruisce il vettore Y contenente i valori dell'uscita. 
%
%     Si veda anche LSIM
%
% Alessandro Giua, Novembre 2002
ni=nargin;
no=nargout;

if ((ni <1) & (ni>3)),
    disp('Errore: controllare il numero degli argomenti')
    return
end
if (length(den) ~= length(y0)+1),
    disp('Errore: se il polinomio al denominatore ha n+1 coefficienti Y0 deve avere lunghezza n')
    return
end
if (den(1)==0),
    disp('Errore: il coefficiente di grado maggiore in DEN deve essere non nullo')
    return
end
if (length(y0)==0),
    disp('Il sistema e'' istantaneo')
    return
end

% Scelgo una rappresentazione in variabili di stato che usa lo spazio di fase:
%   x_1 = y; x_2 = y'; ... x_n = y^(n-1)
% vedi libro (7.5)
n=length(den)-1;
for i=1:n,
    alpha(i)=-den(n+2-i)/den(1);
end
if n==1,
    A=alpha(1);
    B=1/den(1);
    C=1;
else
    A=[zeros(n-1,1),eye(n-1); alpha];
    B=[zeros(n-1,1); 1/den(1)];
    C=[1 zeros(1,n-1)];
end
D=0;
sys=ss(A,B,C,D);

% Se T non e' stato definito scelgo la scala dei tempi pari a 6 volte la
% massima costante di tempo
if ni==2,
    p=abs(real(roots(den)));
    [i,j,v]=find(p);
    if length(v)>0,
        tau=1/min(v);
    else 
        tau=1;
    end
    T=0:6*tau/200:6*tau;
end

% Calcolo lo stato inziale
[a b]=size(y0);
if a < b, 
    x0 = y0';
else
    x0 = y0;
end

% calcolo la risposta libera numericamente
if no==0,
    lsim(sys,0*T,T,x0);
else
    y=lsim(sys,0*T,T,x0);
end

% calcolo la risposta libera analiticamente
syms t
y_elle = C*expm(A*t)*x0
return