% Demo showing the modes of a compound pendulum with equal masses,
% set an equal distance apart.  Graphics look satisfactory for n <= 6
% Based on Sturm System derivation given by Gantmacher and Krein, 
% "Oscillation Matrices", AMS, Section III.2

 n = 4;       % number of masses  (1 <= n <= 6 looks best)

 m = 1;       % mass of each body
 l = .5;      % length between masses
 g = 1;       % gravitational constant

% set up system matrices

 C = (m/2)*eye(n);
 A =  (g*m)/(2*l)*(diag([1:2:2*n-1])-diag(ones(n-1,1),1)-diag(ones(n-1,1),-1));

% compute modes U and eigenvalues diag(D)
 [U,D] = eig((2/m)*A);    % equivalent to inv(C)*A

% sort modes so the ones with lowest frequency appear first
 [junk,indx] = sort(diag(D));
 U = U(:,indx);
 D = D(indx,indx);

%%%%%% everything else is graphical simulation

% set up stationary graphical elements (top block, vertical line, title)
 figure(1), clf                
 text(.5,0,'hi hi hi hi') 
 set(gcf,'doublebuffer','on')  % keep figure from flashing in older MATLABs
 for k=1:n
    subplot(1,n,k)
    fill([-1 1 1 -1],[n*l n*l (n+.5)*l (n+.5)*l],.5*[1 1 1]); hold on
    plot([0 0], l*[-.5 n],'-','color',.7*[1 1 1])
    title(sprintf('p = %f', sqrt(D(k,k))))
    axis([-1 1 -.5*l (n+.5)*l])                                % fix the axes
 end
 tax = axes('position', [.46 .035 .1 .05])                    % print the time 
 ttxt = text(0,0,sprintf('t = %f',0)), axis off

% iterate through time slices
 fprintf('\n -- press <control-c> to stop simulation -- \n\n');
 t = 0;
 while 1                  % infinite loop
    for k=1:n
        subplot(1,n,k)
        q = .5*U(:,k)*sin(sqrt(D(k,k))*t+pi/2);                    % horizontal displacements
        if t>0, delete(plt1(k)), delete(plt2(k)), end              % delete old chain
        plt1(k) = plot([q;0],(0:n)*l,'k-','linewidth',1); hold on  % draw new one (black line)
        plt2(k) = plot(q, (0:n-1)*l,'b.','markersize',40);         % draw new one (blue masses)
        axis([-1 1 -.5*l (n+.5)*l])                                % fix the axes
    end
    axes(tax)                                                      % update time...
    delete(get(gca,'Children'))                                    % delete old time
    text(0,0,sprintf('t = %f',t)), axis off                        % print new time
    drawnow
    t = t+.1;
 end