Probability for random rotations of the orthogonal vectors frame Let $\{A\}$ be the frame of $3$ orthogonal unit vectors $i,j,k$. The frame is rotating freely in 3D space. Origin of the frame - point $(0,0,0)$ - is fixed and endpoints of vectors are always on  the surface of the unit sphere. The movement of the frame is  random so probability of setting the frame in particular localization is the same for all possible orientations.    
Question: 

What is probability that
  the distance of at least one of the endpoints of vectors
  $i,j,k$ from $Oxy$ plane is less than $d$ ($0<d<1$),
  for example $d$ can be set $=0.5$.
  Opposite event would be if all endpoints of vectors
  $i,j,k$ are distanced greater than $d$ (maybe this event is simpler for consideration).

Discussion: 
If we would like to calculate probability for only a single rotating vector ( for example $i$) the calculations could be based simply on describing areas  of the sphere surfaces  assigned to events $E_i$ when $d_i < d$ or $\bar{E_i}$ when $d_i \geq{d}$. But here we have $3$ vectors constrained in the way characteristic for columns of rotation matrix.
What would be the best strategy for solution of this problem ? What kind of representation of the frame should be taken (Euler angles or axis-angle representation of rotation matrix or maybe  quaternions  ?
It is rather certain that for some values of $d$ for example $d= 0.9$ probability would be equal to $1$, for $d$ near to $0$ for example $d= 0.1$ would be rather small.. between these two values function of probability for defined event would be,  can we assume ?, monotonic...
So we have two possible cases of events with every endpoint of one of the vectors $i,j,k$, but they are not independent ones but constrained in a geometrical way. How to use these geometrical constraints in the solution?
 A: I'm going to compute the "opposite" probability because it is easier.

Let $\Omega$ be the probability that all $i, j, k$ are at a distance $> d$ from the plane $Oxy$.
Let $x,y,z$ be the orthogonal unit vectors for the canonical frame.
Let $R \in SO(3)$ be the rotational matrix that bring $x,y,z$ to $i,j,k$.
For any unit vector $n$, let $B(n)$ be the set of unit vectors at an angle smaller than $\theta = \cos^{-1}(d)$ to $n$. i.e.
$$B(n) = \big\{ m \in S^2 : \angle(m,n) < \theta \big\} = \big\{ m \in S^2 : m\cdot n > d \big\}
$$
For any unit vector $n = Rm\in S^2$, the condition that $n$ is at a distance $> d$ to the plane $Oxy$ is equivalent to
$$\begin{align}
n \in B(z) \cup B(-z) 
\iff & n \in B(z) \lor n \in B(-z)\\
\iff & m \in B(R^{-1}z) \lor  m \in B(R^{-1}(-z)) \iff m \in B(R^{-1}z) \lor -m \in B(R^{-1}z)\\
\iff & R^{-1}z \in B(m) \cup B(-m)
\end{align}$$
So the condition that all $i = Rx, j = Ry, k = Rz$ are at a distance $> d$ to the plane $Oxy$ is equivalent to
$$R^{-1}z \in ( B(x)\cup B(-x) )\cap (B(y) \cup B(-y)) \cap (B(z)\cup B(-z))$$
Using symmetry and taking expectation over $R$, we find
$$\begin{align}\Omega 
&= \frac{1}{4\pi}
\verb/Area/(( B(x)\cup B(-x) )\cap (B(y) \cup B(-y)) \cap (B(z)\cup B(-z))\\
&= \frac{2}{\pi}\verb/Area/(B(x)\cap B(y) \cap B(z))
\end{align}$$
Let $T = B(x) \cap B(y) \cap B(z)$.


*

*For $d > \frac{1}{\sqrt{3}}$, $T = \emptyset$ and  $\Omega = 0$.

*For $d < \frac{1}{\sqrt{3}}$, $T$ forms a spherical triangle with vertices at
$$A = (\sqrt{1-2d^2}, d, d),\quad B = (d,\sqrt{1-2d^2},d)\quad\text{ and }\quad C = (d,d,\sqrt{1-2d^2})$$
The vertices are joined by three small circular arcs with 


*

*length $\ell = \sqrt{1-d^2}\left[\frac{\pi}{2} - 2\sin^{-1}\left(\frac{d}{\sqrt{1-d^2}}\right)\right]$.

*geodesic curvature $k_g = \frac{d}{\sqrt{1-d^2}}$.


Let $\varphi_A, \varphi_B, \varphi_C$ be the external angles at $A, B, C$ respectively. By symmetry, they are equal to each other.
Imagine we walk along $\partial T$ in positive orientation.
At vertex $A$, the tangent vector changes from
$\displaystyle\;\frac{1}{\sqrt{1-d^2}}(d,0,-\sqrt{1-2d^2})$ to
$\displaystyle\;\frac{1}{\sqrt{1-d^2}}(-d,\sqrt{1-2d^2},0)$. This implies
$$\varphi_A = \varphi_B = \varphi_C = \cos^{-1}\left(\frac{-d^2}{1-d^2}\right) = \pi - \cos^{-1}\left(\frac{d^2}{1-d^2}\right)$$
Apply Gauss-Bonnet formula, we find
$$\begin{align}
\Omega = \frac{2}{\pi}\verb/Area/(T) 
&=
\frac{2}{\pi}\left[2\pi - (\varphi_A + \varphi_B + \varphi_C + 3 k_g \ell)\right]\\
&= \frac{2}{\pi}\left[
3\cos^{-1}\left(\frac{d^2}{1-d^2}\right) + 
6d \sin^{-1}\left(\frac{d}{\sqrt{1-d^2}}\right)
-\frac{\pi}{2}(2+3d)
\right]
\end{align}
$$
