function angles = getangles(M)
%GETANGLES Returns a matrix of angles between the columns of M.
%	This function treats the matrix M as a series of columns.
%	All possible angles between the columns of M are computed.
%	The value of "999" is treated as a blank in the column, and
%	that position is neglected in computing an angle with another
%	column (999 can be changed to an arbitrary value).  The function 
%	returns a square matrix where each position represents an angle
%	between the appropriate columns.	For example, after typing:
%		angles = getangles(M)
%	The value of angles(2,6) will be the angle between the second 
%	and sixth columns of M.
%
% Written by Finny Kuruvilla	 v1.6.0	Last updated 2-16-2002

blank_marker=999;				% Change to another value if desired.
[rows,cols]=size(M);
angles=zeros(cols,cols);

for k1=1:(cols-1)
   for k2=(k1+1):cols
      common=(M(:,k1)~=blank_marker).*(M(:,k2)~=blank_marker);
      col1=M(:,k1).*common;
      col2=M(:,k2).*common;
      len1=sum(col1.*col1);
      len2=sum(col2.*col2);
      dotproduct=sum(col1.*col2);
      theta=180*(acos(dotproduct/sqrt(len1*len2)))/pi;
      angles(k1,k2)=theta;
      % If "in progress" update statements are not desired, comment out
      %	the following line.
      disp(sprintf('Angle between column %d and column %d is %d',k1,k2,theta));
   end
end

% Since an angle computation is commutative, instead of recomputing, simply add 
%	the transpose of the upper triangular portion to get the full matrix.
angles=angles+angles';