Vector direction for gravity in a circular orbit

Question

I am currently working on a project in C# where i play around with planetary gravitation, which i know is a hardcore topic to graps to it's fullest but i like challenges. I've been reading up on Newtons laws and Keplers Laws, but one thing i cannot figure out is how to get the correct gravitational direction.

In my example i only have 2 bodies. A Satellite and a Planet. This is to make is simplify it, so i can grasp it - but my plan is to have multiple objects that dynamically effect each other, and hopefully end up with a somewhat realistic multi-body system.

When you have an orbit, then the satellite has a gravitational force, and that is ofcourse in the direction of the planet, but that direction isn't a constant. To explain my problem better i'll try using an example:

let's say we have a satellite moving at a speed of 50 m/s and accelerates towards the planet at a speed of 10 m/s/s, in a radius of 100 m. (all theoretical numbers) If we then say that the framerate is at 1, then after one second the object will be 50 units forward and 10 units down.

As the satellite moves multiple units in a frame and about 50% of the radius, the gravitational direcion have shifted alot, during this frame, but the applied force have only been "downwards". this creates a big margin of error, especially if the object is moving a big percentage of the radius.

In our example we'd probably needed our graviational direction to be based upon the average between our current position and the position at the end of this frame.

How would one go about calculating this?

I have a basis understanding of trigonometry, but mainly with focus on triangles. Assume i am stupid, because compared to any of you, i probably am.

(I made a previous question but ended up deleting it as it created some hostility and was basicly not that well phrased, and was ALL to general - it wasn't really a specific question. i hope this is better. if not, then please inform me, i am here to learn :) )

Just for reference, this is the function i have right now for movement:

foreach (ExtTerBody OtherObject in UniverseController.CurrentUniverse.ExterTerBodies.Where(x => x != this))
{
    double massOther = OtherObject.Mass;

    double R = Vector2Math.Distance(Position, OtherObject.Position);

    double V = (massOther) / Math.Pow(R,2) * UniverseController.DeltaTime;

    Vector2 NonNormTwo = (OtherObject.Position - Position).Normalized() * V;

    Vector2 NonNormDir = Velocity + NonNormTwo;
    Velocity = NonNormDir;

    Position += Velocity * Time.DeltaTime;
}

If i have phrased myself badly, please ask me to rephrase parts - English isn't my native language, and specific subjects can be hard to phrase, when you don't know the correct technical terms. :)

I have a hunch that this is covered in keplers second law, but if it is, then i'm not sure how to use it, as i don't understand his laws to the fullest.

Thank you for your time - it means alot!

(also if anyone see multi mistakes in my function, then please point them out!)

Answer 1

I am currently working on a project in C# where i play around with planetary gravitation

This is a fun way to learn simulation techniques, programming and physics at the same time.

One thing I cannot figure out is how to get the correct gravitational direction.

I assume that you are not trying to simulate relativistic gravitation. The Earth isn't in orbit around the Sun, the Earth is in orbit around where the sun was eight minutes ago. Correcting for the fact that gravitation is not instantaneous can be difficult. (UPDATE: According to commentary this is incorrect. What do I know; I stopped taking physics after second year Newtonian dynamics and have only the vaguest understanding of tensor calculus.)

You'll do best at this early stage to assume that the gravitational force is instantaneous and that planets are points with all their mass at the center. The gravitational force vector is a straight line from one point to another.

Let's say we have a satellite moving at a speed of 50 m/s ... If we then say that the framerate is one frame per second then after one second the object will be 50 units right and 10 units down.

Let's make that more clear. Force is equal to mass times acceleration. You work out the force between the bodies. You know their masses, so you now know the acceleration of each body. Each body has a position and a velocity. The acceleration changes the velocity. The velocity changes the position. So if the particle starts off having a velocity of 50 m/s to the left and 0 m/s down, and then you apply a force that accelerates it by 10 m/s/s down, then we can work out the change to the velocity, and then the change to the position. As you note, at the end of that second the position and the velocity will have both changed by a huge amount compared to their existing magnitudes.

As the satellite moves multiple units in a frame and about 50% of the radius, the gravitational direcion have shifted alot, during this frame, but the applied force have only been "downwards". this creates a big margin of error, especially if the object is moving a big percentage of the radius.

Correct. The problem is that the frame rate is enormously too low to correctly model the interaction you're describing. You need to be running the simulation so that you're looking at tenths, hundredths or thousanths of seconds if the objects are changing direction that rapidly. The size of the time step is usually called the "delta t" of the simulation, and yours is way too large.

For planetary bodies, what you're doing now is like trying to model the earth by simulating its position every few months and assuming it moves in a straight line in the meanwhile. You need to actually simulate its position every few minutes, not every few months.

In our example we'd probably needed our graviational direction to be based upon the average between our current position and the position at the end of this frame.

You could do that but it would be easier to simply decrease the "delta t" for the computation. Then the difference between the directions at the beginning and the end of the frame is much smaller.

Once you've got that sorted out then there are more techniques you can use. For example, you could detect when the position changes too much between frames and go back and redo the computations with a smaller time step. If the positions change hardly at all then increase the time step.

Once you've got that sorted, there are lots of more advanced techniques you can use in physics simulations, but I would start by getting basic time stepping really solid first. The more advanced techniques are essentially variations on your idea of "do a smarter interpolation of the change over the time step" -- you are on the right track here, but you should walk before you run.

Answer 2

I'll start with a technique that is almost as simple as the Euler-Cromer integration you've been using but is markedly more accurate. This is the leapfrog technique. The idea is very simple: position and velocity are kept at half time steps from one another.

The initial state has position and velocity at time t₀. To get that half step offset, you'll need a special case for the very first step, where velocity is advanced half a time step using the acceleration at the start of the interval and then position is advanced by a full step. After this first time special case, the code works just like your Euler-Cromer integrator.

In pseudo code, the algorithm looks like

void calculate_accel (orbiting_body_collection, central_body) {
    foreach (orbiting_body : orbiting_body_collection) {
        delta_pos = central_body.pos - orbiting_body.pos;
        orbiting_body.acc =
            (central_body.mu / pow(delta_pos.magnitude(),3)) * delta_pos;
    }
}

void leapfrog_step (orbiting_body_collection, central_body, delta_t) {
    static bool initialized = false;
    calculate_accel (orbiting_body_collection, central_body);
    if (! initialized) {
        initialized = true;
        foreach orbiting_body {
            orbiting_body.vel += orbiting_body.acc*delta_t/2.0;
            orbiting_body.pos += orbiting_body.vel*delta_t;
        }
    }
    else {
        foreach orbiting_body {
            orbiting_body.vel += orbiting_body.acc*delta_t;
            orbiting_body.pos += orbiting_body.vel*delta_t;
        }
    }
}

Note that I've added acceleration as a field of each orbiting body. This was a temporary step to keep the algorithm similar to yours. Note also that I moved the calculation of acceleration to it's own separate function. That is not a temporary step. It is the first essential step to advancing to even more advanced integration techniques.

The next essential step is to undo that temporary addition of the acceleration. The accelerations properly belong to the integrator, not the body. On the other hand, the calculation of accelerations belongs to the problem space, not the integrator. You might want to add relativistic corrections, or solar radiation pressure, or planet to planet gravitational interactions. The integrator should be unaware of what goes into those accelerations are calculated. The function calculate_accels is a black box called by the integrator.

Different integrators have very different concepts of when accelerations need to be calculated. Some store a history of recent accelerations, some need an additional workspace to compute an average acceleration of some sort. Some do the same with velocities (keep a history, have some velocity workspace). Some more advanced integration techniques use a number of techniques internally, switching from one to another to provide the best balance between accuracy and CPU usage. If you want to simulate the solar system, you need an extremely accurate integrator. (And you need to move far, far away from floats. Even doubles aren't good enough for a high precision solar system integration. With floats, there's not much point going past RK4, and maybe not even leapfrog.)

Properly separating what belongs to whom (the integrator versus the problem space) makes it possible to refine the problem domain (add relativity, etc.) and makes it possible to easily switch integration techniques so you can evaluate one technique versus another.

Answer 3

So i found a solution, it might not be the smartest, but it works, and it's pretty came to mind after reading both Eric's answer and also reading the comment made by marcus, you could say that it's a combination of the two:

This is the new code:

foreach (ExtTerBody OtherObject in UniverseController.CurrentUniverse.ExterTerBodies.Where(x =>  x != this))
{
double massOther = OtherObject.Mass;

double R = Vector2Math.Distance(Position, OtherObject.Position);

double V = (massOther) / Math.Pow(R,2) * Time.DeltaTime;

float VRmod = (float)Math.Round(V/(R*0.001), 0, MidpointRounding.AwayFromZero);
if(V > R*0.01f)
{
    for (int x = 0; x < VRmod; x++)
    {
        EulerMovement(OtherObject, Time.DeltaTime / VRmod);
    }
}
else
    EulerMovement(OtherObject, Time.DeltaTime);

}

public void EulerMovement(ExtTerBody OtherObject, float deltaTime)
    {

            double massOther = OtherObject.Mass;

            double R = Vector2Math.Distance(Position, OtherObject.Position);

            double V = (massOther) / Math.Pow(R, 2) * deltaTime;

            Vector2 NonNormTwo = (OtherObject.Position - Position).Normalized() * V;

            Vector2 NonNormDir = Velocity + NonNormTwo;
            Velocity = NonNormDir;



            //Debug.WriteLine("Velocity=" + Velocity);
            Position += Velocity * deltaTime;
    }

To explain it:

I came to the conclusion that if the problem was that the satellite had too much velocity in one frame, then why not seperate it into multiple frames? So this is what "it" does now.

When the velocity of the satellite is more than 1% of the current radius, it seperates the calculation into multiple bites, making it more precise.. This will ofcourse lower the framerate when working with high velocities, but it's okay with a project like this.

Different solutions are still very welcome. I might tweak the trigger-amounts, but the most important thing is that it works, then i can worry about making it more smooth!

Thank's everybody that took a look, and everyone who helped be find the conclusion myself! :) It's awesome that people can help like this!