//===============================================================
// Scilab como herramienta en Sistemas y Control
//===============================================================

//===============================================================
// PARTE 1
//===============================================================

//
// Definiendo una funcion de transferencia en general
//
// Funcion de transferencia de lazo abierto del sistema
// H(s) =   1/ s(s+2)(s+4)
//

// Definiendo H en funcion de sus polos

H1 = poly([ 0 -2 -4],"s", ["roots"]) 

H = 1/H1

pause

// Definiendo H en forma simbolica

s=poly(0,"s") // Polinomio generico s

H =   1/ (s*(s+2)*(s+4))

pause

// Otra forma
H1 =   1/ s

H2=   1/ (s+2)

H3=   1/(s+4)

H=H1*H2*H3 

pause

// Definiendo el sistema lineal cuya funcion de transferencia es H

[sist_test]=syslin('c',H)   

// Hallando el LGR (Lugar Geometrico de las Raices)

evans(sist_test)

pause

evans(sist_test,100)
pause


g=kpure(sist_test) // Maxima ganancia para tener raices sobre el eje imaginario
evans(sist_test,max(g))

pause

//===============================================================
// Respuesta en el dominio de la frecuencia
//===============================================================


//===============================================================
// Diagrama de Bode 
//===============================================================

[ frq repf]=repfreq(sist_test,0.001, 1000)

[db,phi]=dbphi(repf)

subplot(2,1,1)
plot2d(frq,db,logflag="ln")  
xtitle('Magnitud')

subplot(2,1,2)
plot2d(frq,db,logflag="ln")  
xtitle('Angulo')

pause


bode(sist_test)

pause

//===============================================================
// Diagrama de Nyquist 
//===============================================================

Ro = abs(repf);
X = real(repf);
Y= imag(repf);

subplot(2,1,1)
plot2d(X,Y);
xtitle('Diagrama de Nyquist')
xgrid

subplot(2,1,2)
nyquist(frq, repf,'Directamente')



//===============================================================
// PARTE 2 - VARIABLES DE ESTADO Y SIMULACIONES EN EL TIEMPO
//===============================================================


//
// Definiendo un sistema sabiendo su representacion de estado
//

A = [ -0.7 -0.05; 1 0]
B = [  1 0]'
C = [  0  1]
D = [ 0]

[sist_1]=syslin('c',A,B,C,D)   

pause


//
// Verificamos la correcta entrada de las matrices de estado
//

sist_1("A")
sist_1("B")
sist_1("C")
sist_1("D")
sist_1("X0")   // Valor inicial

pause

//
// Hallando la funcion de transferencia
//

[h]=ss2tf(sist_1) 

Num = numer(h)

Den = denom(h)


//
// Caso contrario de querer pasas de funcion de transferencia a 
// representacion de estado usar tf2ss
//

pause

//
// Autovalores (polos) del sistema 
//

spec(A)

// Verificamos polos

roots(Den)

pause

//
// Simulacion en el tiempo
//

t=0:0.05:100

// Respuesta a una funcion escalon

U = ones(1,length(t));

[y1 X]=csim(U,t,sist_1);

plot2d(t, [U', y1'])
//legend(['u(t)';'Salida(t)'],a=2) Corregir esto


pause

xbasc

// Respuesta a una funcion escalon

[y2 X]=csim("step",t,sist_1);
plot2d(t, y2)

pause

xbasc


// Respuesta a una funcion senoidal

U = sin(t);
[y3 X]=csim(U,t,sist_1);
plot2d(t, [y3' X(1,:)' X(2,:)'], style=[color("red"),color("green"),color("blue")])

// Otras herramientas

// Matriz de controlabilidad

Cc=cont_mat(A,B)

pause

Cc=cont_mat(sist_l)

pause

// Matriz de observabilidad
Obs = obsv_mat(A,C)   


[O]=obsv_mat(sl)