% poleplacement.m
%
% Script models a standard third-order SISO system
%   x'''(t) + a_2 x"(t) + a_1 x'(t) + a_0 x(t) = u(t),
% shows state response with u(t)=0 from initial cdt
%   x(0)=1, x'(0)=0, x"(0)=0
% Various eigenvalue options are shown all together.

% Philip D Loewen, 2012-09-28

% Eigenvalues at  s+iw, s-iw, r
% Make comparison plot for a list of up to 5 choices

Slist = [-0.10;-0.11;-0.12;-0.10; -0.11];
Rlist = [-0.10;-0.20;-0.21;-0.22; -0.21];
Wlist = [ 1.00; 1.00; 2.00; 0.50;  0.21];
Clist = ['r';'g';'b';'m';'c'];
Hlist = zeros(5,1);

% Here my favourite three choices supercede definitions just given

Slist = [-0.10;-0.10;-0.10];
Rlist = [-0.10;-0.20;-0.20];
Wlist = [ 1.00; 1.00; 0.20];
Clist = ['r';'g';'b';'m';'c'];
Mlist = ['x';'o';'+';'*';'-'];
Hlist = zeros(5,1);

lw = 2; % LineWidth for trajectory plots

% BUILD FIGURE 1 -- STATE TRAJECTORIES
figure('Name','Responses');

for k = 1:length(Slist),
  s = Slist(k);
  r = Rlist(k);
  w = Wlist(k);
  colour = Clist(k);
  marker = Mlist(k);
  
  A = [0,1,0; 0,0,1; (w^2+s^2)*r -(w^2+s^2+2*r*s) (r+2*s)];
  B = [0;0;1];
  F = [0,0,0];

  % Calculate and print some confirming info
  [V,D] = eig(A);
  DD=diag(D);
  disp([colour,marker,': ','lambda = [',cfmt(DD(1)),', ',cfmt(DD(2)),', ',cfmt(DD(3)),'].']);

  dx = @(t,x) (A + B*F)*x;  % Closed-loop system
  x0 = [1;0;0];

  % Solve dynamics on interval [0,60] with initial value x0
  [TTT,XXX ] = ode45(dx,[0,60],x0);
  Hlist(k) = plot(TTT,XXX(:,1),colour,'LineWidth',lw);
  hold on;
  end;

% Add axes to plot
plot([0,0],10*[-1,100],'k:');
hold on;
plot(100*[-1,1],[0,0],'k:');

% Set plot window
set(gca,'YLim',1.1*[-0.5,1.1]);
plot(100*[-1,0],[1,1],'k','LineWidth',lw);
set(gca,'XLim',[-10,60]);
xlabel('t');
ylabel('x(t)');


% BUILD FIGURE 2 -- POLE LOCATIONS
figure('Name','Complex Plane');
plot([0,0],10*[-1,1],'k:');       % The imaginary axis
hold on;
plot([-10,1],[0,0],'k:');         % The real axis
for k=1:length(Slist),
  s = Slist(k);
  r = Rlist(k);
  w = Wlist(k);
  colour = Clist(k);
  cm = [colour,Mlist(k)];

  plot(s, w,cm,'MarkerSize',14,'LineWidth',lw);
  plot(s,-w,cm,'MarkerSize',14,'LineWidth',lw);
  plot(r, 0,cm,'MarkerSize',14,'LineWidth',lw);
  end;

set(gca,'YLim',1.1*[-1,1],'XLim',[-1,0.1]);
axis equal;
set(gca,'YLim',1.1*[-1,1],'XLim',[-1,0.1]);

xlabel('Re(s)'); ylabel('Im(s)');