- Home /
circular orbit weird with physics
I've been looking through a lot of old threads on how people manage to simulate realistic orbits and gravity using newton's law of gravitation as well as the velocity required to circularize an orbit. I found a lot of information but when i try to use it i encounter a few errors. I am using rigitbody 2d and applyForce.
1) orbits decay. I use the formula F = Gm / r^2 but instead of using a gravitational constant and mass i use a float (0.001f) for testing, in a fixedupdate, which i apply using addforce. The orbits end up nice if I initially add a rough circularization force, but slowly the orbit's apogee increases after each orbit. If i start at a distance of 2000 from the barycenter by the time it returns to the same angle it will be around 2002, and it continues to increase slowly for each orbit. I do not know why it does this and whether it has to do with floating point errors but I was expecting the apogee to remain 2000 consistently.
2) circularization force can't be calculated. I use v = Sqrt(0.001f * r) and in big orbits the force generated is too small while closer, smaller orbits, the circularization force is too big. I don't know if its some inconsistency with my units (its supposed to be G*mass but i just use a custom gravity float)
I feel like some of this isn't possible unless there's an inconsistency in unity's physics. I've read in various threads that people also end up with elliptical orbits and blame unity's physics, but others say it was patched recently
$$anonymous$$aybe the orbit decay is a result of the time step size? You could try setting Time.fixedDeltaTime to a smaller value and see if that makes it better.
Answer by Bunny83 · Sep 09, 2019 at 12:08 AM
Well, nothing you have experienced was unexpected or is a bug or anything like that. When you apply Newtons law of gravitation every frame you essentially performing dead reckoning in discrete time steps. This will never be "accurate" or "stable" over time since even the slightest error will accumulate.
Keep in mind that the real world does not "run" in discrete time intervals and also doesn't have any numerical precision limit on any involved figures. Also keep in mind that orbits are almost always elliptical and not circular.
If you want a stable orbit simulation you can not just use newton's law and integrate it over time with singe precision floating point values. You would need to specify your exact orbit and use absolute figures (like orbital elements)to specify your position on that orbit.
Note that almost all time discrete physics systems will have this issue since any change between two time steps is considered linear. Any non linear process will accumulate errors over time. This error gets smaller if you do smaller time steps, however the realworld does actual continuous integration. So in theory infinitely small timesteps. Though even when you have unlimited time to calculate infinite steps you would still accumulate errors due to number precision limits of floating point numbers.
ps: 2 units off when having a 2000 unit radius is only an error of 0.1% after how many iterations?
Answer by sjhalayka · Sep 09, 2019 at 05:16 AM
Which integration method are you using? I find that when using Euler integration, strange orbit behaviour occurs, like non-constant speed when in a circular orbit.
Here is some C# code for doing symplectic or Runge-Kutta integration:
void proceed_euler(ref Vector3 pos, ref Vector3 vel, Vector3 ang_vel, float dt)
{
vel = vel + acceleration(pos, vel, ang_vel) * dt;
pos = pos + vel * dt;
}
void proceed_symplectic4(ref Vector3 pos, ref Vector3 vel, Vector3 ang_vel, float dt)
{
float cr2 = Mathf.Pow(2.0f, 1.0f / 3.0f);
float[] c = new float[4]
{
1.0f / (2.0f*(2.0f - cr2)),
(1.0f - cr2) / (2.0f*(2.0f - cr2)),
(1.0f - cr2) / (2.0f*(2.0f - cr2)),
1.0f / (2.0f*(2.0f - cr2))
};
float[] d = new float[4]
{
1.0f / (2.0f - cr2),
-cr2 / (2.0f - cr2),
1.0f / (2.0f - cr2),
0.0f
};
pos += vel * c[0] * dt;
vel += acceleration(pos, vel, ang_vel) * d[0] * dt;
pos += vel * c[1] * dt;
vel += acceleration(pos, vel, ang_vel) * d[1] * dt;
pos += vel * c[2] * dt;
vel += acceleration(pos, vel, ang_vel) * d[2] * dt;
pos += vel * c[3] * dt;
// vel += acceleration(pos, vel, ang_vel) * d[3] * dt; // last element d[3] is always 0
}
void proceed_rk4(ref Vector3 pos, ref Vector3 vel, Vector3 ang_vel, float dt)
{
const float one_sixth = 1.0f / 6.0f;
Vector3 k1_velocity = vel;
Vector3 k1_acceleration = acceleration(pos, k1_velocity, ang_vel);
Vector3 k2_velocity = vel + k1_acceleration * dt * 0.5f;
Vector3 k2_acceleration = acceleration(pos + k1_velocity * dt * 0.5f, k2_velocity, ang_vel);
Vector3 k3_velocity = vel + k2_acceleration * dt * 0.5f;
Vector3 k3_acceleration = acceleration(pos + k2_velocity * dt * 0.5f, k3_velocity, ang_vel);
Vector3 k4_velocity = vel + k3_acceleration * dt;
Vector3 k4_acceleration = acceleration(pos + k3_velocity * dt, k4_velocity, ang_vel);
vel += (k1_acceleration + (k2_acceleration + k3_acceleration) * 2.0f + k4_acceleration) * one_sixth * dt;
pos += (k1_velocity + (k2_velocity + k3_velocity) * 2.0f + k4_velocity) * one_sixth * dt;
}
Your answer
![](https://koobas.hobune.stream/wayback/20220612211842im_/https://answers.unity.com/themes/thub/images/avi.jpg)
Follow this Question
Related Questions
calculate the future position of an object in orbit 1 Answer
Natural rotation of orbiting object 3 Answers
Collision and physics: asteroid 1 Answer
Artificial Gravity 2 Answers