%
% Written by:
% -- 
% John L. Weatherwax                2009-04-21
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

clc; close all; clear; 

M_q = -0.568;
M_alpha = -17.98; % negated from books number
L_alpha_over_V = 1.237;
M_dE = -0.175; % negated from books number
L_dE_over_V = 0.001;

w_n = sqrt( -M_alpha )

xi = ( 0.5/(w_n) ) * ( L_alpha_over_V - M_q )

% the poles of the original transfer function (roots when k=0)
s1 = -xi*w_n - w_n * sqrt( xi^2 - 1 )
s2 = -xi*w_n + w_n * sqrt( xi^2 - 1 )


Nk     = 1000;
kArray = linspace(-100,+100,Nk);
S1     = zeros(1,Nk);
S2     = zeros(1,Nk);
for ki = 1:Nk,
  % form the coefficients of the polynomial we want the roots of:
  k = kArray(ki);
  A = 1;
  B = 2 * xi * w_n - k * L_dE_over_V;
  C = w_n^2 + k * L_dE_over_V * ( M_q - M_dE / L_dE_over_V );
  coef = [ A, B, C ];
  rt = flipud( roots( coef ) ); % order these roots so that they are the same as 
  
  S1(ki) = rt(1);
  S2(ki) = rt(2);
end

indsN = find( kArray<0 );
indsP = find( kArray>=0 );

figure; 
phs1 = plot( real(S1(indsN)), imag(S1(indsN)), '.b' ); hold on;
phs1 = plot( real(S1(indsP)), imag(S1(indsP)), 'xb' ); 
plot( real(s1), imag(s1), '.k', 'markersize', 30 );

phs2 = plot( real(S2(indsN)), imag(S2(indsN)), '.r' ); 
phs2 = plot( real(S2(indsP)), imag(S2(indsP)), 'xr' ); 
plot( real(s2), imag(s2), '.k', 'markersize', 30 );

legend( [phs1,phs2], {'root 1(k)','root 2(k)'} );
xlabel('real part'); ylabel('imag part');
grid on;

%saveas( gcf, '../../WriteUp/Graphics/Chapter2/sect_6_prob_1_part_b_root_locus.eps', 'epsc' );