//======|_FUNÇÕES_DE_TRANSFERÊNCIA-RESPOSTA_EM_FREQUÊNCIA_|=======//
clear() //limpa as variáveis anteriores
//----------------------------parâmetros--------------------------//
///parâmetros hidrodinâmicos
Xu = 82.30 //coeficiente linear de arrasto em avanço [kg/s]
Xuu = 309.70 //coeficiente quadrático de arrasto em avanço [kg/m]
Yv = 8.50 //coeficiente linear de arrasto em deriva [kg/s]
Yvv = 505.45 //coeficiente quadrático de arrasto em avanço [kg/s]
Kr = 18.21 //coeficiente linear de arrasto em guinada [kg/s]
Krr = 94.72 //coeficiente quadrático de arrasto em guinada [kg/m]
//parâmetros do ROV
mu = 450.87 //massa virtual de avanço [kg]
mv = 558.92 //massa virtual de deriva [kg]
Jr = 239.44 //inércia virtual de guinada [kg*m^2]
//----------------------velocidades_nominais----------------------//
un = 0.5 //velocidade de avanço [m/s]
vn = 0.5 //velocidade de deriva [m/s]
rn = 0 //velocidade de guinada [rad/s]
//-----------------------sistema_linearizado----------------------//
//d(deltaX)/dt = A*deltaX + B*deltaU
//deltaY = C*deltaX + D*deltaU
//A(nxn): matriz dinâmica
A = [0 , 1 , 0 , 0 , 0 , 0 ;..
0 , -(Xu+2*Xuu*abs(un))/mu , 0 , rn*mv/mu , 0 , vn*mv/mu ;..
0 , 0 , 0 , 1 , 0 , 0 ;..
0 , -rn*mu/mv , 0 , -(Yv+2*Yvv*abs(vn))/mv , 0 , -un*mu/mv ;..
0 , 0 , 0 , 0 , 0 , 1 ;..
0 , vn*(mu-mv)/Jr , 0 , un*(mu-mv)/Jr , 0 , -(Kr+2*Krr*abs(rn))/Jr ]
//B(nxp): matriz de entrada
B = [ 0 , 0 , 0 ;..
1/mu , 0 , 0 ;..
0 , 0 , 0 ;..
0 , 1/mv , 0 ;..
0 , 0 , 0 ;..
0 , 0 , 1/Jr ]
//C(mxn): matriz de saídas
C = [ 0 , 1 , 0 , 0 , 0 , 0 ;..
0 , 0 , 0 , 1 , 0 , 0 ;..
0 , 0 , 0 , 0 , 0 , 1 ]
//D(mxp): matriz de incidência direta
D = [ 0 , 0 , 0 ;..
0 , 0 , 0 ;..
0 , 0 , 0 ]
//---------------------funções_de_transferência-------------------//
//definindo a variável complexa s
s = poly(0,'s')
//matriz resolvente (nxn)
PHI_s = inv(s*eye(A)-A)
//eye(A) retorna a matriz identidade com as dimensões de A
//denominador (polinômio de grau n)
E = det(s*eye(A)-A)
//polos (n polos)
polos = roots(E)
//matriz de transferência (mxp); G(i,j) = Y(i)/U(j)
G = syslin('c',C*PHI_s*B)
//-----------------exibindo_informações_relevantes----------------//
disp('un = '+string(un)+'; vn = '+string(vn)+'; rn = '+string(rn))
disp('E:')
disp(E)
disp('polos:')
disp(polos)
disp('G:')
disp(G)
//------------------------------plots-----------------------------//
xset('window',1)
plzr(G)
//-------------------------diagramas_de_Bode----------------------//
if G(1,1) ~= 0 then
    xset('window',2)
    bode(G(1,1))
    xtitle('Diagramas de Bode: u/F_u')
end
if G(1,2) ~= 0 then
    xset('window',3)
    bode(G(1,2))
    xtitle('Diagramas de Bode: u/F_v')
end
if G(1,3) ~= 0 then
    xset('window',4)
    bode(G(1,3))
    xtitle('Diagramas de Bode: u/T_r')
end
if G(2,1) ~= 0 then
    xset('window',5)
    bode(G(2,1))
    xtitle('Diagramas de Bode: v/F_u')
end
if G(2,2) ~= 0 then
    xset('window',6)
    bode(G(2,2))
    xtitle('Diagramas de Bode: v/F_v')
end

if G(2,3) ~= 0 then
    xset('window',7)
    bode(G(2,3))
    xtitle('Diagramas de Bode: v/T_r')
end
if G(3,1) ~= 0 then
    xset('window',8)
    bode(G(3,1))
    xtitle('Diagramas de Bode: r/F_u')
end

if G(3,2) ~= 0 then
    xset('window',9)
    bode(G(3,2))
    xtitle('Diagramas de Bode: r/F_v')
end
if G(3,3) ~= 0 then
    xset('window',10)
    bode(G(3,3))
    xtitle('Diagramas de Bode: r/T_r')
end