%  PLOT POLES OF CLOSED LOOP SYSTEM

%  FILENAME: PLOTPOLE.M
%  DATE: October 1994
%  VERSION: 02
%  COPYRIGHT 1993, 1994 M.Gretzinger & T.Marlin
%            McMaster University   

% This m-file is called from STABLE.M. It plots the poles of the closed
%       loop system specified by the process and tuning data. This method
%       is only applicable when no deadtime exists in the process.

if theta > 0
        clc
        disp('Dead time exists in process, plotting of poles not possible.')
        disp('')
        disp('Press ENTER to continue')
        pause
end

% Define coefficients of characteristic equation polynomial.
% (The characteristic equation differs depending on whether or not the
%  integral mode is included in the controller, to ensure there is no
%  division by zero.)
if ti == 0
    coeff_s5 = 0 ;
    coeff_s4 = tau1*tau2*tau3^2 ;
    coeff_s3 = 2*xi*tau1*tau2*tau3 + (tau1 + tau2)*tau3^2 ;
    coeff_s2 = tau1*tau2 + 2*xi*(tau1+ tau2)*tau3 + tau3^2 + kc*kp*taulead*td ;
    coeff_s1 = tau1 + tau2 + 2*xi*tau3 + kc*kp*(taulead + td) ;
    coeff_s0 = kc*kp + 1 ;
else
    coeff_s5 = tau1*tau2*tau3^2 ;
    coeff_s4 = 2*xi*tau1*tau2*tau3 + (tau1 + tau2)*tau3^2 ;
    coeff_s3 = tau1*tau2 + 2*xi*(tau1+ tau2)*tau3 + tau3^2 + kc*kp*taulead*td ;
    coeff_s2 = tau1 + tau2 + 2*xi*tau3 + kc*kp*(taulead + td) ;
    coeff_s1 = kc*kp*(taulead + ti)/ti + 1 ;
    coeff_s0 = kc*kp/ti ;
end

chareqn = [coeff_s5 coeff_s4 coeff_s3 coeff_s2 coeff_s1 coeff_s0 ];

% Calculate roots of characteristic equation (poles)
poles = roots(chareqn);

%  scale the plot as 2 times the magnitudes
xscale=2*max(abs(real(poles)));
yscale=2*max(abs(imag(poles)));

% Plot solutions on complex plane
clg ; axis('auto')
hold off
plot(real(poles),imag(poles),'*g')
xlabel('Real')
ylabel('Imaginary')
title('POLES OF CLOSED LOOP SYSTEM')

% Plot lines for the x and y axis in order to easily visualize location
% of poles

hold on
%                                  show real and imag axes
plot(xscale*[-1 1],[0 0],'r', ... 
      [0 0],yscale*[-1 1],'-r')
hold off
figure (1)
pause
close all

% Display the numerical values of the poles
clc
disp('The poles are:')
disp(poles)
disp('Press ENTER to continue')
pause

% Check whether poles are correct, sub roots into characteristic equation,
% making sure fcn = 0.
% This is commented out. To carry out the check, remove % from the following
% 11 lines of code.
%for i = 1:length(poles)
%    if ti == 0
%       fcn = 1 + (kp*(taulead*poles(i) + 1)*kc*(1 + td*poles(i)))...
%              /((tau1*poles(i)+1)*(tau2*poles(i)+1)*...
%              (tau3^2*poles(i)^2 + 2*xi*tau3*poles(i) + 1))
%    else
%       fcn = 1 + (kp*(taulead*poles(i) + 1)*kc*(1 + 1/(ti*poles(i))+ td*poles(i)))...
%              /((tau1*poles(i)+1)*(tau2*poles(i)+1)*...
%              (tau3^2*poles(i)^2 + 2*xi*tau3*poles(i) + 1))
%    end
%end