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

Postby Julio Jerez » Sat Mar 16, 2019 4:03 pm

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
Moderator
 
Posts: 11040
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

Re: Observing odd behavior with orbital dynamics [SOLVED]

Postby misho » Sat Mar 16, 2019 4:27 pm

Ah, excellent, thanks!
Misho Katulic
CTO, FSX SpacePort
TerraBuilder
www.terrabuilder.com
misho
 
Posts: 552
Joined: Tue May 04, 2010 10:13 am

Re: Observing odd behavior with orbital dynamics

Postby misho » Tue Aug 13, 2019 3:05 am

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 :mrgreen:

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[0][0]);
   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[0]);

   curGravTime = dGetTimeInMicroseconds();
   double timeStep = (double)(curGravTime - m_dPrevGravTime) / 1000000.0f;
   m_dPrevGravTime = curGravTime;

   dvNewPos = dvCurrentPos + dvCurrentVelocity.Scale(timeStep);
   double toOrigin[3] = { 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[0], mass * dGravity * toOrigin[1], mass * dGravity * toOrigin[2], 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 Katulic
CTO, FSX SpacePort
TerraBuilder
www.terrabuilder.com
misho
 
Posts: 552
Joined: Tue May 04, 2010 10:13 am

Re: Observing odd behavior with orbital dynamics

Postby Julio Jerez » Tue Aug 13, 2019 1:12 pm

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
Moderator
 
Posts: 11040
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

Re: Observing odd behavior with orbital dynamics

Postby misho » Tue Aug 13, 2019 1:35 pm

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 Katulic
CTO, FSX SpacePort
TerraBuilder
www.terrabuilder.com
misho
 
Posts: 552
Joined: Tue May 04, 2010 10:13 am

Re: Observing odd behavior with orbital dynamics

Postby Julio Jerez » Tue Aug 13, 2019 2:21 pm

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
Moderator
 
Posts: 11040
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

Re: Observing odd behavior with orbital dynamics

Postby misho » Tue Aug 13, 2019 2:58 pm

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 * dt
dvNewPos = 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 Katulic
CTO, FSX SpacePort
TerraBuilder
www.terrabuilder.com
misho
 
Posts: 552
Joined: Tue May 04, 2010 10:13 am

Re: Observing odd behavior with orbital dynamics

Postby Julio Jerez » Tue Aug 13, 2019 3:24 pm

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
Moderator
 
Posts: 11040
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

Re: Observing odd behavior with orbital dynamics

Postby misho » Tue Aug 13, 2019 3:51 pm

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[0], dMass * dGravity * toOrigin[1], dMass * dGravity * toOrigin[2], 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 Katulic
CTO, FSX SpacePort
TerraBuilder
www.terrabuilder.com
misho
 
Posts: 552
Joined: Tue May 04, 2010 10:13 am

Re: Observing odd behavior with orbital dynamics

Postby Julio Jerez » Tue Aug 13, 2019 5:08 pm

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
Moderator
 
Posts: 11040
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

Re: Observing odd behavior with orbital dynamics

Postby misho » Tue Aug 13, 2019 7:16 pm

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[0][0]);
   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 :roll: : 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 Katulic
CTO, FSX SpacePort
TerraBuilder
www.terrabuilder.com
misho
 
Posts: 552
Joined: Tue May 04, 2010 10:13 am

Re: Observing odd behavior with orbital dynamics

Postby Julio Jerez » Wed Aug 14, 2019 10:30 am

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
Moderator
 
Posts: 11040
Joined: Sun Sep 14, 2003 2:18 pm
Location: Los Angeles

Re: Observing odd behavior with orbital dynamics

Postby misho » Wed Aug 14, 2019 12:47 pm

Sounds good, thanks Julio!

Misho
Misho Katulic
CTO, FSX SpacePort
TerraBuilder
www.terrabuilder.com
misho
 
Posts: 552
Joined: Tue May 04, 2010 10:13 am

Re: Observing odd behavior with orbital dynamics

Postby misho » Thu Aug 29, 2019 6:40 pm

Hi Julio!

Have you had a chance to take a look at that formula?

thanks,
Misho
Misho Katulic
CTO, FSX SpacePort
TerraBuilder
www.terrabuilder.com
misho
 
Posts: 552
Joined: Tue May 04, 2010 10:13 am

Re: Observing odd behavior with orbital dynamics

Postby misho » Wed Oct 30, 2019 5:16 pm

Hi Julio,

I am still struggling with this - I am observing small changes in orbital parameters when the object is not under any influence other than the force of gravity - and that shouldn't be happening. A body's orbital path (an ellipse) should be a constant, given no other outside influences.

Just to recap: In the course of this discussion, you suggested I should use Newton's law of universal gravitation formula. I tried it, and it didn't work, so I reverted back to my original method of using mass*gE (gE being gravitational acceleration for Earth, adjusted for altitude). Since then, I looked at it again, and found that I made a stupid sign error :roll: , and now it works - the code I'm using implements Newton's law of universal gravitation formula, as such:

Code: Select all
void   TBEntity::GetGravityForce(dVector &dvGravityForce)
{
   // Force of gravity portion: Newton's law of universal gravitation

   double dBodyMass = GetCurrentBodyMass();
   dVector dvCurrentPos = GetPosition();
   dVector dvCurrentVel = GetVelocity();
   dVector dvNewPos;

   dvNewPos = dvCurrentPos + dvCurrentVel.Scale(SystemVars->m_SimTimeStep);
   double dRadius = vLength(dvNewPos);
   vNormalize(dvNewPos);

   double dtotalForce = -G*EARTH_MASS*dBodyMass / (dRadius*dRadius);
   dvGravityForce = dvNewPos.Scale(dtotalForce);
}


BTW, I tested both methods and the discrepancy between the returned magnitudes of force vectors is very small, only 0.088%. I'm of course going with what you suggested.

However, I am still getting small errors in orbital path shape. This error gets progressively large when I use rocket engine to shape the orbital path into highly elliptical with a large apogee, usually needed to reach other celestial object (Moon, in my case). What I observe is, after I cut the rocket engine, orbital path keeps changing - the apogee (farthest tip of the orbital ellipse) keeps decreasing, and that shouldn't be happening.

Is there anything I'm doing wrong in the above method? If everything looks good, I could prepare a demo test case for demoSandBox for you to take a look at.

Thanks,
Misho
Misho Katulic
CTO, FSX SpacePort
TerraBuilder
www.terrabuilder.com
misho
 
Posts: 552
Joined: Tue May 04, 2010 10:13 am

PreviousNext

Return to Bugs and Fixes

Who is online

Users browsing this forum: No registered users and 2 guests