Some Matlab code for generating random normal vectors uniformly distributed in angle

randomUnitVectors.m

%% First a method using quaternions
qtn = randrot(1000,1);

pts = rotatepoint(qtn, [1 0 0]);

figure()
scatter3(pts(:,1), pts(:,2), pts(:,3))
axis equal


%% Second method (works only in the half sphere)
pts = zeros(3,500);

for ii = 1:length(pts)
[Q,~] = qr(randn(3));
pts(:,ii) = Q*[1 0 0]';
end

figure()
scatter3(pts(1,:), pts(2,:), pts(3,:))
axis equal