%-----------------------------------------------------------------------
% Example 3.9.4: Rotating Masses and Torsional Springs 
%-----------------------------------------------------------------------

% Initialize

   clc                  % clear screen 
   clear                % clear variables
   n = 4;               % order of system 
   J1 = 1.5;            % moment of inertia 1 
   J2 = 3.8;            % moment of inertia 2 
   k1 = 5.1;            % spring constant 1 
   k2 = 8.3;            % spring constant 2 
   A = zeros (n,n); 
   
% Compute A 

   fprintf ('Example 3.9.4: Rotating Masses and Torsional Springs\n');
   A(1,2) = 1;
   A(2,1) = -(k1 + k2)/J1;
   A(2,3) = k2/J1;
   A(3,4) = 1;
   A(4,1) = k2/J2;
   A(4,3) = -k2/J2;  
   show ('A',A)
   
% Find eigenvalues and oscillation periods 

   lambda = eig(A);
   show ('lambda',lambda)
   t = abs(trace(A) - sum(lambda));
   show ('Trace check',t) 
   show ('T1 (sec)',abs(2*pi/imag(lambda(1))))
   show ('T2 (sec)',abs(2*pi/imag(lambda(3))))
%-----------------------------------------------------------------------