Looks somewhat similar to what I did a few years ago in a now-abandoned project of mine: https://github.com/rikusalminen/libkepler/blob/master/src/ke...
Did you by any chance, happen to look at my code for reference? Or are they somewhat similar just out of coincidence and similar naming conventions? If you did, in fact, look at my project for reference, I'm glad to have helped!
Just for comparison, here's the orbit determination from initial conditions code from a newer project of mine. It doesn't compute the classical orbital directly, because that tends to add numerical error (ie. first you compute inclination using acos() and a moment later you compute the orientation matrix using cos(inclination)).
https://github.com/rikusalminen/twobody/blob/master/include/...
You can, of course, get the classical elements from the orbit, they are useful for presenting human-readable numbers.