genuniD(N,D,RandSeed)

function U = genuniD(N,D,RandSeed)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% U = genUniD(N,D,RandSeed) RandSeed=-1 for random
%% generate N D dimensional vectors of unit magitude;
%%
%% Oct. 29, 2002 Released with NESim
%% Feb. 25, 2001 Consolidation
%% Jun. 18, 2003 Implements PC algorithm (implemented by John Conklin)
%%
%% The PC algorith is described in: 
%% "An Algorithm for Uniform Random Sampling of 
%% Points In and On a Hypersphere."
%% Authors: Gerald Guralnik, Charles Zemach, Tony Warnock
%% Information Processing Letters, Volume 21, 
%% Number 1, 10 July 1985, p. 17-21
%%
%% Copyright 2002 C. H. Anderson ([e-mail]) and
%% C. Eliasmith ([e-mail])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(D>10)
	%error('genUniD cannot work higher than 10 dimensions\n');
end
if(N>5000)
	error('genUniD does not work with more than 5000 neurons\n');
end
if(nargin==3)
	if(RandSeed>0)
		rand('state',RandSeed);
	end
end
if(D==1)
   U=ones(1,N);
   U((floor((N+1)/2)+1):N)=-1;
   U=U';
else 
   % Case D = 2d = even
   if mod(D,2)==0
       d=D/2;
       theta=2*pi*rand(N,d);
       u=rand(N,d-1);
       row=ones(N,1)*(1:d-1);
       rhobar=u.^(1./(D-2*row));
       rho=ones(N,d);
       rho(:,1:d-1)=sqrt(1-rhobar.^2);
       x=cos(theta).*rho;
       y=sin(theta).*rho;
       for i=1:d-1
           x(:,i+1)=x(:,i+1).*prod(rhobar(:,1:i),2);
           y(:,i+1)=y(:,i+1).*prod(rhobar(:,1:i),2);
       end
   else
   % Case D = 2d+1 = odd
       d=floor(D/2);
       theta=2*pi*rand(N,d);
       u=rand(N,d);
       row=ones(N,1)*(1:d);
       rhobar=u.^(1./(D-2*row));
       rho=ones(N,d+1);
       rho(:,1:d)=sqrt(1-rhobar.^2);
       rho(:,d+1)=2*floor(2*rand(N,1))-1;
       x=rho;
       x(:,1:d)=cos(theta).*x(:,1:d);
       y=sin(theta).*rho(:,1:d);
       for i=1:d-1
           x(:,i+1)=x(:,i+1).*prod(rhobar(:,1:i),2);
           y(:,i+1)=y(:,i+1).*prod(rhobar(:,1:i),2);
       end
       x(:,d+1)=x(:,d+1).*prod(rhobar(:,:),2);
   end
   U=cat(2,x,y);
end