%function [Nxr] = phasedensity(M,r);

%% [nr] = phasedensity(M);
%%
%%  Examples  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  
%%  Input and output arguments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%  M   (matrix) N x R matrix of N points in r-dimensional space
%%  r   (scalar) radius of sphere-counting
%%
%%  Nxr (vector) N x 1 vector of counts for sphere of radius r
%%
%%
%%  Authors  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Erin Miller   (emiller@csss.santafe.edu)
%%  Nicky Grigg   (grigg@csss.santafe.edu)
%%  Aaron Clauset (clauset@csss.santafe.edu)
%%
%%  Santa Fe Institute, Complex Systems Summer School 2003
%%  date created : 06.25.03
%%  last modified: 06.26.03
%%  MATLAB file (not necessarily OCTAVE compatible)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% path(path,'/amd/poptart/root/research/aaron/m_files/');

%% check dimensionality of input matrix
[N, R] = size(M);

%% calculate euclidean distances between all points
for i=1:N
    for j=i+1:N
      D(i,j) = sqrt(sum(diff([M(i,:); M(j,:)]).^2));
    end
end
Nxr = sum((triu(D) + tril(transpose(D)))<r - eye(size(D)));
Nxr = Nxr./sum(Nxr);

