I am creating a rotation matrix capable of converting the reference system $A=(x_a, y_a, z_a)$ into a second reference system $B=(x_b, y_b, z_b)$ where the $x_b$ and $y_b$ axes lie in a plane with normal $n$ and the $z_b$ axis is parallel to the normal $n$.
I am using the following formula in order to convert a $P_A=$ point to a $P_B$ point:
$$P_B= P_A \pmatrix{ {n_y\over\sqrt{n_x^2+n_y^2}}& {-n_x\over\sqrt{n_x^2+n_y^2}}& 0\\ {n_xn_z\over\sqrt{n_x^2+n_y^2}}& {n_yn_z\over\sqrt{n_x^2+n_y^2}}& -\sqrt{n_x^2+n_y^2}\\ n_x& n_y& n_z }= \pmatrix{ p_x*{n_y\over\sqrt{n_x^2+n_y^2}} + p_y*{-n_x\over\sqrt{n_x^2+n_y^2}} + p_z * 0\\ p_x*{n_xn_z\over\sqrt{n_x^2+n_y^2}} + p_y*{n_yn_z\over\sqrt{n_x^2+n_y^2}} + p_z * -\sqrt{n_x^2+n_y^2}\\ p_x*n_x + p_y*n_y + p_z * n_z\\ }$$
however, trying this formula with a javascript code I get the following result (in red the transformed $x$ axis, in green the transformed $y$ axis, in blue the transformed $z$ axis and in yellow the normal $n$ to the gray plane below):
Does the rotation matrix seem correct to you?
The original coordinate frame is specified as follows
$ p = R q $
where $p$ is the world coordinates and $q$ is the coordinates with respect to the axes of your initial reference frame. These axes have unit vectors that are the three columns of the $3 \times 3$ matrix $R$, where the third column is your vector $n$. From here
$q = R^T p $
If you rotate the axes about $n$, then the new axes will be given by
$ R' = R_z(\theta) $
where
$R_z(\theta) = \begin{bmatrix} \cos(\theta) && - \sin(\theta) && 0 \\ \sin(\theta) && \cos(\theta) && 0 \\ 0 && 0 && 1 \end{bmatrix} $
Hence if the new coordinates of $q$ are $q'$ then
$q = R' q'$
Thus $q' = R'^T q = R'^T R^T p $
And finally, $ p = R R' q' $
So the new axes are specified by the matrix $R R'$