Orbital physics - Maintaining a circular orbit

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().

Code (CSharp):

1. float r = Vector3.Distance(star.Position, planet.Position);
2. float totalForce = -(Constants.G * star.Mass * planet.Mass) / (r * r);
3. Vector3 force = (planet.Position - star.Position).normalized * totalForce;
4. 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<Rigidbody>().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.

Code (CSharp):

1. float initV = Mathf.Sqrt(2f * Constants.G * star.Mass / planet.Position.magnitude);
2. float escapeV = Mathf.Sqrt(4f * Constants.G * star.Mass / planet.Position.magnitude);
3. 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)

 

Hi,

Your use of vis-visa is fine and as you have determined there must be something else going on.

I created the following script and attached it to a RigidBody with a mass of one. (I used the raw position with star assumed at (0,0,0) and G=1 to save myself some typing.)

Code (CSharp):

1. public class UnityPlanet : MonoBehaviour {
2.  
3.     // orbits around a "star" at the origin with fixed mass
4.     public float starMass = 1000f;
5.  
6.     // Use this for initialization
7.     void Start () {
8.         float initV = Mathf.Sqrt(starMass / transform.position.magnitude);
9.         GetComponent<Rigidbody>().velocity = new Vector3(0, initV, 0);
10.     }
11.  
12.     // Update is called once per frame
13.     void FixedUpdate () {
14.         float r = Vector3.Magnitude(transform.position);
15.         float totalForce = -(starMass) / (r * r);
16.         Vector3 force = (transform.position).normalized * totalForce;
17.         GetComponent<Rigidbody>().AddForce(force);
18.     }
19. }

and I get the elliptical orbit you are seeing. (Edited)

This kinda thing is something I have been thinking about a LOT, since I am a week away from submitting a Gravity Physics asset that does elliptical orbits for bodies (showing the orbit in the scene view first), gravitational particles and a ton of other cool gravitational stuff. I'll announce it on nbodyphysics.com and the Asset Forum once it has been released.

EDIT: (Blush) Now that I add an object at (0,0,0) I can see the orbit is actually NOT a circle, but the ellipse you are seeing. As an experiment I put my Gravity Engine asset in the same scene. With the same initial velocity (0,10,0) it DOES move in a circle. This means that the Unity engine physics has gone astray somewhere. I tried all four of the ForceMode settings for AddForce but nothing resolved it.

Cheers,

Peter

 

Thanks for your reply!

It's interesting that I'm not the only one having this problem.

I tried to make the rigidbody kinematic and implement my own, simple physics for this one (apply force for veloctiy, velocity for position, nothing fancy), and I get the same, wrong, result.

I also had a physicist friend of mine look over the code and he couldn't find anything either.

Could you figure out why it works with your Gravity Engine and maybe give a hint? That would be wonderful!

Cheers,
gartoks

 

I tried two other things:
- apply the force in start() as well
- add the velocity in the first fixed update
Neither had an impact.

I do not understand why the Unity physics behaves in this way.

My gravity engine does the acceleration and velocity calculations each FixedUpdate cycle and then updates the positions of the bodies.

A "pro tip" for this is not to use the obvious physics update (e.g. Vnew = Vold + a * dt, Xnew = Xold + v*dt), which is known as Euler's method, but rather use an algorithm that is reversible and conserves energy. Gravity Engine offers Leapfrog (see Wikipedia:Leapfrog Integration) as well as a variable time-step Hermite method. There is also a special purpose integrator for three body solutions.

I have started to get the tutorials ready for asset submission. You can see them at https://www.youtube.com/channel/UCxH9ldb8ULCO_B7_hZwIPvw.

 

The asset "Gravity Engine" is now available in the asset store:
https://www.assetstore.unity3d.com/en/#!/content/62432

 

I know this thread is quite old but as its the only one I could find through Google I would just like to comment that at least in 5.5.1f1, gartoks method works for a circular orbit so the physics issue must have been fixed in an update.

 

Here is my "Force Push On Start" script that has a circularization toggle in the inspector. I use this in Orbital Dogfight to create an asteroid ring.

Code (CSharp):

1. using System.Collections;
2. using System.Collections.Generic;
3. using UnityEngine;
4.  
5. public class ForcePushOnStart : MonoBehaviour {
6.     public float forceAmount = 1000;
7.     public Vector3 pushVec;
8.     private Rigidbody objectToPush;
9.     public float _G = 6674f; // increased gravitational constant
10.     public float _M = 100000f; // temp number : will be replaced by Dominant Gravity Source's Mass
11.     public float _r = 10000f; // temp number : will be replaced by distance between mass centers
12.     public float _F;
13.  
14.     public bool _Circularizing = false; // set this to TRUE in the inspector to create a circular orbit
15.     private GravityObject m_GO; // this is a public class that keeps track of the dominant gravity source / mass // etc for a given object
16.  
17.     public bool randomRotation = false; // can ignore the random raotation functions for this
18.     private float xseed;
19.     private float yseed;
20.     private float zfeed;
21.  
22.     // Use this for initialization
23.     void Start () {
24.  
25.         m_GO = GetComponent<GravityObject>();
26.         objectToPush = gameObject.GetComponent<Rigidbody>();
27.  
28.         StartCoroutine("WaitingForInfo");
29.     }
30.  
31.     private void StartPush()
32.     {
33.         objectToPush.mass = objectToPush.mass * transform.localScale.x;
34.         Vector3 gravity = m_GO.m_DominantGravitySource.transform.position - transform.position;
35.         _M = m_GO.m_DominantGravitySource.GetComponent<Rigidbody>().mass; //the Mass attracting
36.         Vector3 pushDirection = Vector3.ProjectOnPlane(transform.forward, -gravity.normalized).normalized; //gets the correct direction to push for a circular orbit
37.  
38.         _r = gravity.magnitude;
39.         //Debug.Log("Has Radius");
40.  
41.         //Debug.Log("Has Gravity Direction");
42.  
43.         pushVec = pushDirection;
44.  
45.         if (_Circularizing)
46.         {
47.             GetCircularizingForce();
48.         }
49.         else
50.         {
51.             pushVec = pushVec * forceAmount;
52.             objectToPush.velocity = pushVec;
53.         }
54.  
55.         if (randomRotation) // this is can be ignored ..
56.         {
57.             xseed = Random.Range(2f, 89f);
58.             yseed = Random.Range(2f, 89f);
59.             zfeed = Random.Range(2f, 89f);
60.             transform.rotation = Quaternion.Euler(xseed, yseed, zfeed);
61.         }
62.     }
63.  
64.     private IEnumerator WaitingForInfo() //get the dominant gravity source : this is a two body equation
65.     {
66.         yield return new WaitUntil(() => m_GO.m_DominantGravitySource != null); // get the largest mass exerting force on this
67.  
68.         StartPush();
69.     }
70.  
71.     void GetCircularizingForce()
72.     {
73.         _F = Mathf.Sqrt((_G * _M) / _r); // here is the trick .. only the larger Mass's value is included in the equation
74.         pushVec = pushVec * _F;
75.         objectToPush.velocity = pushVec;
76.     }
77. }
78.  

and some more orbit play

asteroid sinusoidal waveform collapse (dont ask) each asteroid is in a circular orbit