%%*************************************************************
% code writer: EOF
% code date:2014.03.18
% e-mail: jasonleaster@gmail.com
% code purpose : 
%           I just want to share with someone who is interesting
% in adaptive control. This code is to help people to understand
% second order system.
%%**************************************************************
clear all
clc
syms s f t m;
K01 = 1;
K02 = 1;
K0 = (K01*K02)./(1+K01*K02);
hold on;
figure(1);
T0 = 1;
for epsilon_0 = 0.1:0.2:2
    T = T0./(1+K0);
%     epsilon_0 = 0.5*(1/(K01*K02*T0));
    epsilon = epsilon_0./(1+K0);
    K = K0/(1+K0);
    f = (K./((T.^2).*(s.^2)+2*epsilon_0.*T.*s+K0)).*(1./s);
    m = ilaplace(f);
    ezplot(m,[0,120]);
    axis([0 60 0 1.2]);
end
legend(‘0.1‘,‘0.3‘,‘0.5‘,‘0.7‘,‘0.9‘,‘1.1‘,‘1.3‘,‘1.5‘,‘1.7‘,‘1.9‘);
hold off;
figure(2);
hold on;
for epsilon_0 = 0.1:0.2:2
    if epsilon_0 <1
        plot(-epsilon_0.*(1./T),(1./T).*sqrt(1-(epsilon_0).^2),‘*‘);
        plot(-epsilon_0.*(1./T),-(1./T).*sqrt(1-(epsilon_0).^2),‘*‘);
    elseif abs(epsilon_0-1.1) < 0.1  
        plot(-epsilon_0.*(1./T)+(1./T).*sqrt((epsilon_0).^2-1),0,‘*‘,‘Color‘,‘r‘);
        plot(-epsilon_0.*(1./T)-(1./T).*sqrt((epsilon_0).^2-1),0,‘*‘,‘Color‘,‘r‘);
    elseif abs(epsilon_0-1.3) < 0.1 
        plot(-epsilon_0.*(1./T)+(1./T).*sqrt((epsilon_0).^2-1),0,‘*‘,‘Color‘,‘g‘);
        plot(-epsilon_0.*(1./T)-(1./T).*sqrt((epsilon_0).^2-1),0,‘*‘,‘Color‘,‘g‘);
    elseif abs(epsilon_0-1.5) < 0.1 
        plot(-epsilon_0.*(1./T)+(1./T).*sqrt((epsilon_0).^2-1),0,‘*‘,‘Color‘,‘b‘);
        plot(-epsilon_0.*(1./T)-(1./T).*sqrt((epsilon_0).^2-1),0,‘*‘,‘Color‘,‘b‘);
    elseif abs(epsilon_0-1.7) < 0.1 
        plot(-epsilon_0.*(1./T)+(1./T).*sqrt((epsilon_0).^2-1),0,‘*‘,‘Color‘,‘y‘);
        plot(-epsilon_0.*(1./T)-(1./T).*sqrt((epsilon_0).^2-1),0,‘*‘,‘Color‘,‘y‘);
    elseif abs(epsilon_0-1.9) < 0.1 
        plot(-epsilon_0.*(1./T)+(1./T).*sqrt((epsilon_0).^2-1),0,‘*‘,‘Color‘,‘k‘);
        plot(-epsilon_0.*(1./T)-(1./T).*sqrt((epsilon_0).^2-1),0,‘*‘,‘Color‘,‘k‘);
    end
end
axis([-6 1 -3 3]);
hold off;
figure(3);
hold on;
Wn = (1./T);
for epsilon_0 = 0.1:0.2:2 
    Wd = Wn.*(1-(epsilon_0).^2);
    if epsilon_0 < 1
        belta = acos(epsilon_0);
        tr = (pi-belta)./Wd;
        plot(epsilon_0*10,tr,‘*‘);
    end
end
legend(‘epsilon = 0.1‘,‘epsilon = 0.3‘,‘epsilon = 0.5‘,‘epsilon = 0.7‘,‘epsilon = 0.9‘);