function ThreeMassSystem(m1,m2,m3,k12,k23)
% function ThreeMassSystem(m1,m2,m3,k12,k23)

if nargin==0
   m1 = 1; m2 = 0.1; m3 = 0.1;
   k12 = 0.1; k23 = 0.1;
end

M  = [m1 0  0
      0  m2 0
      0  0  m3];
K  = [ k12 -k12      0
      -k12  k12+k23 -k23
       0   -k23      k23];
[V,D] = eig(K/M);
PHI = M\V;
W   = sqrt(abs(diag(D)));
[W,ind] = sort(W);

x_offset = [2 5 8];

t = 0;
A = [0 1 0 1
     1 0 1 0
     0 1 0 1
     1 0 1 0];
while 1
   for i=1:3
      w = W(i);
      phi = PHI(:,ind(i));
      phi = phi/max(abs(phi));
      t = t+0.1;
      x = phi*cos(w*t);

      C1 = [x(1)+x_offset(1)-0.5 0
            x(1)+x_offset(1)-0.5 1
            x(1)+x_offset(1)+0.5 1
            x(1)+x_offset(1)+0.5 0];
      C2 = [x(2)+x_offset(2)-0.5 0
            x(2)+x_offset(2)-0.5 1
            x(2)+x_offset(2)+0.5 1
            x(2)+x_offset(2)+0.5 0];
      C3 = [x(3)+x_offset(3)-0.5 0
            x(3)+x_offset(3)-0.5 1
            x(3)+x_offset(3)+0.5 1
            x(3)+x_offset(3)+0.5 0];
      eval(['subplot(31' num2str(i) ')'])
      gplot(A,C1)
      hold on
      gplot(A,C2)
      gplot(A,C3)
      plot([x(1)+x_offset(1)+0.5 x(2)+x_offset(2)-0.5],[0.5 0.5])
      plot([x(2)+x_offset(2)+0.5 x(3)+x_offset(3)-0.5],[0.5 0.5])
      axis equal
      axis([0 10 -1 2])
      
      if phi(1)~=phi(2),alpha12 = -phi(2)/(phi(1)-phi(2)); else alpha12=0;end
      if phi(2)~=phi(3),alpha23 = -phi(3)/(phi(2)-phi(3)); else alpha23=0;end
      if alpha12>0 & alpha12<=1
         plot((alpha12*(x(1)+x_offset(1))+(1-alpha12)*(x(2)+x_offset(2)))*[1 1],[0 1],'r')
      end
      if alpha23>0 & alpha23<=1
         plot((alpha23*(x(2)+x_offset(2))+(1-alpha23)*(x(3)+x_offset(3)))*[1 1],[0 1],'r')
      end
      
      hold off
   end
   pause(0.01)

end