CodeGuru Home VC++ / MFC / C++ .NET / C# Visual Basic VB Forums Developer.com
Results 1 to 5 of 5
  1. #1
    Join Date
    Jan 2011
    Posts
    9

    Exclamation 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

  2. #2
    Join Date
    Jan 2011
    Posts
    9

    Re: Calculating planetary coordinates [Physics algorithm]

    I just realized 2 mistakes...
    1) I defined G as being 8.99E-9, when it should have been 6.67E-11, and
    2) I have pluto and neptune's mass equal - when pluto should be 6.58E+23

    still though, this doesn't seem to affect overly much... The sun is seemingly moving without any return.

  3. #3
    Join Date
    Jan 2011
    Posts
    9

    Re: Calculating planetary coordinates [Physics algorithm]

    Gah.. **** anti post-editing.
    I found another error, and yet it still doesn't affect my output overly much.

    I realize in the loop I put i++ on the j loop, and I fixed that small error (because before it was not calculating the force on the sun), and yet, the sun is still moving to distances over 10000m away from the origin.
    Is this normal for the sun?

  4. #4
    Join Date
    Jan 2011
    Posts
    9

    Re: Calculating planetary coordinates [Physics algorithm]

    ...If a moderator could please close/delete this thread, I'd appreciate it.
    I can't believe I just posted 4 times because I couldn't edit, only to figure out that the sun's initial velocity was what was messing everything up in the end. o_o'

  5. #5
    Join Date
    Jan 2006
    Location
    Fox Lake, IL
    Posts
    15,007

    Re: Calculating planetary coordinates [Physics algorithm]

    Post the solution, for others, use THREAD TOOLS at the top of this post to CLOSE THE THREAD YOURSELF.

    And, go back and add
    Code:
    // Code Tags like this
    David

    CodeGuru Article: Bound Controls are Evil-VB6
    2013 Samples: MS CODE Samples

    CodeGuru Reviewer
    2006 Dell CSP
    2006, 2007 & 2008 MVP Visual Basic
    If your question has been answered satisfactorily, and it has been helpful, then, please, Rate this Post!

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
  •  





Click Here to Expand Forum to Full Width

Featured