|
-
January 4th, 2011, 08:21 PM
#1
Calculating planetary coordinates [Physics algorithm]
I have a very quick question, and I figured because this is an algorithmic question, perhaps someone would be able to answer this for me here.
I've done a small algorithm in c++ to calculate the coordinates of planets based on the force applied to them, which gives the acceleration, velocity, and finally position (all in components)
I'll explain my issue first and attach my code at the bottom.
I'm not sure if this is an issue or not, but the sun is moving (in meters) quite far away from the origin. I'm aware that the sun is supposed to move a small distance due to the pull of the larger planets, but the coordinates are reaching into the tens of thousands of meters - which does not seem to be logical to me.
The processing for my algorithm is as follows:
-Reset force to 0;
-Get the sum of all forces by F += GM(x[j]-x[i])/r^3
-calculate specific planet's acceleration by the equation a = F/m
-calculate specific planet's velocity by V = Vinitial + a*t
-calculate specific planet's position by X = Xinitial + v*t
-loop
To me, this sounds like this should make sense - and so I figured it was an algorithmic error.
Anyone who's good in physics or algorithm have any idea?
Below is the actual C++ code:
Code:
double Fx[10], Fy[10], Fz[10]; // Array of force vectors
float gravity_constant = 8.99E-9; // Kepler's gravity constant
// Set the forces of all the planets to 0
for( int i = 0; i < 10 ; i++)
{
Fx[i] = 0;
Fy[i] = 0;
Fz[i] = 0;
}
for( int i = 0; i < 10 ; i++)
{
for( int j = 0; i < 10 ; i++)
{
if(j!=i) // if not itself (added to ensure no division by zero error)
{
// calculate the force at x
Fx[i] += (objects[j].x - objects[i].x) / sqrt(
pow((objects[j].x - objects[i].x),2) +
pow((objects[j].y - objects[i].y),2) +
pow((objects[j].z - objects[i].z),2))
* (gravity_constant * objects[i].mass * objects[j].mass)
/ (
pow((objects[j].x - objects[i].x),2) +
pow((objects[j].y - objects[i].y),2) +
pow((objects[j].z - objects[i].z),2));
// calculate the force at y
Fy[i] += (objects[j].y - objects[i].y) / sqrt(
pow((objects[j].x - objects[i].x),2) +
pow((objects[j].y - objects[i].y),2) +
pow((objects[j].z - objects[i].z),2))
* (gravity_constant * objects[i].mass * objects[j].mass)
/ (
pow((objects[j].x - objects[i].x),2) +
pow((objects[j].y - objects[i].y),2) +
pow((objects[j].z - objects[i].z),2));
// calculate the force at z
Fz[i] += (objects[j].z - objects[i].z) / sqrt(
pow((objects[j].x - objects[i].x),2) +
pow((objects[j].y - objects[i].y),2) +
pow((objects[j].z - objects[i].z),2))
* (gravity_constant * objects[i].mass * objects[j].mass)
/ (
pow((objects[j].x - objects[i].x),2) +
pow((objects[j].y - objects[i].y),2) +
pow((objects[j].z - objects[i].z),2));
}
}
}
time += interval; // Increase time by the interval
// Loop for calculating acceleration, velocity, and position
for( int i = 0; i < 10 ; i++)
{
// calculate acceleration
objects[i].ddx = Fx[i] / objects[i].mass;
objects[i].ddy = Fy[i] / objects[i].mass;
objects[i].ddz = Fz[i] / objects[i].mass;
// calculate velocity
objects[i].dx += objects[i].ddx * interval;
objects[i].dy += objects[i].ddy * interval;
objects[i].dz += objects[i].ddz * interval;
// calculate position
objects[i].x += objects[i].dx * interval;
objects[i].y += objects[i].dy * interval;
objects[i].z += objects[i].dz * interval;
// Output the objects
drawPlanet(i, objects[i].x, objects[i].y, objects[i].z);
}
where objects is a struct:
Code:
struct tagObjects
{
double x, dx, ddx; // x component of position, velocity, acceleration
double y, dy, ddy; // y component of position, velocity, acceleration
double z, dz, ddz; // z componont of position, velocity, acceleration
double mass; // mass of the object (planet)
double radius; // radius of the object (for deciding width of spheres)
}objects[11];
and initialized with starting positions/values
Code:
objects[0].mass = 1.99E+30;
objects[0].radius = 109;
objects[0].texture = NULL;
objects[0].x =
objects[0].y =
objects[0].z = 0;
objects[0].dx = 9.303626;
objects[0].dy =-11.739118;
objects[0].dz =-5.243946;
// mercury
objects[1].mass = 3.34E+23;
objects[1].radius = 0.382512033;
objects[1].texture = NULL;
objects[1].x =-20665696392;
objects[1].y =-59636889090;
objects[1].z =-29712844059;
objects[1].dx = 3.66E+04;
objects[1].dy =-9.51E+03;
objects[1].dz =-8.88E+03;
// venus
objects[2].mass = 4.87E+24;
objects[2].radius = 0.948840564;
objects[2].texture = NULL;
objects[2].x =-1.07478E+11;
objects[2].y =-6401037913;
objects[2].z = 3920831085;
objects[2].dx = 8.82E+02;
objects[2].dy =-3.19E+04;
objects[2].dz =-1.45E+04;
// earth
objects[3].mass = 5.98E+24;
objects[3].radius = 1;
objects[3].texture = NULL;
objects[3].x =-26516914541;
objects[3].y = 1.32754E+11;
objects[3].z = 57555479554;
objects[3].dx =-2.98E+04;
objects[3].dy =-4.78E+03;
objects[3].dz =-2.06E+03;
// mars
objects[4].mass = 6.40e23;
objects[4].radius = 0.532603753;
objects[4].texture = NULL;
objects[4].x = 2.08092E+11;
objects[4].y = 1150108018;
objects[4].z =-5098849339;
objects[4].dx = 1.30E+03;
objects[4].dy = 2.39E+04;
objects[4].dz = 1.09E+04;
// jupiter
objects[5].mass = 1.90E+27;
objects[5].radius = 11.2089807;
objects[5].texture = NULL;
objects[5].x = 5.94749E+11;
objects[5].y = 4.1426E+11;
objects[5].z = 1.63068E+11;
objects[5].dx =-7.89E+03;
objects[5].dy = 1.02E+04;
objects[5].dz = 4.54E+03;
// saturn
objects[6].mass = 5.69E+26;
objects[6].radius = 9.44920901;
objects[6].texture = NULL;
objects[6].x = 9.48999E+11;
objects[6].y = 9.31359E+11;
objects[6].z = 3.4387E+11;
objects[6].dx =-7.44E+03;
objects[6].dy = 6.12E+03;
objects[6].dz = 2.85E+03;
// uranus
objects[7].mass = 8.67E+25;
objects[7].radius = 4.00730625;
objects[7].texture = NULL;
objects[7].x = 2.17596E+12;
objects[7].y =-1.85516E+12;
objects[7].z =-8.43327E+11;
objects[7].dx = 4.66E+03;
objects[7].dy = 4.29E+03;
objects[7].dz = 1.81E+03;
// neptune
objects[8].mass = 1.03E+26;
objects[8].radius = 3.88266098;
objects[8].texture = NULL;
objects[8].x = 2.5472E+12;
objects[8].y =-3.41681E+12;
objects[8].z =-1.46188E+12;
objects[8].dx = 4.50E+03;
objects[8].dy = 2.91E+03;
objects[8].dz = 1.08E+03;
// pluto
objects[9].mass = 1.03E+26;
objects[9].radius = 0.185008075;
objects[9].texture = NULL;
objects[9].x =-1.42071E+12;
objects[9].y =-4.20637E+12;
objects[9].z =-8.84208E+11;
objects[9].dx = 5.28E+03;
objects[9].dy =-1.96E+03;
objects[9].dz =-2.19E+03;
Please, if anyone has any ideas of what I'm doing wrong - assuming these values are wrong as I think they are - please send reply and tell me what I am doing wrong
Tags for this Thread
Posting Permissions
- You may not post new threads
- You may not post replies
- You may not post attachments
- You may not edit your posts
-
Forum Rules
|
Click Here to Expand Forum to Full Width
|