I would like to compose/decompose an $n$-dimensional orthogonal rotation matrix (restricting to simple planar rotations, which rotates in the specified plane of rotation, and fixes in the plane orthogonal to the plane of rotation) from/to the two basis vectors of the plane of rotation, and the angle of rotation.

The most common method in 3 dimensions is composing/decomposing the rotation matrix from/to an axis and an angle, but this doesn't work in higher dimensions.

For example in $\mathbb{R}^3$ the well-known rotation matrix

$R_{xy}=\begin{bmatrix}\cos\theta & -\sin\theta & 0 \\ \sin\theta & \cos\theta & 0 \\ 0 & 0 & 1 \end{bmatrix}$

it's obvious that the plane of rotation is the $xy$-plane spanned by basis vectors $b_0 = (1,0,0)$ and $b_1=(0,1,0)$ and the angle of rotation is $\theta$. However composing/decomposing it mathematically is rather challenging.

What is the solution for a general (restricted to a single, but arbitrary plane) rotation matrix in $\mathbb{R}^n$?