Struggling to convert cartesian elements into orbital elements in Unity2D + keplerian solution to the two-body problem

Im building a game where i simulate planetary orbits in Unity 2D C#. I previously used newtons laws of gravitation but quickly discovered that rounding errors in floats can lead to massive orbital drift after hundreds or thousands of orbits. This is a big problem for me as I need the simulation to be as consistent as possible every time. Rather than massively complicate things by trying to overwrite the physics system to use doubles instead of floats, I did some research into Keplerian Orbits and the 2 body problem.

Keplerian orbits are far more complex to implement than newtonian and I am struggling to find resources on how to:

• a) accurately simulate a 2D orbit for extended periods of time using Keplerian elements.
• b) convert a rigidbodys current position and velocity (Cartesinal elements) into Keplerian orbital elements so I can make changes to the orbit simulation based on external factors such as adding force to the rigidbody, collisions or another gravitational body attracting the rigidbody (like a moon or something).

Im not very good at understanding long and complex equations intuitively so id appreciate it if you could use meaningful variable names in the answer and try and explain what each block of code does. I want to be able to understand the code rather than just copy and paste a solution as its likely im going to have to come back and make tweaks in the future.

There's a lot going on here, and I'm afraid I don't have enough information about your specific use case to give you a clear-cut answer.

I will specifically address the "for extended periods of time" part of your question (a).

It sounds like all your investigation into Keplerian orbits and the 2-body problem is for the purpose of creating an orbit simulation that does not accrue error over time. Changing all your coordinates to Kepler elements won't solve this problem for you, nor would using doubles instead of floats.

Unless each body in your simulation will only ever be affected by the gravitational force of a single other body (e.g., modeling the solar system exclusively considering the Sun's gravitational force on the planets, and ignoring their forces on the Sun and each other), this won't get you what you're looking for. It sounds like that kind of constraint would be a no-go for you. If you decide that you're OK with this constraint, then you can use an exact analytical solution of Keplerian orbits.

Without that constraint, you're pretty much stuck using numerical solutions. Numerical methods have error. That's an inescapable fact. Real-time physics engines (such as PhysX in Unity) are, in computational physics lingo, mostly glorified Euler method simulations.
There are several ways to go about mitigating the error inherent to Euler simulations. The easiest is to decrease the simulation step size by increasing Unity's physics framerate/fixed timestep. More complex approaches like the Runge-Kutta methods can be implemented by adding a corrective term to the gravitational force you apply to a Rigidbody each frame.

tl;dr: First, decide whether you're OK with the constraints you have to abide by to let you use an exact solution. If you're not, then start looking at error correction in numerical simulation methods.