plotting stars in cartessian coordinates using the stars ra and dec produces the wrong location of the stars
A photograph of stars is essentially a
gnomonic projection
of part of the
celestial sphere.
If $\vec P$ is computed for the photograph center (RA, Dec) in the same way as $\vec S$,
and the north celestial pole is
$\hat{\mathrm k} = (0, 0, 1)^{\mathrm T}$,
then the photograph has basis vectors perpendicular to $\vec P$ and each other,
$$\hat{\mathrm u} = \frac{\vec P \times \hat{\mathrm k}}{\| \vec P \times \hat{\mathrm k} \|},
\hat{\mathrm v} = \frac{\hat{\mathrm u} \times \vec P}{\| \hat{\mathrm u} \times \vec P \|}$$
and you can project $\vec S$ onto the photograph as
$$\left( \frac{\vec S \cdot \hat{\mathrm u}}{\vec S \cdot \vec P},
\frac{\vec S \cdot \hat{\mathrm v}}{\vec S \cdot \vec P} \right)$$
scaling as needed.