Thanks,JimCalculating forces between stars in 3D space

[bkelly]

Hello,
What is the best way to determine direction, force, and velocity in three dimensional Cartesian coordinates?

Explanation: I am writing a program to do some star motion simulations. I will edit in some stars, give them their mass, along with initial position and velocity, and see what happens. My plan is to track the position and the velocity with the X, Y, and Z Cartesian values. The distance is simple but I am uncertain as to calculating the force and resultant accelerations and velocity in 3 D.

Will it work our right if I calculate each of the X, Y, and Z values separate from each other? By that I mean work each axis one at a time, do the force, acceleration, and velocity calculations in each axis, and ignore the other two? That sounds like it would work but seems too simplistic.

 

[Werg22]

I take it that you the said forces and accelerations are due to the interaction (gravity) of the stars among themselves. If that's the case, given a star A and the entire collection of the remaining stars, B, C, D..., the resultant force on A at some instant is the vector sum of the forces of gravity of B, C, D, ... on A. However, to accurately describe the motion of the system is a near impossible task (see http://en.wikipedia.org/wiki/N-body_problem#Three-body_problem").

What I would suggest is treat time in quanta. Say you divide time into quanta of 1 second, then starting at t = 0 with some initial values for velocity and acceleration for A, at t = 1 A will have a velocity and position calculated by taking the acceleration as constant (and in the same direction) through t = 0 and t = 1. Doing this for the remaining stars, the acceleration is recalculated at t=1 according to their new positions, and so on. The motion of the system resembles more and more the continuous time ideal as you chose smaller divisions.

 

[bkelly]

You are correct. I did not want to go into that aspect is it gets complicated quickly. In essence, there will be a linked list of stars, call them A,B,C, and D. Start with A. It compares itself with B, calculates the force imparted by B and direction of force and saves them. Then a compares itself with C, calaculates that force and direction, sums it with the force from B, and does the sam with D. Then B compares itself with A, C, D, etc. After making a round robin of all the comparisons, I calculate, say, one hours worth of time for velocity and position, then update all the velocities and positions. Don't take that all that literally, just conceptually. It does get complicated.

The question is, can I calculate each of X, Y, and Z independantly of each other?

Example, if the distance between two stars is 4, 5, and 6 hundred million meters along the x, y, and z axes respectfully, the straighline line distance between the two is about 8.7. But if I calculate the distance and resultant force for 4 hundred million meters along the x axis, and do the same for y and z independantly, will that work out. Logically it seems like it would work, but I'm a math lightweight and doubt that it will work in practice.

I have found and worked through the equations for straight line force, and am now trying to determine how they will work in 3D.

How would you work this arithmetic?

 

[bkelly]

BTW: The little 3D graphics in the page you linked to is exactly what I want to do. I will be able to edit in the mass, position, and velocity of each star. There is more, but let's keep the side effects short.

 

[Werg22]

I see. Yes it's possible to get X, Y and Z independently. The magnitude of the force along each axis is taken as being proportional to the component of the distance separating the two objects along that same axis. For example, say star A is centered at (0, 0 , 0) and star B at (3, 4, 12). Then the components of the distance separating them are 3, 4 and 12. The force due to gravity along the X axis is equal to 3k, and along the Y and Z axis 4k and 12k respectively, for some constant k. Then k is calculated by referring to the equality ((3k)^2 + (4k)^2 + (5k)^2)^1/2 = GP/169, or, simplified, 13k = GP/169 (where P stands for the product of the stars masses). In general, k = GP/R^3.

Hope that helps.

 

[bkelly]

Application

What does your k stand for?

How do I apply this to my calculations?

We can work out the earth’s orbit around the sun in two dimensions. Start at t+0 with the sun at the center, x = 0 and y = 0, and the Earth at about 149,476,000 kilometers from the sun, or in meters, 149*e9 meters. The earth’s velocity along the x-axis is about 53,000 kph, or 14,872 m/sec along the x-axis (positive for simplicity )and at this instance, 0 along the y axis. (For simplicity, let's calculate as though the Earth has no effect upon the sun and the Earth is in the quadrent such that motion along the y-axis toward the sun is expressed as negative values.)

My calculations tell me that the gravitation force between the sun and the Earth is 38e21 Newtons. Divide that by the mass of the Earth and the acceleration on the Earth is 0.00622 meters/sec. (Hmm, that seems to work out to 22.392 meters per second per hour. If you want to see my calcuations I can post, but you probably already know all this.)

Note: I am using Microsoft Excel for quick calculations. Plugging in the earth’s mass, radius, and 1 kg weight at the surface, I get 9.8 m/sec^2 acceleration. That looks pretty good. Substituting in the mass and distances for Earth and sun, I get the values I just mentioned.)

How do I apply the acceleration of 0.00622 meters/sec^2 to the earth’s motion to get the correct values for velocity along x and y axes at t+3600 ( After 1 hour or 3600 seconds of time? I expect the velocity along the x-axis to be reduced some small amount and the velocity along the y-axis to move in the negative direction some amount.

 

[Werg22]

Given that v^2/r = GP/r^2, v^2 = GP/r. Here v stands for the velocity (which is constant) of the Earth along its circular orbit. The angular velocity of the earth, say w, is v/r. After time 3600, assuming you're taking the angle to be 0 at t = 0, the angle shifts to 3600*v/r. To get the velocity along the X and Y axis at this time, you calculate the horizontal and vertical components of the vector with magnitude v and angle 3600*v/r. That would be v*cos(3600*v/r) and v*sin(3600*v/r).

On k: Given two objects A and B, say you want to calculate the force on A imparted by B along each axes. Say A is positioned at (a, b, c) and B at (d, e, f). Then the force along the X axis is k*(d - a), along the Y axis k*(e - b) and along the Z axis k*(f - c). Note that this gives the magnitude and the direction of the force (you can get positive and negative values depending on the position of B relative to A).

About the circular orbit: A circular orbit would imply that the acceleration is constantly changing (not in magnitude, but in direction). Isn't this something that won't be possible in your program? If you decide to update the velocities, positions and accelerations at intervals of time, rather than going about a circular orbit, the objects will fall onto each other in a motion describing a spiral made up of line segments (the spiral could have many loop though, so it might not be that they'll fall onto each other before a while).

 

[bkelly]

Calculating forces

Hello Werg,

Regarding the circular orbit: I think that if two stars are positioned such that a circular orbit should ensue, the calculations should work out such that it does. That will be my first test. Put a star about the size of the sun at coordinates 0,0,0 and a much smaller one about the size of the Earth at 150,000,000 kilometers, give the small one the right velocity, and see if it orbits in something close to a circular orbit.

After reading your post, I have these procedures to calculate the force between the stars and to determine the forces along each of the x, y, and z axis. Do you think these will work?

For each pair of stars, calculate the distance between them with a three dimensional Pythagoream as follows:

void CStar::Calculate_Distance_Tther_Star()
{
delta_x_position = x_position - other_star_x_position;
delta_y_position = y_position - other_star_y_position;
delta_z_position = z_position - other_star_z_position;
double temp = ( delta_x_position * delta_x_position ) +
( delta_y_position * delta_y_position ) +
( delta_z_position * delta_z_position );
distance_to_other_star = sqrt( temp );
};

Contrary to many “C” programmers, I like to use several variables in simple statements rather than complex statements. It doesn’t slow the process down one bit and it provides better insight into the code when using the debugger.

The total force between the two stars is calculated here.

void CStar::Calculate_Force_Towards_Other_Star()
{
double temp_1 = mass / distance_to_other_star;
temp_1 *= Gravitational_Constant;
double temp_2 = other_star_mass / distance_to_other_star;
attraction_between_stars = temp_1 * temp_2;
};
These are huge numbers so calculations are broken up so as to reduce the magnitude of the values thereby reducing floating point errors.

Then for each of the three axes, use the cosine value of adjacent over hypotenuse times the force to give the force along that axis. The angle used has its apex at the other star and the opposite is the difference between the two stars along the X axis. When delta_x goes towards zero, the force along X goes to zero.

sine = delta_x_position / distance_to_other_star;
x_axis_force = sine * attraction_between_stars;
x_axis_force -= x_axis_force;

The same equation will work for Y and Z axis.

That’s probably too much for one post so I stop here. Does this look like it will work so far.