## Observing odd behavior with orbital dynamics

Report any bugs here and we'll post fixes

Moderators: Sascha Willems, Thomas

### Re: Observing odd behavior with orbital dynamics

You do not have to do any complex math.
All you have to understand is the idea of backward Euler integration. Which is simply figure out the derivation at a time in the future and use it as the derivative at present time.

In the case of an orbist is quite simple.
1 in a temporary step integrate the particle by the current velocity and position.
2 use this new position to calculate the ideal orbit location, wich you know because it must be of some known radius r +h,
3 With the new position calculate the gravity force and pass that as the force in the force callback.

The effect of this is a tiny correction that does not seem much, but over long time it is what compensate for the numerical errors accumulated by the forward integrator in the engine.
Julio Jerez
Moderator Posts: 10954
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

### Re: Observing odd behavior with orbital dynamics [SOLVED]

Ah, excellent, thanks!
misho

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

### Re: Observing odd behavior with orbital dynamics

Hi Julio!

It's been a while since I visited this (been busy with demo installer) but I finally had a chance to implement this. Thanks again for the suggestion! I had to read it a few times but in the end it made perfect sense Now, while in orbit, the orbital parameters are rock steady (as they should be) EXCEPT, apogee and perigee seem to be slowly decreasing/increasing, at a rate of about 5m every second.

Here is my code - I think I got everything right.

Code: Select all
`void   TBEntity::GetGravityForce(dVector &dvGravityForce){   dFloat Ixx;   dFloat Iyy;   dFloat Izz;   dFloat mass;   dMatrix bMatrix;   unsigned64 curGravTime = 0;   NewtonBodyGetMatrix(nBody, &bMatrix);   NewtonBodyGetMass(nBody, &mass, &Ixx, &Iyy, &Izz);   dVector dvCurrentPos = dVector(bMatrix.m_posit.m_x, bMatrix.m_posit.m_y, bMatrix.m_posit.m_z);   dVector dvCurrentVelocity;   dVector dvNewPos;   NewtonBodyGetVelocity(nBody, &dvCurrentVelocity);   curGravTime = dGetTimeInMicroseconds();   double timeStep = (double)(curGravTime - m_dPrevGravTime) / 1000000.0f;   m_dPrevGravTime = curGravTime;   dvNewPos = dvCurrentPos + dvCurrentVelocity.Scale(timeStep);   double toOrigin = { dvNewPos.m_x, dvNewPos.m_y, dvNewPos.m_z };   double dAltitude = length(toOrigin) - EARTH_RAD_E;   normalize(toOrigin); // this gives a vector pointing to origin, of unit length.   // calculate gravity constant adjusted for altitude:   double dGravity = -9.80665*pow((EARTH_RAD_E / (EARTH_RAD_E + dAltitude)), 2);   dvGravityForce = dVector(mass * dGravity * toOrigin, mass * dGravity * toOrigin, mass * dGravity * toOrigin, 1.0f);}`

So I was wondering, shouldn't apogee/perigee values be constant? I thought the above approach would fix those small errors as well...

Misho
misho

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

### Re: Observing odd behavior with orbital dynamics

I had to look for the definition of apogee and perigee in Google, for what I can see they are just points in an orbit that are given special meanings.
In the case of Earth they are the point father away and the point closer around the sun orbits.

Are you saying the orbits is decain 5 meter per second? That sounds like a huge error, not sure how those names play a role.

In any case if you get an stable orbits, that's all that count, the orbits does not has to be circular.
Julio Jerez
Moderator Posts: 10954
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

### Re: Observing odd behavior with orbital dynamics

Julio Jerez wrote:I had to look for the definition of apogee and perigee in Google, for what I can see they are just points in an orbit that are given special meanings.
In the case of Earth they are the point father away and the point closer around the sun orbits.

Hi Julio!

Sorry, yes, Apogee and perigee are astronomical terms equivalent to what apoapsis and periapsis are in mathematics terms. Basically, they are the farthest and closest points on the ellipse, relative to one of the ellipse's focii. They are VERY important to astrodynamics, because that is where the spacecraft fire their engines when they change the shapes of orbits. For example, Apollo had to wait until it reached apogee to fire its engines and change the shape of their Earth orbit to intersect it with Moon.

Julio Jerez wrote:Are you saying the orbits is decain 5 meter per second? That sounds like a huge error, not sure how those names play a role.

In any case if you get an stable orbits, that's all that count, the orbits does not has to be circular.

No, the orbit is stable, it does not decay, but apoapsis and periapsis change, which means that the orbit keeps changing shape (from circular to egg-shaped, then back to circular). For example, if the orbit is a perfect cirlce, apoapsis and periapsis are the same values. However, in spaceflight, that's almost never the case, and in my case, apoapsis is 402 km, and periapsis is 397 km. So, for example, as the spacecraft orbits, the apoapsis drops to 400 km and periapsis rises to 399 km... and then back to 402 and 397 km.

This oscillation is what I am worried about... I'm almost certain that those values shouldn't be changing, and that there is still an error somewhere in the calculations...

Misho
misho

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

### Re: Observing odd behavior with orbital dynamics

However, in spaceflight, that's almost never the case, and in my case, apoapsis is 402 km, and periapsis is 397 km. So, for example, as the spacecraft orbits, the apoapsis drops to 400 km and periapsis rises to 399 km... and then back to 402 and 397 km.

but that's not really a problem. in fact that's a cool result.
basically you are achieving that orbit with not extract information other that the instant velocity of the space craft.

I did not check the function, but I asume your code asume the predict velocity is based of the radios of the orbit one time in the future to be circular.
the error you are getting is a numerical integration typical of backed Euler.
essentially you are trying to integrate a circular orbits, and the orbits is osculating around that ideal path, but stable on average but is not exact at each step.
the length of one orbit is more or less 20 km, and for what you say the error is about +- 2 meters around the fix points. that's a remarkably small error for such long trajectory.
that's a remarkably small error for such a huge length, considering that the spacecraft is simple drifting

to see if what I am saying is even in the ball park of been correct,
try increasing the simulation sub step and see if the error decreases.

this is not passing small time step to newton update, there is a function that set the internal sub steps
Julio Jerez
Moderator Posts: 10954
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

### Re: Observing odd behavior with orbital dynamics

Julio Jerez wrote:but that's not really a problem. in fact that's a cool result.
basically you are achieving that orbit with not extract information other that the instant velocity of the space craft.

Correct - and also, I am calculating the orbital parameters using a state vector, using this method.

Julio Jerez wrote:I did not check the function, but I asume your code asume the predict velocity is based of the radios of the orbit one time in the future to be circular.
the error you are getting is a numerical integration typical of backed Euler.
essentially you are trying to integrate a circular orbits, and the orbits is osculating around that ideal path, but stable on average but is not exact at each step.

Correct- I am using this code to predict POSITION (not velocity), based on what you suggested:

Code: Select all
`// x1 = x0 + v0 * dtdvNewPos = dvCurrentPos + dvCurrentVel.Scale(timeStep);`

where timeStep is one physics cycle. I then use dvNewPos to calculate more accurate gravitational constant.

Julio Jerez wrote:the length of one orbit is more or less 20 km, and for what you say the error is about +- 2 meters around the fix points. that's a remarkably small error for such long trajectory.
that's a rammable small error for such a huge length, considering that the spacecraft is simple drifting

to see if what I am saying is even in the ball park of been correct,
try increasing the simulation sub step and see if the error decreases.

this is not passing small time step to newton update, the is a function that set the internal sub steps

Ok, to put correct values into perspective, the apogee and perigee values are expressed as kilometers above earth surface. In reality, the orbit is a slightly flattened circle, whose size is EARTH_RAD + 400 +/- 4 km, so, the orbit SHAPE oscillates between

Apoapsis (Earth rad + Apogee): 6370 + 402 = 6772 km
Periapsis (Earth rad + Perigee): 6370 + 397 = 6767 km
(a bit more elliptical)

and

Apoapsis (Earth rad + Apogee): 6370 + 400 = 6770 km
Periapsis (Earth rad + Perigee): 6370 + 399 = 6769 km
(a bit less elliptical)

And this oscillation's period is about 90 minutes - the time to complete a full orbit at this altitude. If these errors are within the numerical integrator limits, that's perfectly fine with me. However - I was observing the same errors BEFORE I put in the suggested improvement (backwards Euler) - so, this extra step didn't seem to improve things - and that's why I was wondering if I understood the method correctly and if my code was correct...

Misho
misho

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

### Re: Observing odd behavior with orbital dynamics

backward euler (implicit) is draw numerical errors from in integration step because if uses the derivation one step in the future, in this case because you deal with an orbits, one step in the future can be predicted by using the current velocity take a tentative step and rejecting the position onto the orbit. all those assumption aren't 100% correct,
I assume the orbits is circular, the prediction is linear, and so on. but that assumption sound that the error introduced is small.

forward euler (explicit) use the state (velocity ans position) to calculate the next position the error introduce by this method is much larger because essential the objective muve tangential to the trajectory.

I look at your function and I do no know if this is right
Code: Select all
` // calculate gravity constant adjusted for altitude:   double dGravity = -9.80665*pow((EARTH_RAD_E / (EARTH_RAD_E + dAltitude)), 2);`

I though that the force is just F(t+dt) = G * M1 * M2 / r(t +_ dt)^2
which we are assuming to be circular

where
G = 6.67408e-11N m2/kg2
M1 is mass of the planet
M2 is the mass of the rocket
r (r + dt) is the estimated distance from the center of earth to the center of the craft at t + dt

so all the exercise reduces to to calculate r (t + dt)

where are you getting that value of G from? it seems a fitting approximation function.
Julio Jerez
Moderator Posts: 10954
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

### Re: Observing odd behavior with orbital dynamics

Julio Jerez wrote:I look at your function and I do no know if this is right
Code: Select all
` // calculate gravity constant adjusted for altitude:   double dGravity = -9.80665*pow((EARTH_RAD_E / (EARTH_RAD_E + dAltitude)), 2);`

I though that the force is just F(t+dt) = G * M1 * M2 / r(t +_ dt)^2
which we are assuming to be circular

where
G = 6.67408e-11N m2/kg2
M1 is mass of the planet
M2 is the mass of the rocket
r (r + dt) is the estimated distance from the center of earth to the center of the craft at t + dt

so all the exercise reduces to to calculate r (t + dt)

where are you getting that value of G from? it seems a fitting approximation function.

dGravity is not force, it is a gravitational constant calculation that takes altitude into account. I used the theory/formulae from this source.

The gravity FORCE calculation is the next equation:

Code: Select all
`dvGravityForce = dVector(dMass * dGravity * toOrigin, dMass * dGravity * toOrigin, dMass * dGravity * toOrigin, 1.0f);`

I understood that the formula

F(t+dt) = G * M1 * M2 / r(t +_ dt)^2

is used only when M1 and M2 are of comparatively same masses, as in 2 celestial bodies. I thought that in the case of Earth and a rocket, rocket has a negligible mass compared to Earth, and the formula is not applicable. So - the only thing I needed was an altitude-adjusted g value and direction, and use that for deriving the force vector using F = ma formula (which is how I calculated dvGravityForce )

Misho
misho

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

### Re: Observing odd behavior with orbital dynamics

is used only when M1 and M2 are of comparatively same masses, as in 2 celestial bodies.

well I do not know who said that. the formula F(t) = G * M1 * M2 / r(t)^2
is in fact a universal law of physics not just for gravity force, but also for optics, acoustics, magnetic fields, electric fields and even particle physics.

In the case of a gravitation field as long as the masses can be represented by a point mass, the equation seem valid.
I know there will be many experts from places like stack overflow and gamedevnet that will come back with quantum physics, general gravitation, vending for space and * like that, to those people I can only say this engine does not deal with that.

let us clarify something. in the equation G * M1 * M2 / r(t)^2

we assume that the orbit is mostly circular, so when this is that case the estimation at f at t = t0 + dt
is simply the projection of that position one step in the future to the circulatory trajectory, it works because you have the knowledge that the trajectory is approximated circular. therefore it only change direction but not much in magnitude

for a generic method what to do is that you expand

f (t+dt) = G * M1 * M2 / r(t + td)^2 as an algebraic polynomial expansion and and you discard the higher order term of dt^n
which is what I suspect the expression you are using is.

please try using the term I suggest, see how that goes. it should work as long as your orbits are close to circular. it will also work for elliptic orbits as a long as you know what orbit you want to accomplish.
for any other initial condition value problems, it will not produce the desired result, you need use he expansion but there is not guarantee of a stable orbit.
Julio Jerez
Moderator Posts: 10954
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

### Re: Observing odd behavior with orbital dynamics

Julio Jerez wrote:please try using the term I suggest, see how that goes. it should work as long as your orbits are close to circular. it will also work for elliptic orbits as a long as you know what orbit you want to accomplish.
for any other initial condition value problems, it will not produce the desired result, you need use he expansion but there is not guarantee of a stable orbit.

Thanks Julio! Now, I'm rather confused. In the first part of this thread (back in March), you suggested that all I need to do is to adjust my gravitational constant g to account for altitude, which I did, and everything got fixed, except for those small inaccuracies with apogee and perigee values. Now, you're suggesting not to use that force formula, but to use F(t) = G * M1 * M2 / r(t)^2... where g is not even used?

But anyway - I did implement that formula, here is the code:

Code: Select all
`void   TBEntity::GetGravityForce2(dVector &dvGravityForce){   dMatrix bodyMatrix;   unsigned64 curGravTime = 0;   NewtonBodyGetMatrix(nBody, &bodyMatrix);   double dMass = GetCurrentBodyMass();   dVector dvCurrentPos = dVector(bodyMatrix.m_posit.m_x, bodyMatrix.m_posit.m_y, bodyMatrix.m_posit.m_z);   dVector dvCurrentVel = GetVelocity();   dVector dvNewPos;   // get time step   curGravTime = dGetTimeInMicroseconds();   double timeStep = (double)(curGravTime - m_dPrevGravTime) / 1000000.0f;   m_dPrevGravTime = curGravTime;   dvNewPos = dvCurrentPos + dvCurrentVel.Scale(timeStep);   double dRadius = vLength(dvNewPos);   vNormalize(dvNewPos);   double dtotalForce = G*EARTH_M*dMass / (dRadius*dRadius);   dvGravityForce = dvNewPos.Scale(dtotalForce);}`

The result is, actually, really bad : The Apogee starts climbing really fast and Perigee starts decreasing. After a few seconds, It shows that I am in a highly elliptical orbit, and the return trajectory will crash me back to Earth...

The orbit I start out with for this scenario is highly circular - however, this formula would not be suited for other cases where the orbit is highly elliptical (such as Beyond Earth Orbit travel, or launch)

So, yeah, I have no idea what I'm doing wrong with the above formula, and I'm willing to try anything you suggest...

Misho
misho

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

### Re: Observing odd behavior with orbital dynamics

The value G*M1/r^2 is 9.8 on the surface of the planet.

But this assume objects move slow enought that r the radios does not changed enough to make a difference, but if you have an object that move fast enought that can go from one side of the planet to the other, the the radios change both in magnitude and direction, you can no longer assume 9.8

I see other thing incorrect in that function, it will edit it lator, to what I think should be.
Julio Jerez
Moderator Posts: 10954
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

### Re: Observing odd behavior with orbital dynamics

Sounds good, thanks Julio!

Misho
misho

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

Previous

Return to Bugs and Fixes

### Who is online

Users browsing this forum: No registered users and 1 guest