% Motor de corrente contínua

R = 1;
b = 1;
c = 1.5;
K = 0.5;
L = 1;
Jc = 0.5;
Jm = 0.2;
s = tf('s');

G = c*K/(((Jc+Jm*c^2)*s+b)*(R+L*s)+K^2*c^2);
G = zpk(G)

% Para atingir saída com regime de 100 rad/s
V0 = 100/0.48;
step(V0*G)

C = 0.25/s;
F = feedback(C*G,1)
r = 100
step(r*F)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Pêndulo Invertido 
% Malha aberta
M = 10; m = 2; ell = 1; g = 9.8;
K = [10 50 -300 -100];
A = [0 1 0 0; 0 0 (m*g)/M 0; 0 0 0 1; 0 0 (M+m)*g/(M*ell) 0];
B = [0;1/M ;0; 1/(M*ell)];

% Malha fechada
Af = A+B*K;
x0 = [0 0 -pi/4 0]';
D = 0;

% Posição
Cx = [1 0 0 0];
Fx = tf(ss(Af,x0,Cx,0))
[xt,t] = impulse(Fx);
figure
plot(t,xt)

% Ângulo phi
Ct = [0 0 1 0];
Ft = tf(ss(Af,x0,Ct,0))
[th,t] = impulse(Ft);
figure
plot(t,(th+pi/2)*180/pi)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Sobrelevação versus coeficiente de amortecimento (em %). 

xi = 0:0.01:1;
psi = exp(-xi*pi./(sqrt(1-xi.^2)))*100;
plot(xi,psi)

% Tempo de estabilização.
% y(t) = 1-(exp(-xi*wn*t)/(sqrt(1-xi^2))).*cos(wd*t-atan(xi/sqrt(1-xi^2)));
close all
s = tf('s');
xi = 0.8;
wn = sqrt(5);

F = wn^2/(s^2+2*xi*wn*s+wn^2);

tfinal = 14

t = 0:0.01:tfinal;
ys = 1+(exp(-xi*wn*t)/(sqrt(1-xi^2)));
yi = 1-(exp(-xi*wn*t)/(sqrt(1-xi^2)));

plot(t,ys,'-k')
hold on 
plot(t,yi,'-k')
step(F)

plot([0 tfinal],[1.02 1.02],'-r')
plot([0 tfinal],[0.98 0.98],'-r')

% Aproximação 
yas = 1+(exp(-xi*wn*t));
yai = 1-(exp(-xi*wn*t));

plot(t,yas,'-g')
hold on 
plot(t,yai,'-g')



