%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  [V, D] = spring_plot( <masses> , <spring const> , <cos coeff> , <sin coeff> )
%
%    Plots the position of oscillating masses for line of masses  
%    with the given cos and sin coefficients
%   _______
%      |
%      |   c_1
%     m_1
%      |   c_2
%     m_2 
%      |   c_3
%     m_3
%      |   c_4
%   ___|___ 
%    \ \ \ \
%
%    Note: Modes of oscillation are sorted by eigenvalue
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [Vout,Dout] = spring_plot(m, c, a, b)

  if ( nargin < 3 )
    fprintf('\n    Use:  [V, D] = spring_plot( <masses>, <spring const>, <cos coeff>, <sin coeff> )\n\n')
    return
  end

  if ( min(c) < 0 )
    error('Spring constants cannot be negative!')
  end

  if ( min(m) <= 0 )
    error('Masses cannot be zero or negative!')
  end

  if ( length(c) < length(m) )
    error('Not enough spring constants!')
  end

  if ( length(c) > length(m)+1 )
    error('Too many spring constants!')
  end

  if ( length(c) == length(m) )
    c = [c 0];
  end

  if (exist('b') == 0)
    b = 0;
  end

  if ( length(a) > length(b) )
    b = [b zeros(1,length(a)-length(b))];
  elseif ( length(b) > length(a) )
    a = [a zeros(1,length(b)-length(a))];
  end

  if ( length(a) > length(m) )
    error('Too many coefficients!')
  end




  K = diag( c(1:end-1) + c(2:end) ) - ...
      diag( c(2:end-1),  1 )        - ...
      diag( c(2:end-1), -1 );

  Minv   = diag( 1./m );
  [V, D] = eigs( Minv * K , length(m) ); 


  for ( k = 1:length(m) )
    if ( V(1,k) < 0 )
      V(:,k) = -V(:,k);
    end
  end


  freq = sqrt( diag(D) );

  if ( D(1,1) > 0.0001 )
    U = @(i,t)  a(1) * V(i,1) * cos( freq(1) * t ) + ...
                b(1) * V(i,1) * sin( freq(1) * t );
  else
    U = @(i,t)  a(1) * V(i,1)                      + ...
                b(1) * V(i,1) * t;
  end

  for ( k = 2:length(a) )
    U = @(i,t)  U(i,t) +                             ...
                a(k) * V(i,k) * cos( freq(k) * t ) + ...
                b(k) * V(i,k) * sin( freq(k) * t );
  end

  center = [0,  min( U( 1, 0:.01:10*pi) - U(2, 0:.01:10*pi) ) - .5];

  for( k = 3:length(m) )
    center(k) = center(end) + ...
                min( U(k-1,0:.01:10*pi) - U(k, 0:.01:10*pi) ) - .5;
  end

  figure
  hold on;

  for( mass = 1:length(m) )
    fplot( @(t) U(mass,t) + center(mass),  [0, 10*pi] );
  end

  if (nargout > 0)
    Vout = V;
  end

  if (nargout == 2)
    Dout = D;
  end

end
