Hi everyone!

I’m currently working on a 2-body orbital physics simulation and came across a problem involving the orbital physics:

I am using Newtons 2nd Law (F = -G * M *m / r^2) to calculate the gravitational forces for each planet (currently only one for simpler testing, but works with multiple). Then I apply them to the planet’s rigidbodies using the standard AddForce(…) method in FixedUpdate().

```
float r = Vector3.Distance(star.Position, planet.Position);
float totalForce = -(Constants.G * star.Mass * planet.Mass) / (r * r);
Vector3 force = (planet.Position - star.Position).normalized * totalForce;
planet.GetComponent<Rigidbody>().AddForce(force);
```

This works fine, producing the expected results.

My problem surfaces here:

I am trying to calculate the initial velocity required to keep the planets on a circular orbit by using the Vis-Viva equation:

v = sqrt(G * M * (2/r - 1/a))

(derived from setting centripetal force equal to gravitational force)

Simplified for a circular orbit:

v = sqrt(G * M / r).

Using this equation my planet is too slow to keep a circular orbit though.

Note: Position is just a simplification for transform.position. Mass a simplification for GetComponent().mass

Now: The highest velocity allowing for any orbit is the escape velocity v_e = v * sqrt(2).

After a bit of trying I figured out that I need a velocity of sqrt(2) * v for a circular orbit. Basically the escape velocity.

My escape velocity is scaled by sqrt(2) too, making it sqrt(4 * G * M / r).

If I use these scaled values everything works as intended. I also tested with different masses, radii, gravitational constants, etc. Same results.

Also: All parameters are scaled to reasonable values, so no rounding errors from too large or too small numbers. But that shouldn’t affect the result itself if the equations are right.

For completeness: star.Mass = 1000, planet.Mass = 1, r = 10 (star is at 0,0,0; planet at 10,0,0), G = 1

The initial velocity is set in the Start() method once.

```
float initV = Mathf.Sqrt(2f * Constants.G * star.Mass / planet.Position.magnitude);
float escapeV = Mathf.Sqrt(4f * Constants.G * star.Mass / planet.Position.magnitude);
planet.GetComponent<Rigidbody>().velocity += new Vector3(0, initV, 0);
```

I know I’m currently orbiting in the x-y plane, so no error there

Question: Why is it scaled by sqrt(2)? I could just scale it but I want to get it right. Anyone have an idea? I am relatively new to Unity so I’m not sure if it’s maybe something with Unity’s physics engine that creates that problem.

It’s also a possibility that I’m just stupid and overlooked something obvious

Thanks in advance,

gartoks

EDIT:

OK, it MUST be something with Unity. I calculated the eccentricity from the initial position and velocity. With my scaled version (which is actually the escape velocity) I get e = 1, with the velocity as it is supposed to be I get 0, and 0 is the correct answer.

Formula used btw: **e** = ((v^2 - G * M / r)***r** - (**r** * **v**) * **v**) / (G * M)

(Bold letters are vectors, same letter not-bold is the magnitude of that vector)