function [n,u,v]=RandAxis(N)
% RandAxis Generate 3D Axes in random direction n.

% revision history:
% 04/13/98 Mark D. Shattuck <mds> RandAxis.m  
% 04/03/02 mds changed to cross for clarity
% 09/21/14 mds noticed u,v are not uniform (fixed)
%              Switch t->u s->v

n=randn(3,N);       % generate random direction
r=sqrt(sum(n.^2));
n=n./repmat(r,3,1); % normalize

g=2*pi*rand(1,N);  % generate random angle for perpendicular unit vector

sg=sin(g);
cg=cos(g);
st=sqrt(1-n(3,:).^2);

u=[(sg.*n(2,:)-cg.*n(3,:).*n(1,:))./st;...
  -(sg.*n(1,:)+cg.*n(3,:).*n(2,:))./st;...
           (cg-cg.*n(3,:).*n(3,:))./st;];
         
v=cross(n,u);


% old way u is denser on the poles and v is denser in the equator
% s=sin(th);
% c=cos(th);
% r=sqrt(1-(n(1,:).*c-n(2,:).*s).^2);
% s=s./r;
% c=c./r;
% 
% t=[n(3,:).*s;n(3,:).*c;-n(1,:).*s-n(2,:).*c]; % random norm perp
% s=cross(n,t);

% % t2=cross(n,t1) perp to both
% t2=[t1(2,:).*n(3,:)-t1(3,:).*n(2,:);...
%     t1(3,:).*n(1,:)-t1(1,:).*n(3,:);...
%     t1(1,:).*n(2,:)-t1(2,:).*n(1,:)];
