Let the point $p$ be given in cartesian coordinates by $(\sin\theta\cos\phi, \sin\theta\sin\phi, \cos\phi)$, so that $\theta$ and $\phi$ are the usual spherical coordinates on $S^2$. Then the matrix that rotates $n=(1,0,0)$ to $p$, and also preserves the plane passing through both of them and the origin, is given by
$M=\left(\begin{array}{ccc}\sin^2\phi + \cos\theta\cos^2\phi & (-1+\cos\theta)\sin\phi\cos\phi & \sin\theta\cos\phi\\
(-1+\cos\theta)\sin\phi\cos\phi & \cos^2\phi+\cos\theta\sin^2\phi & \sin\theta\sin\phi\\
-\sin\theta\cos\phi & -\sin\theta\sin\phi & \cos\theta
\end{array}\right).$
It's straightforward but tedious to check that this is indeed an orthogonal matrix with determinant $1$ that sends $n$ to $p$ and fixes $u=n\times p=(-\sin\theta\sin\phi,\sin\theta\cos\phi,0)$, and hence the desired transformation.
I produced this matrix by considering that, for each fixed $\phi$, what you desire is a one-parameter subgroup of $SO(3)$: the group of rotations that fix $u$. It therefore is the family of exponentials of an infinitesimal transformation, i.e. an element of the Lie algebra $\frak{so}\mathrm{(3)}$. For example, a tiny rotation tilting the $z$-axis toward the $x$-axis differs from the identity by
$\left(\begin{array}{ccc}
& & \epsilon\\
& & \\
-\epsilon & &
\end{array}\right)$, and
$\exp\left(\begin{array}{ccc}
& & \theta\\
& & \\
-\theta & &
\end{array}\right) = \left(\begin{array}{ccc}
\cos\theta & & \sin\theta\\
& 1 & \\
-\sin\theta & & \cos\theta
\end{array}\right).$
Similarly, the Lie algebra element corresponding to an infinitesimal rotation of the $z$-axis in the direction of the $y$-axis is given by $\left(\begin{array}{ccc}
& & \\
& & \epsilon\\
& -\epsilon&
\end{array}\right)$, and
$\exp\left(\begin{array}{ccc}
& & \\
& & \theta\\
& -\theta &
\end{array}\right) = \left(\begin{array}{ccc}
1 & & \\
& \cos\theta & \sin\theta\\
& -\sin\theta & \cos\theta
\end{array}\right).$
We want the Lie algebra element which has components $\cos\phi$ toward the $x$-axis and $\sin\phi$ toward the $y$-axis; thus we desire the exponential
$M=\exp\left(\begin{array}{ccc}
& & \theta\cos\phi\\
& & \theta\sin\phi\\
-\theta\cos\phi & -\theta\sin\phi &
\end{array}\right)=\exp(A).$
This matrix $A$ is diagonalizable, so its exponential $M$ can be computed by hand, which is how I arrived at my answer above.