Star Grav. Forces

[Shadow]

Hello all,

 

I'll get straight to the point. I've been trying to write a program that would simulate the movement of stars, and I've been having slight problems in thinking up the most fundamental part of this idea: the way gravitation interacts with those stars. Here's some basic info on what I need, how I need it, and how I've been trying to solve the problem.

 

Every star (I'll be calling them objects) has it's given mass, position, and velocity. We are in Cartesian space. Since I'm using C as a programing language, I have no choice but to represent the velocities and their directions as X and Y speeds. You probably have no idea what I mean, so let me explain on an example.

 

 

Here we have and object situated at the coordinates [math] x = 4[/math] and [math] y = -3[/math]. As you can see, the blue line is the velocity vector of that object. In plain language, it say "In one unit of time, this object will be displaced by this much in this direction". However, I had to be able to represent BOTH the size and the direction of the vector using nothing but dimensions. I can't use angles, or anything like that. The only way I could think of doing that was...well... the way I did it, since I can't find the right words to describe it. But that's the point of an image, though, isn't it? The direction of the X and Y vector can be expressed by it's positivity or negativity.

 

So, the vector [math]\vec{v}[/math] can be represented as [math]\vec{v} = \sqrt{x^2 + y^2} \hspace {2mm} where \hspace {2mm} \vec{x} = 3, \hspace {2mm} \vec{y} = -2[/math]

 

However, the squaring of X and Y would destroy our sense of direction, since squared numbers are always positive. This would mean that we'd have to store the direction [math](\pm 1)[/math] in a separate variable, and multiply the result by this variable in the end. But there is no need for that, since I am automatically representing everything in the "X and Y" way, so there is really nothing to be squared. To sum it up, in plain speak I am saying "This object will move by this much in the X direction (left, right, left is negative, right is positive) and by this much in the Y direction (up, down, up is positive, down is negative) in one unit of time".

 

 

Now that you know this part, I can go on to explain the actual problem I'm encountering. Say we have two objects, a and b, with the following definitions:

 

object[a].coordX = 2

object[a].coordY = 1

object[a].velX = 2

object[a].velY = 3

object[a].mass = 5 kg

 

object.coordX = 3

object.coordY = 3

object.velX = -2

object.velY = 1

object.mass = 7 kg

 

Here is an image representing the objects and their velocities.

 

 

Now, here is my problem; I need to calculate the change in the X/Y velocities of both object due to gravity. Here is how I TRY do it.

 

 

 

The effects of gravity on object A

 

diffX = object.coordX - object[a].coordX = 1 -> The distance gravity has to "travel" on the X plane

 

diffY = object.coordY - object[a].coordY = 2 -> The distance gravity has to "travel" on the Y plane

velX = object[a].velX = 2 -> The X velocity of the object

velY = object[a].velY = 3 -> The Y velocity of the object

r -> distance between the two objects

[math]F_g[/math] -> gravitational force

x -> a temporary variable

gravX -> The X value of gravity

gravY -> The Y value of gravity

finVelX -> The final amount of units by which the object will be displaced on the X plane

finVelY -> The final amount of units by which the object will be displaced on the Y plane

 

 

[math]r = \sqrt{diffX^2 + diffY^2}[/math]

 

[math]r = \sqrt{5}[/math]

 

 

[math]F_g = G \frac{m_1m_2}{r^2}[/math]

 

[math]F_g = 46.69 \cdot 10^{-11}[/math]

 

 

[math]x = \frac{F_g}{\sqrt{diffX^2+diffY^2}}[/math]

 

[math]x = 9.2 \cdot 10^{-11} \cdot \sqrt{5}[/math]

 

 

[math]gravX = diffX \cdot x = 9.2 \cdot 10^{-11} \cdot \sqrt{5} [/math]

 

[math]gravY = diffY \cdot x = 18.4 \cdot 10^{-11} \cdot \sqrt{5}[/math]

 

 

[math]finVelX = gravX + velX = 2 + (9.2 \cdot 10^{-11} \cdot \sqrt{5})[/math]

[math]finVelY = gravY + velY = 3 + (18.4 \cdot 10^{-11} \cdot \sqrt{5})[/math]

 

 

The effects of gravity on object B

 

 

 

diffX = object[a].coordX - object.coordX = -1 -> The distance gravity has to "travel" on the X plane

 

diffY = object[a].coordY - object.coordY = -2 -> The distance gravity has to "travel" on the Y plane

velX = object.velX = -2 -> The X velocity of the object

velY = object.velY = 1 -> The Y velocity of the object

r -> distance between the two objects

[math]F_g[/math] -> gravitational force

x -> a temporary variable

gravX -> The X value of gravity

gravY -> The Y value of gravity

finVelX -> The final amount of units by which the object will be displaced on the X plane

finVelY -> The final amount of units by which the object will be displaced on the Y plane

 

 

[math]r = \sqrt{diffX^2 + diffY^2}[/math]

 

[math]r = \sqrt{5}[/math]

 

 

[math]F_g = G \frac{m_1m_2}{r^2}[/math]

 

[math]F_g = 46.69 \cdot 10^{-11}[/math]

 

 

[math]x = \frac{F_g}{\sqrt{diffX^2+diffY^2}}[/math]

 

[math]x = 9.2 \cdot 10^{-11} \cdot \sqrt{5}[/math]

 

 

[math]gravX = diffX \cdot x = -9.2 \cdot 10^{-11} \cdot \sqrt{5} [/math]

 

[math]gravY = diffY \cdot x = -18.4 \cdot 10^{-11} \cdot \sqrt{5}[/math]

 

 

[math]finVelX = gravX + velX = -2 + (-9.2 \cdot 10^{-11} \cdot \sqrt{5})[/math]

[math]finVelY = gravY + velY = 1 + (-18.4 \cdot 10^{-11} \cdot \sqrt{5})[/math]

 

This would mean that in the time between [math]t_0[/math] and [math]t_1[/math], object[a] will be displaced by [math]2 + (9.2 \cdot 10^{-11} \cdot \sqrt{5})[/math] to the right, and by [math]3 + (18.4 \cdot 10^{-11} \cdot \sqrt{5})[/math] up, while object will be displaced by [math]2 + (9.2 \cdot 10^{-11} \cdot \sqrt{5})[/math] to the left, and by [math]1 - (18.4 \cdot 10^{-11} \cdot \sqrt{5})[/math] up. To calculate the amount object[a/b] will be displaced between [math]t_1[/math] and [math]t_2[/math] I would have to redo all these calculations using the exact same values except for r.

 

Could someone please verify this and tell me if it works, and, if not, propose a better (working ) solution?

 

I thank anyone who even bothers to read this in advance, and may God bless those who actually understand it

 

Shadow

[h4tt3n]

You're on the right track! I have to go now, but I'll post a solid reply to your questions later tonight.

 

Cheers, Michael

[D H]

So, the vector [math]\vec{v}[/math] can be represented as [math]\vec{v} = \sqrt{x^2 + y^2} \hspace {2mm} where \hspace {2mm} \vec{x} = 3, \hspace {2mm} \vec{y} = -2[/math]

 

Whoa! That thing [math]\vec{v}[/math] is not a vector. You do not need to square things here. Much better is to truly use vectors: [math]\vec v = 3\hat x -2 \hat y[/math].

 

object[a].coordX = 2

object[a].coordY = 1

object[a].velX = 2

object[a].velY = 3

object[a].mass = 5 kg

 

object.coordX = 3

object.coordY = 3

object.velX = -2

object.velY = 1

object.mass = 7 kg

Better; now you are using vectors -- but you are missing something. What are your units? You have units for the masses, but not for the positions or for the velocities.

 

Now, here is my problem; I need to calculate the change in the X/Y velocities of both object due to gravity. Here is how I TRY do it.

 

Much better is to use the vector form for gravitational acceleration,

 

[math]\vec a_{1,2}

= \frac{Gm_2}{||\vec r_{1,2}||^3} \vec r_{1,2}[/math]

 

where [math]\vec a_{1,2}[/math] is the acceleration of object 1 resulting from the gravitational force exerted on object 1 by object 2 and [math]\vec r_{1,2}[/math] is the displacement vector from object 1 to object 2.

 

In terms of your example, assuming you are using SI units (in which case [math]G = 6.673 \cdot 10^{-11} \, \text{m}^3 / \text{kg} / \text{s}^2[/math])

 

[math]\vec r_{a,b}

= (\hat x + 2 \hat y)\,\text{meters}[/math]

 

[math]\vec a_{a,b}

= (4.178 \cdot 10^{-11}\hat x + 8.356 \cdot 10^{-11} \hat y) \,\text{meters}/\text{second}^2[/math]

 

[math]\vec r_{b,a} = -\vec r_{a,b}

= -(\hat x + 2 \hat y)\,\text{meters}[/math]

 

[math]\vec a_{b,a}

= -(2.984 \cdot 10^{-11}\hat x + 5.969 \cdot 10^{-11} \hat y) \,\text{meters}/\text{second}^2[/math]

 

Note that I was very explicit regarding units. Doing so helps catch errors.

 

Now for the issue of propagating the state over time. This problem is an example of an initial value problem -- a set of ordinary differential equations that describe the evolution of the state of some system coupled with a complete description of the state at some initial point. Analytic solutions (e.g., Kepler's equations) exist for a system of two masses. For multiple bodies one must resort to numerical techniques that yield approximations of the true solution. There are many ways to solve such problems numerically. You are using something called the Euler method. This is a good start, but beware that the Euler method almost always yields very lousy results, and this is definitely the case for central force problems. Once you get the Euler integration working you should move on to a better integrator. The most commonly used technique is fourth-order Runge-Kutta integration. You might also want to look into the slightly simpler verlet integrator, velocity verlet, or Heun's method.

[h4tt3n]

Ok,

 

As you've probably noticed by now, there isn't much programming related help to get on this forum. I recommend you to take a look at these physics programming forums instead:

 

http://www.gamedev.net/community/forums/

 

http://www.bulletphysics.com/Bullet/phpBB3/

 

Regarding your post: Your proposed method does not work, and this is because you do not seem to understand the concepts of movement in classical mechanics.

 

First, you need to know Newtons three laws of motion, here in the quck and dirty version:

 

(1) Every object stays in the same state of motion (or at rest) unless an external force is applied to it.

 

(2) When a force is applied to an object it accellerates according to this rule: Force = Mass * Acceleration

or

Acceleration = Force / Mass

 

(3) For every action (Force) there is an opposite and equal reaction (also a force).

 

When calculating how two or more bodies influence each other gravitationally, you need to update these values every time step:

 

Force -> acceleration -> velocity -> position

 

All of these are vectors, and I recommend that you define an x and y value for each in order to keep track of them.

 

You've got the calculation of force right, so lets move on to acceleration. When two bodies attract each other they are influenced by the same force, except in opposite direction (Newton's 3rd law). But the resulting acceleration may be different!

 

The acceleration of M1 and M2 due to gravitational force F is:

 

A(M1) = F / M1

 

A(M2) = F / M2

 

So, if M1 is twice as big as M2 then M1 is only accelerated half as much as M2.

 

Now comes the tricky part... we need to figure out the new velocity and position of the celestial bodies as a function of the gravitational force applied on them. This is called integrating with respect to time, and we do this with a so called numerical integration algorithm.

 

There are many to choose from, but for a simple simulation like yours I recommend the symplectic Euler integration method:

 

New velocity = old velocity + acceleration * timestep

New position = old position + new velocity * timestep

 

It is very important to understand tha you can calculate gravitational force and acceleration accurately, but that you can't do the same with velociy and position. It is still an unsolved matemathical problem to accurately predict the motion of 3 or more bodies due to graviational interaction. The different integration algorithms around try to "predict" movement as accurately as possible, but they can never be exact. You can increase accuracy by making the timestep value smaller.

 

(Once you get more advanced you'll need to implement the 4th order Runge-Kutta algorithm!)

 

I'll write some more on how to apply all this in a program later...

 

Cheers, Michael

[Shadow]

Hey mike,

 

I've long since found out that the calculations I posted above are rubbish While I didn't know half the things you wrote about in your previous post, I did manage to make the program work, I did it I guess half my way half your way. You mentioned, in the other post, that you'd give me some constructive feedback on my source. That'd be brilliant !! If you could send me you're ICQ/Skype or something by PM, I'd love to go over the source with you. It's exactly this kind of help that I need the most and that I have none of )

 

Let me know and thank you again,

 

Shadow

[h4tt3n]

You're welcome. IMHO this is exactly what forums like this are for. I'll be away for the weekend, but I'll pm you my email sometime sunday. Anyway, for further reference, why not post the code here? That way our discussion may help others in a similar situation...

 

cheers, Michael

[Shadow]

The code is pretty long, and also pretty incomplete I think it'd be more productive if I post the code after it's finished and...well...presentable :D but thanks, I'll be waiting for your address