function M = randrotmat(sz)
	% R = randrotmat(sz)
	%	generate random rotation matrices R
	%
	% Input:
	%	sz	array	number of matrices to create
	% Output:
	%	R	3 x 3 x sz(1) x sz(2) ...
	%				uniform random rotation matrix
	%
	% References:
	%	https://en.wikipedia.org/wiki/Rotation_matrix#Uniform_random_rotation_matrices
	%	https://archive.org/details/isbn_9780124096738/page/118
	
	x	= rand([3 1 sz]);
	R	= bsxfun(@times, [1 0 0; 0 1 0; 0 0 0], cos(2*pi*x(1,1,:))) ...
		+ bsxfun(@times, [0 1 0;-1 0 0; 0 0 0], sin(2*pi*x(1,1,:))) ...
		+ [0 0 0; 0 0 0; 0 0 1];
	v	= cat(1, cos(2*pi*x(2,1,:)) .* sqrt(x(3,:,:)), ...
				 sin(2*pi*x(2,1,:)) .* sqrt(x(3,:,:)), ...
				 sqrt(1 - x(3,:,:)));
	H	= eye(3) - 2*mmx('mult', v, v, 'nt');	
	M	= -mmx('mult', H, R);
	M	= reshape(M, [3 3 sz]);
end