# Orbital physics - Maintaining a circular orbit

Discussion in 'Physics' started by gartoks, May 8, 2016.

1. ### gartoks

Joined:
May 7, 2016
Posts:
2
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;
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 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)

Last edited: May 8, 2016
2. ### PeterMu

Joined:
Oct 1, 2014
Posts:
47
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;
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

Last edited: May 11, 2016
3. ### gartoks

Joined:
May 7, 2016
Posts:
2

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

4. ### PeterMu

Joined:
Oct 1, 2014
Posts:
47
I tried two other things:
- apply the force in start() as well
- add the velocity in the first fixed update

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.

Joined:
Oct 1, 2014
Posts:
47
6. ### Krazysh01

Joined:
Apr 28, 2014
Posts:
1
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.

7. ### ristophonics

Joined:
May 23, 2014
Posts:
32
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;
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

Last edited: May 8, 2020