Rotating two vector randomly keeping the relative orientation between them unchanged Let's consider we have two vectors $A(x_1, y_1, z_1)$ and $B(x_2, y_2, z_2)$. 
Now want to rotate these two vectors in $3D$ space (such that the relative orientation between them is always same).
How can I do that?
PS. I have tried to rotate the two vectors independently by angles $\phi$, $\theta$, and $\psi$. But I did not get the uniform spherical distribution. But I was getting a higher intensity near the pole of the sphere.
 A: Use the axis-angle representation for the rotation.
To generate the axis, I recommend generating random vectors within $x = -1 .. +1, y = -1 .. +1, z = -1 .. +1$, accepting the first one whose length squared is at least $1/4$ but not more than $1$, and normalizing it to unit length:
Function RandomUnitVector():
    Do
        x = random number from -1 to +1, inclusive
        y = random number from -1 to +1, inclusive
        z = random number from -1 to +1, inclusive
        rr = x*x + y*y + z*z
    While (rr < 0.25 || r > 1)

    r = sqrt(rr)
    x = x / r
    y = y / r
    z = z / r

    Return (x, y, z)
End Function

Also generate an uniformly random angle $\varphi$ between $-180°$ and $+180°$ ($-\pi$ to $\pi$), including one excluding one of the limits.
If $(x, y, z)$ is the unit axis vector, $c = \cos \varphi$, $v = 1 - c$, and $s = \sin \varphi$, the rotation matrix $\mathbf{R}$ is
$$\mathbf{R} = \left [ \begin{matrix}
c + x^2 v    &  x y v - z s  &  x z v + y s \\
x y v + z s  &  c + y^2 v    &  y z v - x s \\
x z v - y s  &  y z v + x s  &  c + z^2 v
\end{matrix} \right ]$$
Because the angle of rotation is uniformly distributed (across all unique angles), and the axis is uniformly distributed within the surface of an unit sphere, this rotation matrix produces uniformly distributed results.
