% matlab script to solve problem 3 in HW 7
% Prabir Barooah
clear all;

num = 50*conv([1 1],[1 2]);
den = conv([1 5 40],[1 0.03 0.06]);
G = tf(num,den);


%open loop rise time and peak overshoot
time = [0:0.001:30]';
y = step(G,time);

%it is seen from G that the Dc gain is
DC_gain_G  =280/2.4;
y_ss = DC_gain_G*1; %steady state value of the output for unit step input
S = stepinfo(y,time,y_ss)

%rise time ~ 4 sec, peak overshoot ~ 220
%clearly not good


%% lead controller D(s) k(s+3)/(s+10)

k = 3 %works
D = k*tf([1 3],[1 10]);

Gcl = feedback(D*G,1);

time_cl = [0:0.01:20]';
step_resp_cl = step(Gcl,time_cl);

figure;
plot(time_cl,step_resp_cl,'b');

%easy to see from Gcl that
DC_gain_Gcl  =2567/2648;
y_cl_ss = DC_gain_Gcl*1; %steady state value of the output for unit step input
S = stepinfo(step_resp_cl,time_cl,y_cl_ss);

rise_time= S.RiseTime
Mp = S.Overshoot/100


%% damping ratio and nat freq of each complex conjugate pole pair:
k = 3;
D = k*tf([1 3],[1 10]);
Gcl = feedback(D*G,1);

[z,p,k] = zpkdata(Gcl,'v');
complexpoles = p(find(imag(p)~=0));

%this is bad programming: may not always work:
p1 = complexpoles(1);
p2 = complexpoles(3);

wn1 = sqrt(p1*p1');
zeta1 = -(p1 + p1')/2/wn1;

wn2 = sqrt(p2*p2');
zeta2 = -(p2 + p2')/2/wn2;