% define the tf of a LTI system (in this case, II order system with complex poles) 
% and then evaluate the step response

clear all
close all
clc

% the LTI system is described by the following W(s):
% W(s)=N(s)/D(s)=4/(s^2+s+2) 

% define the N(s) coefficients
num1=4;
% define the D(s) coefficients
den1=[1 1 2 ]; 
% define W by using tf function
sys=tf (num1,den1);


% compute zita and omega_n: D(s)=s^2+a1*s+a0, 
% in terms of zita and omega_n, D(s)= s^2+2*zita*omega_n*s+omega_n^2
% then omega_n=sqrt(a0) and zita=a1/2/sqrt(a0)
omega_n=sqrt(2);
zita=1/2/omega_n;

% compute the static gain, i.e. k=W(0)
k=dcgain(sys);
% compute the steady-state value of the step response (y_ss)
U0=1; %amplitude of the step signal
y_ss=k*U0;

% compute zita and omega_n by using damp function (when zita <1)
[Wn,Z] = damp(sys);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the following formulas are valid when zita is small (zitta<<1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% evaluate the settling time, ts5%
ts5=3/Z(1)/Wn(1);

% evaluate the peak time, tp
tp=pi/Wn(1)/sqrt(1-Z(1)^2);

% evaluate the overshoot, 
over_s=exp(-pi*Z(1)/sqrt(1-Z(1)^2));

% evaluate the peak value
ymax=yss*(1+over_s);

% evaluate the oscillation period, T
T=2*pi/Wn(1)/sqrt(1-Z(1)^2);

% evaluate rise time
tr=1/omega_n;



% another way to define W
s=tf('s');
W=4/(s^2+s+2);


% evaluate the step response
figure (1)
step(W)
hold on
grid on

figure(2)
step(sys)
hold on
grid on

% return the parameters characterizing the step response
% by default the settling is computed by assuming an error (yss-y_t) within
% 2%
stepinfo(W)