Inverse of a skew-symmetric matrix For $a, x, y, z \in \mathbb R$, let $$M= \left(
  \begin{array}{cccc}
    \cos(a)  & \sin(a) \, x & \sin(a)\,  y & \sin(a) \, z \\
    -\sin(a) \, x & \cos(a) & \sin(a) \,z & -\sin(a)\, y \\
   -\sin(a) \, y & -\sin(a) \, z & \cos(a) & \sin(a) \, x \\
  -\sin(a) \, z & \sin(a) \, y & -\sin(a)\, x & \cos(a)
  \end{array}
\right).$$ 
Without doing much calculation, why the matrix $M-I_4$ is invertible and why its inverse is given $$(M-I_4)^{-1}=\,\frac{-1}{2} I_{4} - \frac{\cot(\sqrt{x^2+y^2+z^2} )}{2\sqrt{x^2+y^2+z^2}} A,$$
where $A$ is the skew-symmetric matrix given by 
$$A=\left(
  \begin{array}{cccc}
    0 & x & y &  z \\
    -x & 0 & z & -y \\
    -y & -z & 0 & x \\
    -z & y & -x & 0
  \end{array}
\right).$$
Thank you in advance
 A: Here's an approach which is relatively light on computation.
As in the comments, we observe that $A^2 = -(x^2 + y^2 + z^2)I$.  It follows that any rational function of $A$ can be written in the form $p I + q A$ for some $p,q$.  In other words, there necessarily exist numbers $p,q$ such that
$$
(M-I)^{-1} = ([\cos(a) - 1]I + \sin(a)A)^{-1} = p I + q A
$$
Now, expand the product $(M-I)(pI + qA),$ and solve for the $p$ and $q$ that cause this product to equal $I$.  In particular, we have
$$
(M - I)(pI + qA) = ([\cos(a) - 1]I + \sin(a)A)(pI + qA)\\
= p(\cos(a) - 1)I + (q[\cos(a) - 1] + p\sin(a))A + \sin(a)qA^2\\
= \left[p(\cos(a) - 1) - \sin(a)(x^2 + y^2 + z^2)q \right]I +  \left[q(\cos(a) - 1) + p\sin(a)\right]A.
$$
Thus, it suffices to find $p,q$ that solve the system of equations
$$
\begin{array}{ccccccc}
(\cos(a) - 1) &p &- &\sin(a)(x^2 + y^2 + z^2)&q &= &1\\
\sin(a) &p &+ &(\cos(a) - 1)&q &= &0.
\end{array}
$$
Solving this system yields
$$
p = \frac{1 - \cos(a)}{D}, \quad q = \frac{\sin(a)}{D}
$$
where $D$ is the determinant of the coefficient matrix of this system, that is
$$
D = \sin^2(a)(x^2 + y^2 + z^2) + (\cos(a) - 1)^2.
$$
I see no general way of simplifying the resulting expression.
A: Here is a proof of the first part (inversibility of $M-I_4$)
and a computation of the inverse (though I do not get a result in the form of yours).
One obtains an efficient simplification by using half-angle formulas (https://en.wikipedia.org/wiki/Tangent_half-angle_substitution) :
$$\cos(a)=\dfrac{1-t^2}{1+t^2} \ \  \text{and} \ \ \sin(a)=\dfrac{2t}{1+t^2}$$ 
Let $N=M-I_4$.
All entries in matrix $N$ have a common factor $\dfrac{2t}{1+t^2}$. Therefore we can set
$$N=\dfrac{2t}{1+t^2}P \ \ \text{with} \ \  
P:=\begin{bmatrix}-t&x&y&z\\-x&-t&z&-y\\-y&-z&-t&x\\-z&y&-x&-t\end{bmatrix}$$
Therefore, an equivalent issue is to show the invertibility of $P$.
Please note that $P=A-tI_4$ with matrix $A$ as you have defined it.
It turns out that the determinant of $P$ is very compact :
$$\det(P)=(t^2 + x^2 + y^2 + z^2)^2\tag{1}$$
which is non-zero, excepted... exceptional cases $x=y=z=0$ and $t=0$ (the last case corresponding to angles $a=k\pi, k \in \mathbb{Z}$). Set apart these cases, $P$ is always invertible.
Now, what is the inverse of $M-I$ ? The inverse of $P$, considered as a quaternion (see remark 1 below) is
$$P^{-1}=\dfrac{1}{n}\begin{bmatrix}-t&-x&-y&-z\\x&-t&-z&y\\y&z&-t&-x\\z&-y&x&-t\end{bmatrix}=\dfrac{1}{n}(-A-tI_4)\tag{2}$$
$$\text{with }n:=t^2+x^2+y^2+z^2\tag{3}$$ 
(one can easily check (2) by multiplying $P$ and $P^{-1}$).
Therefore, as $M-I_4=\dfrac{2t}{1+t^2}P=\sin(a)P$, using (2) and (3), we get
$$(M-I_4)^{-1}=\dfrac{1}{\sin(a)}P^{-1}=\dfrac{1}{n\sin(a)}(A-tI_4)\tag{4}$$
which is different from your formula (as it has been remarked, this formula needs to depend on $a$).
Remarks : 
1) It is not surprizing to get formula (1) : it is to be related to the $4 \times 4$ matrix representation of quaternions (see paragraph "Matrix representation" in https://en.wikipedia.org/wiki/Quaternion) and their so-called "norm". By the way, I wouldn't be that surprised that your issue comes from quaternionic representations...
2) Using (1) with $t=0$, we get 
$$\det(A)=(x^2 + y^2 + z^2)^2\tag{2}$$
3) Formula (1) has to be connected to the so-called Pfaffian. If we have had to compute the determinant of a skew-symmetric matrix of even order, it is good to know that there exists a formula as a square of a polynomial expression in its entries called its associated Pfaffian. See https://en.wikipedia.org/wiki/Pfaffian where you will find the following formula for order $4$ :
$$\det\begin{bmatrix}0&a&b&c\\-a&0&d&e\\-b&-d&0&f\\-c&-e&-f&0\end{bmatrix}=(af-be+dc)^2$$
A: $A$ has the form $A = \begin{pmatrix}J&L\\-L&J\end{pmatrix}$, where $JL+LJ = 0$, $J^2 = -x^2I$, and $L^2 = (y^2+z^2)I$. Therefore,
$$
A^2 = \begin{pmatrix}J&L\\-L&J\end{pmatrix}\begin{pmatrix}J&L\\-L&J\end{pmatrix} = \begin{pmatrix}J^2-L^2&0\\0&J^2-L^2\end{pmatrix} = -r^2I_4,
$$
where $r := \sqrt{x^2+y^2+z^2}$. Hence, for $z\in\mathbb C$ we have $(A-z)(A+z) = A^2-z^2I = -(r^2+z^2)I$ and thus
$$
(A-z)^{-1} = -\frac{A+z}{r^2+z^2}.
$$
Now, as I wrote in the comments, we have $M-I = (\cos(a)-1)I + \sin(a)A$. So, for $a = 2k\pi$ the matrix $M-I$ is singular, whereas for $a=(2k+1)\pi$ we have $M-I = -2I$. If $a$ is not one of these values,
$$
(M-I)^{-1} = \left[\sin(a)\left(A - \frac{1-\cos(a)}{\sin(a)}\right)\right]^{-1} = \frac 1{\sin(a)}\cdot\frac{A+\alpha}{r^2+\alpha^2},
$$
where $\alpha = \frac{1-\cos(a)}{\sin(a)}$.
