Coverting the keplerian elements to cartesain vectors ( posistion and velocity) is relatievly simple by using rotation matrices show in the document here: https://downloads.rene-schwarz.com/download/M001-Keplerian_Orbit_Elements_to_Cartesian_State_Vectors.pdf
Though I was wondering would I be possible to convert the elements into vectors using quaternions instead of matrix rotation? If so how would you do it?
The rotation method uses a 3-1-3 Euler rotation: rotate by Z, then X, then Z again. For each rotation, one axis (in the new frame) is fixed. A quaternion based rotation requires finding a vector which is immobile between the initial frame and the final frame (see here for an excellent visualization and explanation of quaternions). One method to compute the quaternion you need, is to compute the rotation matrix as you normally would and extract the components of the quaternion corresponding to that rotation using the method described here (but ignore the yaw-pitch-roll notation because you're dealing with a 3-1-3 rotation not a 3-2-1 rotation).
The conversion to Keplerian orbital elements to Cartesian orbital element is not typically done using rotation matrices. The reason for this is that it's an incomplete method: it does not work for hyperbolic orbits, and I don't think it works for near circular orbit (because the eccentricity is ill-defined, so the eccentric anomaly ($E$) will also be ill-defined). Instead, the method used by GMAT and Nyx (and surely others) consists in using the semi-parameter to calculate the radius and velocity vectors, cf. this explanation and the associated algorithm.
