I have a really old project of mine dealing with relativistic celestial mechanics. Questions I ask on this site might be related to that. An old paper on that subject matter, which is a bit messy, is available here. In that paper I use the expression:

to simulate orbits. At least I get the Schwarzschild radius right as seen here The innermost stable circular orbit is actually a bit off even if it look ok when you start the simulation from the aphelion of a not too eliptical orbit. You get the expression by plugging:

into:

.

I currently believe that you might get results somewhat more in line with what is expected from the Schwarzschild solution if you use the gamma shown below, but I have not checked it out yet:

Using this "new gammafunction" for pure non-radial motion you get the same expression as the first above. For pure radial motion you instead get:

We can make a superposition of these two solutions that looks like:

I do not know if this expression has any clear benefit, but it least it is based on using an anisotropic velocity of light, just as in GR. I might add some more info on my blog.