A specific setup problem

[misho]

Hi guys,

If I read you correctly, I am already doing that. I have a gravity field defined as a vector pointing to the centre of Earth. All the bodies are orbiting the Earth properly, however, it is a stationary, non-rotating earth.

It is like this: you are standing on the edge of a still merry-go-round, get off it, and start running until you get to a certain speed of, let's say, 10 rpm.

Now, if you were on the edge of a spinning merry go round, and you get off, and start running in the direction of the spin, you will reach those 10 rpm faster, because you already had speed when you got off the merry-go-round.... and conversely, if you started running in the opposite direction, you would have to run a bit longer to reach 10 rpm in the other direction, because you have to overcome the initial speed of the merry-go-round.

So far, are we on the same page?

Now - In the approaches and equations you guys suggested, I don't see anywhere any factoring for the spinning earth, and I don't quite see how this would "resolve itself"...

[Julio Jerez]

I have a gravity field defined as a vector pointing to the centre of Earth. All the bodies are orbiting the Earth properly, however, it is a stationary, non-rotating earth.

maybe my English get in the way of the explanation.

there is not such thing as stationary, this is like saying on the surface of the earth find a place when the gravity is zero.

in your world all bodies are under an acceleration field given by.

Code: Select all

a = -g - 2 * corssProduct (EarthOmega, BodyLinalVelicty) - crossproduct (EarthOmega, crossProduct (EarthOmega, BodyLinalVelicty)


in you call back you calculate that value and you multiplied by the mass of the body, regarless where the body is.
In fact you can put a bunch of particle and apply that force, to each of then and you can see tornadoes for orbit.

[misho]

Code: Select all

a = -g - 2 * corssProduct (EarthOmega, BodyLinalVelicty) - crossproduct (EarthOmega, crossProduct (EarthOmega, BodyLinalVelicty)


Ok - I see. "-g" is the only thing I am applying at the moment, correct? that's my gravity vector pointing
to the earth centre.

Also, EarthOmega, in this particular form, is assumed to be around z-axis (Earth pole-to-pole axis), correct?

If I assumed those correct, I'll put that in my callback and report on the results... next year

[Julio Jerez]

I made a mistake is the equation I wrote correction, the last term, is the position of your space ship relation of the earth origing. here is the correct equation.

Code: Select all

a = -g - 2 * corssProduct (EarthOmega, BodyLinalVelicty) - crossproduct (EarthOmega, crossProduct (EarthOmega, BodyPositionReletiveToErthCenter )

remember g is not the 9.8 it is the general gravitational newton force.

Code: Select all

dir = r / sqrt (dot(r, r));
g is (dir.Scale (G * m1 * m2 / r^2))


yes you need to chose an inertial frame of reference.
you can chose an earth center frame with the y axis aligned to the line the connect the north and south pole. That may be restrictive because it will limit your simulator to low earth orbit objects.
a more general inertial frame can be the sun a the origin and the up vector the perpendicular to the earth orbit around the sum, but in that model you can simulate the seasons, as the planet moves around the sun, and you can even simulate more very complex stuff like focal points and other stuff.
in this model the earth will be a big static object the moves around the sun.
this model may be to complex to start but boy it will be something that even NASA will envy.

let us say you use the first coordinate system.
in this system your code that equation in your force and torque callback, the earth static is a static body with a fixed angular velocity that you integrate manually each frame.

if later you want to move to a more advance system that is also easy but we need the help of a symbolic derivative program like mathematics to calculate the forces acting on the body.

in this model all we nee to do is to come up with the expression that calculate the total kinetic an potential energy of a body in space on a fixed coordinate frame of reference centered on the sun and we calculate the Langragian derivatives to get the forces and torques acting on that body.
BTW the Langragian formulation is the method use by big space agencies like NASA to calculate trajectory or simulation of celestial bodies. but that a different problem.

for now see if you can get model one working and them lots of things will make sense to you that do not seem logical now.


for the fist model if

[misho]

remember g is not the 9.8 it is the general gravitational newton force.

Right - I understand. Here is my "version" of what you wrote, much lengthier but I think it does the same thing

Code: Select all

   NewtonBodyGetMatrix(body, &matrix[0][0]);
   NewtonBodyGetMassMatrix (body, &mass, &Ixx, &Iyy, &Izz);
   double toOrigin[3] = {matrix.m_posit.m_x, matrix.m_posit.m_y, matrix.m_posit.m_z};
   normalize(toOrigin);
   double fRadius = sqrt( matrix.m_posit.m_x*matrix.m_posit.m_x + matrix.m_posit.m_y*matrix.m_posit.m_y + matrix.m_posit.m_z*matrix.m_posit.m_z );
   dFloat dGravity = (dFloat)(-M_GM / (fRadius*fRadius));
   dVector gravityForce  (mass * dGravity * toOrigin[0], mass * dGravity * toOrigin[1], mass * dGravity * toOrigin[2], 1.0f);

you can chose an earth center frame with the y axis aligned to the line the connect the north and south pole. That may be restrictive because it will limit your simulator to low earth orbit objects.

Right, except, my coordinate frame has a Z-axis from pole to pole, but that's trivial. Unfortunately, I HAVE to chose earth center frame because my visual display engine has an upper ceiling of 100,000,000 ft (30480 km) above ground, almost to an earth geostationary orbit - it will not allow positioning of objects any farther than that. The visual display engine was designed as a flight simulator, so the developer didn't see the need to have user go too far into space. I am OK with that, in fact, I have contacted the developer and they acknowledged that they might consider extending or removing this limit, but they want to see how far I can push this project first. So - Earth Center Frame it is, for now.

a more general inertial frame can be the sun a the origin and the up vector the perpendicular to the earth orbit around the sum, but in that model you can simulate the seasons

As far as seasons, they are simulated already, (different land textures and correct sun/star positions) but that is done based on time, not on simulating Earth as an object circling around sun.

For now, I am content with simulating Low to Medium Earth Orbit, with present-day and soon to be space hardware, (and of course, historic, Apollo, Gemini, Vostok). This is NOT a space shooter, but a realistic simulator I have given thought to simulating other celestial bodies (moon, Mars) by re-sizing and re-texturing the Earth into moon or Mars (I have an ability to do so using the SDK that is the part of the Flight Simulator), and letting user perform a trajectory burn towards, let's say, Mars, quitting Earth-mode, re-configuring for Mars-mode and re-starting in Mars-mode with a calculated incoming velocity trajectory from Earth. It would be like a hibernation, time-compression sleep

let us say you use the first coordinate system.
in this system your code that equation in your force and torque callback, the earth static is a static body with a fixed angular velocity that you integrate manually each frame. for now see if you can get model one working and them lots of things will make sense to you that do not seem logical now.

I totally agree. I will implement this shortly!

[Julio Jerez]

cool let us see how far you get with that.

a more general inertial frame can be the sun a the origin and the up vector the perpendicular to the earth orbit around the sum, but in that model you can simulate the seasons

As far as seasons, they are simulated already, (different land textures and correct sun/star positions) but that is done based on time, not on simulating Earth as an object circling around sun.

when I say simulate I mean the effect comes out of an emergent behavior form the equation for motion, yes you can simulate the season by implementing the equations, but with a sun base coordinate system, is different in that system earth become another body, but that's completelly a different simulation.

also with an earth bound system you can simulate the moon, and stuff way beyond the moon orbit, in fact the moon is kind of an exception as satellite, because is huge and yet is so closet to the earth that is in a orbital lock. The reason for that is that the earth act as if it was an elactic body because of the oceans. as the moon goes around the earth the earth deforms and that deformation steal energy from the equation -1.2 G m * m / r

That's another thing you could simulate is which the moon is in a lock orbit around the earth.
for that you only need another energy function that model the ocean deformation and have continues derivative. but that's is another problem altogether.

if you want to understand why I say this things even I have never done, it is because a fundament theorem of the calculus of variations called the Langragian of a system of particle.
you should check out some of the youtube videos
https://www.youtube.com/watch?v=zhk9xLjrmi4&t=64s

That's how I derive almost every single equation uses in the newton engine and how it possible in newton to make almost any joint that can be imagined.

and I do no care what the self appointed experts keep saying about newton, I know, that not withstanding a bugs, newton is correct. while when using impulse non of this stuff can de achieved.
If fact in newton I can even simulate electrical system that can be coupled with mechanical system, imagine that.

[misho]

Hi Julio,

Reading the equation:

Code: Select all

a = -g - 2 * corssProduct (EarthOmega, BodyLinalVelicty) - crossproduct (EarthOmega, crossProduct (EarthOmega, BodyPositionReletiveToErthCenter )

This is basically a summation of acceleration vectors, which multiplied by the mass of the object will finally give the total force acting on the object, correct?

In that equation, how exactly do I read: (since crossproduct result is a vector)

Code: Select all

2 * corssProduct (EarthOmega, BodyLinalVelicty)

Is it same as: (?)

Code: Select all

corssProduct (EarthOmega, BodyLinalVelicty).Scale(2.0)

Thanks,
Misho

[Julio Jerez]

yes that's right.

[misho]

thanks - OK, here are the results... my code:

Code: Select all

const double EARTH_ROT_INCR = 0.000072921150; // Earth Omega, in radians per sec.

NewtonBodyGetMatrix(body, &matrix[0][0]);
NewtonBodyGetMassMatrix(body, &mass, &Ixx, &Iyy, &Izz);
dVector vPosition = matrix.m_posit;
dVector vToOrigin = vPosition;
vNormalize(vToOrigin);

dVector vVelocityVector;
NewtonBodyGetVelocity(body, &vVelocityVector[0]);

double fRadius = sqrt(vPosition%vPosition);
dFloat dGravity = (dFloat)(-M_GM / (fRadius*fRadius));
dVector gravityAccel = vToOrigin.Scale(dGravity);

dVector EarthOmega = dVector(0.0f, 0.0f, EARTH_ROT_INCR);
dVector coriolisAccel = (EarthOmega*vVelocityVector).Scale(-2.0f) - EarthOmega*(EarthOmega*vPosition);

dVector TotalForce = (gravityAccel + coriolisAccel).Scale(mass);

NewtonBodySetForce(body, &TotalForce[0]);

With this in my ApplyForceAndTorqueCallback, the orbit starts decaying (the object starts falling) quite rapidly... from the initial orbit of 390km, it lost 100 km while rotating 65 degrees around the earth, so obviously, it would hit earth before the full orbit is completed.

If I use gravity only, as before, so:

Code: Select all

dVector TotalForce = gravityAccel.Scale(mass);

Everything goes back to normal, orbit is stable. Is there anything that I missed, just by looking at the code?

[Julio Jerez]

in general when a particle move sideway the in the direction of the rotation, the force will pool the body toward the center of rotation. the centripetal acceleration will push the body away.
Coriolis is much stronger that contritely and dominate the motion. this is what cause hurricanes and tornadoes.
here are some hints as to what direction the Coriolis force shoul point toward when the body moves.

The Coriolis effect is the behavior added by the Coriolis acceleration. The formula implies that the Coriolis acceleration is perpendicular both to the direction of the velocity of the moving mass and to the frame's rotation axis. So in particular:
if the velocity is parallel to the rotation axis, the Coriolis acceleration is zero.
if the velocity is straight inward to the axis, the acceleration is in the direction of local rotation.
if the velocity is straight outward from the axis, the acceleration is against the direction of local rotation.
if the velocity is in the direction of local rotation, the acceleration is outward from the axis.
if the velocity is against the direction of local rotation, the acceleration is inward to the axis.

I suggest you use one force at a time until you know all sign are right.
try -2 w cross Vel

and ignore gravity and w x (w x r) and check the direction, alone the angular velocity should pull you in opposite you push you away.
the add the centripetal, and centripetal away push you away.

you may have to use some tuning paramter damper Coriolis a bit. I would imagine the since omega is so small compare to everything else this forces should be very small compared to gravity.

earth omega = 7.2921159 × 10^−5 rad/sec
earth radius = 6.371 x 10^6 m

this mean that centrepeta sod be negligible as it is 7.2921159 × 10^−10 * 6.371 x 10^6

coriolisis should also be very small if the vehicle velocity is about a geo sync orbit speed of 3100 m/s
so Coriolis accel ~ 7.2921159 × 10^−5 * 3100 * 2 should be a small value compared to Gravity which show should be around 7 to 9 m/s2 a low orbital.

you figures should be around those values ballpark numbers.

[misho]

Ok - here are the readouts from the gravity and coriolis acceleration vectors, plus a velocity vector:
[EDIT: "CORL MAG is actually the sum of coriolis and centripetal accelerations]



For simplicity, the object is in equatorial orbit at 390 km (around the altitude of ISS), moving eastwards and is pointing towards Polaris (the sim has ability to draw constellations for educational purposes), and I paused the sim at the zero latitude and zero longitude position (just below west Africa). Positive X axis is coming out of East 90 longitude, and Positive Y axis is coming out of 0 longitude (where the object is in this snapshot). +ve Z axis is coming out of the north pole.

I nthis instance, gravity vector has a magnitude of 8.731, and is pointing along -ve Y axis, as expected, towards the centre of the earth

Coriolis vector has a magnitude of 1.084 (at a first glance, that seems to be too large compared to gravity magnitude), and is also pointing towards the centre of the earth, which seems not to be correct, from what you've quoted (Velocity and direction of rotation are in the same direction):

if the velocity is in the direction of local rotation, the acceleration is outward from the axis.

If I flip the sign of the coriolis vector, orbit altitude starts to increase instead of decreasing... so, that's not the solution either, because the magnitude of the coriolis vector seems too large compared to gravity vector. Any idea why the coriolis magnitude would be so large?

[misho]

Just as a follow up, for clarity, I've separated coriolis (CORL) and centripetal (CPTL) components:



You can see that the centripetal component is much smaller than the coriolis component, so the error is probably somewhere in the coriolis magnitude/direction calculation.

I'm calculating coriolis and centripetal components as:

Code: Select all

   dVector coriolisAccel = (EarthOmega*vVelocityVector).Scale(-2.0f);
   dVector centripetalAccel = EarthOmega*(EarthOmega*vPosition).Scale(-1.0f);

and the final force is calculated as:

Code: Select all

dVector TotalForce = (gravityAccel + coriolisAccel + centripetalAccel).Scale(mass);

[EDIT] Actually - when taking into account that my orbit at 390km is much lower than geosynchronous satellite and orbital velocity is around 7680 m/s, the ballpark values you were quoting in your last post are bang-on with the figures I'm getting in the screenshot above... Do you agree? In that case, I'm either not properly applying these forces, or the formula for coriolis component is incorrect...

[Julio Jerez]

Ok - here are the readouts from the gravity and coriolis acceleration vectors, plus a velocity vector:
[EDIT: "CORL MAG is actually the sum of coriolis and centripetal accelerations]

the Coriolis is the component CoriolisAccel = -2 * omega x velocity
of the three component this is the only force that is not aligned with the radius to the earth only the motion of the ship is perpendicular to the earth angular velocity and aligned parallel to the plane of rotation.

[Julio Jerez]

this values seems right, it look as if your ship is going too fast. you now nee to calculate what is the geo sync velocity that will keep everything in a steady.

here are to other test, you need to do.
as the ship fall toward the center does it should gain speed along the radius.
this new speed will create a new component that tangent to the linear velocity increasing the tangent speed of the ship until it find a steady this will increase the centripetal until the both cancel each other. or the ship goes into spiral orbit and crash.

the secund test is, change the direction of the velocity and see If the force is repulsive.

them you need a orbital computer function to calculate what should be the velocity of you ship to keel in orbit. is just use that function and set in one side and set the next accel to zero, and do rattson iterations as function of velocity,

seccalcial this si as fallow,
change velocity by a small delta velocity and calculate Accel, if acceleration goes up, then change velocity by small value on the opposite direction, is the acceleration keep don on that direction until the acceleration reach near zero or the point you need to keep in the orbit you want.
this method is very simple is called root by bisection, look it up in wiikipedia.

that will give a steady orbit velocity, that will last for a while, by you will nee to adjust periodically using the same function. the function is what control your thrusters form time to time.






if you invert the ship velocity does the void move away for the oprvit? second

[Julio Jerez]

oh there is a bug, sorry about that.
the velocity you plug in the coriolies term is not the total linear velocity of the ship relative to the inertia frame, it should be the linear velocity of the ship relative to the rotating medium.

for example say the you put a coin on a turn table, the coin will not experience any Coriolis forces because the coin is not moving in the turn table, it is only when the coin gain velocity that its espericence these side Coriolis force. If this is the case the velocity you use for Coriolis is

Code: Select all

V = shilpLinearvelocity - crossproduct (earthOmega, distanceToearthCenter);

now if a ship is in perfect orbit, its velocity should be.

Code: Select all

crossproduct (earthOmega, stanceToearthCenter)

and the centripal acceleration should cancel the gravitational pull because they are always aligned.
When move the ship start to move, it will need to change its velocity and this is when tiny coriolis forces act on the ship because the term
(shilpLinearvelocity - crossproduct (earthOmega, distanceToearthCenter)
is no longer zero but a small value that generate a non zero force.
please try that.