function R = sample_random_rotor( threshold )
    %% sample random vector but give fixed amount of rotation
    % strictly speaking when sampling thousands of them this may not be the
    % most elegant solution though
    
    % ax2om A.12 from Rowenhorst et al., MSMSE 23, 2015 
    N012 = [2.0*unifrnd(0.0,1.0)-1.0, 2.0*unifrnd(0.0,1.0)-1.0, 2.0*unifrnd(0.0,1.0)-1.0];
    N012 = N012 ./ (N012(1)^2+N012(2)^2+N012(3)^2)^0.5;
    omega = threshold / 180.0*pi;
    disp(['Threshold ' num2str(omega,'%.4f')]);
    
    c = cos(omega);
    s = sin(omega);
    n1 = N012(1);
    n2 = N012(2);
    n3 = N012(3);
    
    R = [c+(1-c)*n1^2, ...
        (1-c)*n1*n2+s*n3, ...
        (1-c)*n1*n3-s*n2; ...

        (1-c)*n1*n2-s*n3, ...
        c+(1-c)*n2^2, ...
        (1-c)*n2*n3+s*n1; ...

        (1-c)*n1*n3+s*n2, ...
        (1-c)*n2*n3-s*n1, ...
        c+(1-c)*n3^2 ];
    % no rotation, identity
end

