## Observing odd behavior with orbital dynamics

Report any bugs here and we'll post fixes

Moderators: Sascha Willems, Thomas

### Observing odd behavior with orbital dynamics

Hi Julio, et.al.

I have a bit of an odd and specific problem, and I think it is on Newton side.

When I launch my rocket into orbit, I follow the "normal" path, by climbing to the desired altitude, leveling the pitch, and burning engine till the orbit is more or less circular. Then - I cut the engines out.

At this point, I would expect the orbit to be "locked" - there are no external forces acting on the rocket, except gravity. The newton engine takes gravity and the velocity and calculates the state vector, which should be a "steady", if not completely circular, orbit.

Except - I don't see that. My Orbital parameters display (both numerical and graphical) shows orbital elements changing. Perigee decreases while apogee increases, and eccentricity slowly increases as well. I have checked, and there are no external forces acting on the body, except gravity. My display is accurate, and is calculated using object's position, so in no way influences the body's motion.

I have checked with stackexchange's Space Exploration community, and they confirmed that if there are no external forces (engines, atmospheric drag, other massive bodies) acting on the body, the orbit should be steady, and I should look into how the physics engine is handling this. The more detailed discussion can be found here.

So - any idea why this might be happening?
misho

Posts: 515
Joined: Tue May 04, 2010 10:13 am

### Re: Observing odd behavior with orbital dynamics

How much does it drift?
That behavior is to be expected.
But is essay to fix.
Basically what happens is that you are applying at non linear force to the body but a simple semi implicit integrator is not enought to reproduce a large scale integration.
This is why nation has a force and torque callback.
There you need to implement a mini integration step
To advance the model in time, the use that value as a correction force.
The engine does that for angular acceleration, but not for gravitational acceleration, for gravity if G is constant, the Simi implicit integrator is accurate enought.

Post you force and torque callback, I can make that correction.
Bear in mind that this is a numerical integration, over long time it will drift, but it can be much better.
Julio Jerez
Moderator

Posts: 10847
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

### Re: Observing odd behavior with orbital dynamics

Notice this is not cheating at all, this is adding a full implicit integration to the position step.
This is because the rocket move at very high speed, so each velocity change in every step is large enought to produce the drift that you see.
Essenctially the truck is to do a newton raps on step in the force and toque callback, to predict the velocity at t + st and use that as the gravity at time t.

This is a very standard integration method, but is only applicable when you know the expression for the derivative.

This is the case for giro, forces so the engine does it,
For a contain g the factor is zero is the velocity is not to big.
But for g = -G / r^2 when v is significant large the correction factor is not negligible, so over time it will make the object drift.
Julio Jerez
Moderator

Posts: 10847
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

### Re: Observing odd behavior with orbital dynamics

Julio Jerez wrote:How much does it drift?
Post you force and torque callback, I can make that correction.
Bear in mind that this is a numerical integration, over long time it will drift, but it can be much better.

Thanks Julio for the prompt response! In the SpaceExploration post, I posted 3 screenshots of my orbital display:

As you can see (time stamp readout on the bottom), in a span of 10 seconds, apoapsis and periapsis drifted about 10 km (in the same direction). It is as if the whole orbital ellipse is shifting. Eccentricity "e" went slightly up as well, and "a" and "b" ellipse parameters are changing slowly. Inclination angle is rock steady.

Sounds good, I will PM you my force and torque callback, check your mail in a few minutes!
misho

Posts: 515
Joined: Tue May 04, 2010 10:13 am

### Re: Observing odd behavior with orbital dynamics

I sent the PM with code more than 2 hours ago, but it is still in outbox - not sure what is holding it up...
misho

Posts: 515
Joined: Tue May 04, 2010 10:13 am

### Re: Observing odd behavior with orbital dynamics

I think it stays in the outbox until Julio reads your pm, then it moves over to sent
Sweenie

Posts: 478
Joined: Mon Jan 24, 2005 7:59 am
Location: Sweden

### Re: Observing odd behavior with orbital dynamics

Got it, thanks!!
misho

Posts: 515
Joined: Tue May 04, 2010 10:13 am

### Re: Observing odd behavior with orbital dynamics

10 km in 10 secund sound like a too big of an error, the must be a bug is the gravity
I took a quick look, and It does no look to me that G is right.

you have this expression

dFloat dGravity = -GE;
dVector gravityForce(mass * dGravity * toOrigin[0], mass * dGravity * toOrigin[1], mass * dGravity * toOrigin[2], 1.0f);

-Ge I assume Is the math of earth / distance from earth center to the rigid body.

so that will be
Ge = M / (r.dot (r)
dir = normalized(r).
g = dir.Scale (Ge);

is that what you are doing? because I did not see when Ge is calculated,
did you assume that r is more of less constant.
Julio Jerez
Moderator

Posts: 10847
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

### Re: Observing odd behavior with orbital dynamics

Julio Jerez wrote:-Ge I assume Is the math of earth / distance from earth center to the rigid body.

Sorry - GE is simply Earth's gravitational acceleration:

Code: Select all
`const double GE            = 9.80655;`

toOrigin is simply a vector from earth center to body's location, normalized.

And the math for this couldn't be simpler: F=m x a, F and a are vectors, mass is scalar. (Or, a is a vector when toOrigin.scale(a))

Code: Select all
`dVector gravityForce(mass * dGravity * toOrigin[0], mass * dGravity * toOrigin[1], mass * dGravity * toOrigin[2], 1.0f);`

Sure, I could have used

Code: Select all
`dVector gravityForce = toOrigin.Scale(mass*dGravity);`

That would be more readable, I guess... Simply, gravityForce is a gravity force vector that points to the center of the Earth. Unless I am missing something huge - that's how I have the force of gravity set up.
misho

Posts: 515
Joined: Tue May 04, 2010 10:13 am

### Re: Observing odd behavior with orbital dynamics

10 km in 10 seconds sounds like a lot, but the distances involved ARE huge - Earth diameter is 12756 km...

Another thing that is worth mentioning is, I have some unused code (which I could re-activate) that does the following (using Newton):

• Prompts me for Latitude, Longitude, Altitude (above ground) and Heading
• Creates velocity vector of the magnitude necessary for the perfectly circular orbit, tangential to the sphere (Altitude + EarthRad), in the direction of Heading
• Applies that velocity vector to the body
The result is, my space object is in a perfectly circular orbit. I remember testing this a few years back (in order to validate Newton for this project) and the space object kept spinning around Earth, no errors orbital decay, or degradation. I remember leaving it over night a few times, and it counted 20 or so orbits (real time) without any problem. I used the same gravity setup as the one I sent you.
misho

Posts: 515
Joined: Tue May 04, 2010 10:13 am

### Re: Observing odd behavior with orbital dynamics

misho wrote:
Julio Jerez wrote:Sorry - GE is simply Earth's gravitational acceleration:
Code: Select all
`const double GE            = 9.80655;`

That will be wrong, and I am guessing the error is so large because the gravity is too strong,
if you are going to compared to real value.

and the altitude of a Satellite, G is probably around 9.2 or something like that which does not sound much but at the speed you are moving it make a difference.
Orbits are not arbitrary each orbit has a specific velocity. so if you have a realistic velocity that you go form some database, by you G is does not math you rocket will probably drift or away of come down to earth.

Anyway I see if I can writ a function for you this Saturday.
Julio Jerez
Moderator

Posts: 10847
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

### Re: Observing odd behavior with orbital dynamics

Julio Jerez wrote:
misho wrote:
Julio Jerez wrote:Sorry - GE is simply Earth's gravitational acceleration:
Code: Select all
`const double GE            = 9.80655;`

That will be wrong, and I am guessing the error is so large because the gravity is too strong,
if you are going to compared to real value.

and the altitude of a Satellite, G is probably around 9.2 or something like that which does not sound much but at the speed you are moving it make a difference.
Orbits are not arbitrary each orbit has a specific velocity. so if you have a realistic velocity that you go form some database, by you G is does not math you rocket will probably drift or away of come down to earth.

Anyway I see if I can writ a function for you this Saturday.

Thanks Julio! If the problem is adjusting gravitational acceleration constant for altitude, I can certainly do that myself and test! At the altitudes I am at, let's say, 400 km, International Space Station, the constant is about 90% less of what it is on the surface of the Earth. I have the formula from Wiki:

g(h) = g(0)(Re/(Re + h))^2

I can certainly plug that in and see if I am getting what I'd expect...
misho

Posts: 515
Joined: Tue May 04, 2010 10:13 am

### Re: Observing odd behavior with orbital dynamics

Julio,

I implemented the adjusted gravity constant, as such:

Code: Select all
`double dAltitude = length(toOrigin) - EARTH_RAD_E;double dGravity = -GE*pow((EARTH_RAD_E / (EARTH_RAD_E + dAltitude)), 2);`

and it WORKED! Thank you!

The numbers are rock steady, except for both apogee and perigee, which still drift, but only by about 10 meters, every 3 seconds... which is WAY less than before, and probably within the numbers you expected to see when you were explaining the potential cause. So - now that I have it fixed, I would still love to have your fix that corrects this small drift, but no rush. Thanks again!
misho

Posts: 515
Joined: Tue May 04, 2010 10:13 am

### Re: Observing odd behavior with orbital dynamics

The numbers are rock steady, except for both apogee and perigee, which still drift, but only by about 10 meters, every 3 seconds...

ah that's a long the line of the error I was expecting.
from here you can reduce that erro even further by doing this secund trick.

the equation of motionof a particle under the action of gravity and not otherhe forces is
m * a = f
a = dv/dt = f / m
we know f / m = -G
where G ar any give position and time is give by -G0/(r * (x + h)

there are two ways to integrate that equation: forward euler and backward euler
forward euler is what the engine does for the velocity and backward for the position.
this makes sence because G close to earth is constantt and does not chanage direction no matter how fast you move near earth.

foward euler is expresed by this
(v1 - v0) / dt = G(t,x)

v1 = v0 + G(t,x) * dt

this is G is the derivative at time and location t, and you mutiply by the dt and add that to the current velocity. graphically this is like this.
Untitled.png (2.71 KiB) Viewed 310 times

as you can see when g(t,x) is highly non linear or v0 is very high the change of G(t,x) * dt is quite large and as you can see in the picture v1 drift fromr the ideal v1 that is project on the curve.

let us see the backwork euler now.

this is given by using the velocity not at time v0 but at time v(-dt) and project forward by makinge a change of variable, the final equation is this.

v1 = v0 + G(t + dt, x + dx) * dt

this is, using the value of the derivative at one dt in the future as if it was the derivative at time t.

this is far more complex for function of multiple variables because the way G(t +dt, x, dx) is calcuatione is by expanding on a taylor series and negleting the high order terms.
but for the of an or it busy, is quite eassy to estime G at t + dt.

besially we know that not matter how fast a rocket can travel, r does not changes significantly in magnitude so the linear derivative can be ignored.
the direction however does change significanlly at the speed of a satelity, for example an stationary one has 24 houd period, and rocket move even faster.
graphically the backword euler looks like this.
Untitled1.png (4.46 KiB) Viewed 310 times

as you can see in backward euler the velocity at time t + dt end up below the ideal curve, so backward Euler is stable because in loses energy. but that picture is exaggerated, in reality is quite acurate and stable for smooth function like your case.

what remind is how to get G(t +dt, x+dx)

but that is quite eassy, your vehicle is in orbist and has a radius, you know that orbit move on plane, therefore you can stimate G(t+dt, x+dx) by taking the position of the satelite a t0, that is the position of the com of the rigid body.

them get the linar velocity of the rigid body, and add the value v * dt
this is let x0 be the position vector
x1 = x0 + v0 * dt
now with that position, you know the radious is r + h

so you calculate a unit vector for position x1 and multiply by the mangity of r + h, that will be the radius of the rocket at g (t+dt, x+dx)

now you use this value of G to calcuate the force that will be applied to the satellite.

do it right and the error should reduce to a very small value. almost zero.
Julio Jerez
Moderator

Posts: 10847
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

### Re: Observing odd behavior with orbital dynamics

Ok, great, thanks Julio!

The math looks hairy but doable (as soon as I see derivatives and integration I cringe because I haven't done a lot of this in code). I'll see if I can handle it, and let you know if I need any help!

Misho
misho

Posts: 515
Joined: Tue May 04, 2010 10:13 am

Next