How can I improve my Orbit Simulation logic?

[DaveC426913]

Orbit Simulation logic

I'm building an Orbit Simulator (this is what I build every time I learn a new language). I'm looking for the pseudocode for calculating the acceleration.

I know the formula: G*m/(d2-d1)^2, I'm trying to figure out why I'd need the hypotenuse.

You see, I can calculate the accelX and accelY just from the x and y coords. Then I can udpate the x and y velocities, then I can update the x and y positions. Nowhere do I actually need to calculate the hypotenuse.

But this yields incorrect motion.

So, am I supposed to
1] calculate the hypotenuse from the x/y coords,
2] then figure out the acceleration as a single scalar value along the hypotenuse
3] then immediately reconvert the single scalar acceleration back into accX/accY?

 

[Zizy]

Could you please tell which language are you learning now and in which languages did you manage to write sim this way? Cause if I remember correctly, I was warned specifically for orbit sim that trajectory will be a spiral cause of finite precision...
If you set energy to be a constant and calculate from there, trajectory won't be a perfect circle (it will be a bit elliptic), but it will not become a spiral...

 

[DaveC426913]

Could you please tell which language are you learning now and in which languages did you manage to write sim this way? Cause if I remember correctly, I was warned specifically for orbit sim that trajectory will be a spiral cause of finite precision...
If you set energy to be a constant and calculate from there, trajectory won't be a perfect circle (it will be a bit elliptic), but it will not become a spiral...

I've written it in classic VB, Java and now I'm writing it in Flash/AS3, but frankly, it doesn't matter which language, it is simply the logic (i.e. pseudocode) that I'm looking for.

 

[D H]

The expression in the first post gives the magnitude but not the direction of the acceleration. I much prefer the vectorial form of Newton's law of gravitation, in which the acceleration on body j due to body i is

$$\mathbf a_{j \to i} = \frac {Gm_i}{||\mathbf r_{j \to i}||^3}\mathbf r_{j \to i}$$

where $\mathbf r_{2\to1}$ is the vector from body j due to body i:

$$\mathbf r_{j \to i} = \mathbf r_1 - \mathbf r_2$$

Pseudo code:

# Compute force of body j to body i as output
Vect force_j_to_i ( # Returns force vector
   double mu_i,   # Mass of body i
   Vect r_i,      # Position of body i
   double mu_j,   # Mass of body j
   Vect r_j)      # Position of body j
{
   Vect r_j_to_i = r_i - r_j;
   double rsq = r_j_to_i.dot(r_j_to_i);
   double ffact = mu_i*mu_j/(rsq*sqrt(rsq));
   return r_j_to_i * afact;
}
# Compute accelerations of all bodies toward each other
void compute_accels (
   Vect * a,    # Output accelerations
   double * m,  # Input positions
   Vect * r,    # Input positions
   int n)       # Number of bodies
{
   for (int ii=0; ii<n; ii++) {
      a[ii] = 0;
   }
   for (int ii=0; ii<n-1; ii++) {
      for (int jj=ii+1; jj<n; jj++) {
         double force_j_i = force_j_to_i (m[ii], r[ii], m[jj], r[jj]);
         a[jj] += force_j_i/m[jj];
         a[ii] -= force_j_i/m[ii];
      }
   }
}

 

[DaveC426913]

The expression in the first post gives the magnitude but not the direction of the acceleration.

It gives a magntitude in the x direction and a magnitude in the y direction.

 

[D H]

Ohh. Well that's wrong then. You will get wrong answers if you compute the terms individually. The acceleration isn't even going in the correct direction!

You can use the expression in post #1, but you need to use that expression to compute the magnitude of the acceleration, with d₁ and d₂ being vectors, not components. Then compute the unit vector from body 2 to body 1. Finally, scale that unit vector by the acceleration magnitude to yield the acceleration vector.

or

You can use the vectorial form of Newton's law of gravitation that I posted. It cuts out a lot of intermediate calculations.

 

[DaveC426913]

Ohh. Well that's wrong then. You will get wrong answers if you compute the terms individually. The acceleration isn't even going in the correct direction!

You can use the expression in post #1, but you need to use that expression to compute the magnitude of the acceleration, with d₁ and d₂ being vectors, not components.

Yeah, this is what my intuition is telling me. But when I map out the logic, I don't seem to get it. Other than slope (i.e. rise/run i.e. x2-x1 / y2-y1) how else would one compute the accel.? I suppose I could use degrees but I don't recall using that.

You can use the vectorial form of Newton's law of gravitation that I posted. It cuts out a lot of intermediate calculations.

Hm. Yeah. I wonder if that code could be any more hard-to-read? My variables are: thisStar for the outer loop and otherStar for the inner loop, so I always know what I'm computing.

Also, I've never heard of Vect as a data type.

 

[D H]

That's partly the physicist in me, and partly the quality of the pseudo code you get at PF consulting rates.

Physicists use short, terse names because that makes complex equations easier to read. I've had this battle with software engineering weenies, and I have always prevailed in this battle because

• I use long meaningful names where that makes sense,
• I use short meaningful names where that makes more sense,
• For example, I find v_cbody_to_tbody_in_cbody_from_clvlh, or even more tersely, v_cb2tb__cb_clvlh, much easier to read than velocity_from_chaser_body_to_target_body_expressed_in_chaser_body_coordinates_but_viewed_from_chaser_local_vertical_local_horizontal_frame. Long names get in the way when you have a lot of different quantities, different reference points, different reference frames crawling all over. It's even worse with velocities because the frame in which they are observed and the frame in which they are expressed can differ.
• The code is easy to validate via inspection and review because the code looks a whole lot like the equations in the detailed architecture document,
• The design is easy to validate via inspection and review because the equations in the design document look a whole lot like the equations in the referenced texts, journal papers, etc.

I even have a nifty little perl script to generate snippets of C code from my LaTeX equations!

 

[WarPhalange]

I tried making one in C# last summer. It sucked. I wanted it to have graphics, though. I don't know if you're doing that. The graphics killed me. Too many things to learn too fast. =/

 

[DaveC426913]

The expression in the first post gives the magnitude but not the direction of the acceleration. I much prefer the vectorial form of Newton's law of gravitation, in which the acceleration on body j due to body i is

$$\mathbf a_{j \to i} = \frac {Gm_i}{||\mathbf r_{j \to i}||^3}\mathbf r_{j \to i}$$

where $\mathbf r_{2\to1}$ is the vector from body j due to body i:

$$\mathbf r_{j \to i} = \mathbf r_1 - \mathbf r_2$$

Pseudo code:

# Compute force of body j to body i as output
Vect force_j_to_i ( # Returns force vector
   double mu_i,   # Mass of body i
   Vect r_i,      # Position of body i
   double mu_j,   # Mass of body j
   Vect r_j)      # Position of body j
{
   Vect r_j_to_i = r_i - r_j;
   double rsq = r_j_to_i.dot(r_j_to_i);
   double ffact = mu_i*mu_j/(rsq*sqrt(rsq));
   return r_j_to_i * afact;
}
# Compute accelerations of all bodies toward each other
void compute_accels (
   Vect * a,    # Output accelerations
   double * m,  # Input positions
   Vect * r,    # Input positions
   int n)       # Number of bodies
{
   for (int ii=0; ii<n; ii++) {
      a[ii] = 0;
   }
   for (int ii=0; ii<n-1; ii++) {
      for (int jj=ii+1; jj<n; jj++) {
         double force_j_i = force_j_to_i (m[ii], r[ii], m[jj], r[jj]);
         a[jj] += force_j_i/m[jj];
         a[ii] -= force_j_i/m[ii];
      }
   }
}

What is this 'Vect' data type? I don't know how to handle that in AS3.

 

[DaveC426913]

You can use the expression in post #1, but you need to use that expression to compute the magnitude of the acceleration, with d₁ and d₂ being vectors, not components.

I'm afraid I'm stuck here. How do I quantify a direction and a magnitude without quantifying them as x and y components? (I feel silly, I've done this sim several times before).

 

[JaWiB]

I'm afraid I'm stuck here. How do I quantify a direction and a magnitude without quantifying them as x and y components? (I feel silly, I've done this sim several times before).

Your vector contains both x and y components, but then you subtract one position vector from another and take the square of the magnitude of the result.

If whatever language you're using doesn't have vector objects/operations available, then you'll have to explicitly perform operations on the components

 

[DaveC426913]

Your vector contains both x and y components, but then you subtract one position vector from another and take the square of the magnitude of the result.

If whatever language you're using doesn't have vector objects/operations available, then you'll have to explicitly perform operations on the components

That's fine. So I simply keep my x and y coords, as in:

starArray[thisStar].dX, starArray[thisStar].dY
starArray[otherStar].dX, starArray[otherStar].dY

but it's not right.

I'm sure it's not right because I never need to calc the hypotenuse - and I'm sure I should.

 

[JaWiB]

Yep. That's the same as finding the magnitude of the resulting vector in the equation.

 

[Ben Niehoff]

What language are you using? If it's object-oriented, you may wish to define some kind of vector object.

Letting one object be at the origin, Newton's law of gravitation is

$$\vec F = - G m_1 m_2 \frac{\vec r}{|\vec r|^3}$$

You do need to calculate the hypotenuse in your code, because you need the hypotenuse cubed in the denominator. E.g., the x component of the force would be

$$F_x = - G m_1 m_2 \frac {x} {(x^2 + y^2)^{3/2}}$$

I'm still a bit confused, though, because if you've done this a million times before, can't you just look up what you did in your previous code?

Also, as someone else noted, the method you are using to calculate orbits (i.e., finding the acceleration, using it to change the velocity, and using that to change the position) is highly susceptible to rounding errors.

 

[DaveC426913]

What language are you using? If it's object-oriented, you may wish to define some kind of vector object.

Letting one object be at the origin, Newton's law of gravitation is

$$\vec F = - G m_1 m_2 \frac{\vec r}{|\vec r|^3}$$

You do need to calculate the hypotenuse in your code, because you need the hypotenuse cubed in the denominator. E.g., the x component of the force would be

$$F_x = - G m_1 m_2 \frac {x} {(x^2 + y^2)^{3/2}}$$

Ah. Uh. I've never seen these equations before. I'm pretty sure I didn' do it this way last time, but I'm game.

I'm still a bit confused, though, because if you've done this a million times before, can't you just look up what you did in your previous code?

cuz that was a zillion years ago.

Also, as someone else noted, the method you are using to calculate orbits (i.e., finding the acceleration, using it to change the velocity, and using that to change the position) is highly susceptible to rounding errors.

Yeah, that's OK, it's just a toy.

Though I'm open to any methods that work.

 

[D H]

What language are you using? If it's object-oriented, you may wish to define some kind of vector object.

Letting one object be at the origin, Newton's law of gravitation is

$$\vec F = - G m_1 m_2 \frac{\vec r}{|\vec r|^3}$$

You do need to calculate the hypotenuse in your code, because you need the hypotenuse cubed in the denominator. E.g., the x component of the force would be

$$F_x = - G m_1 m_2 \frac {x} {(x^2 + y^2)^{3/2}}$$

Ah. Uh. I've never seen these equations before. I'm pretty sure I didn' do it this way last time, but I'm game.

Sure you have, in post #4:

The expression in the first post gives the magnitude but not the direction of the acceleration. I much prefer the vectorial form of Newton's law of gravitation, in which the acceleration on body j due to body i is

$$\mathbf a_{j \to i} = \frac {Gm_i}{||\mathbf r_{j \to i}||^3}\mathbf r_{j \to i}$$

 

[DaveC426913]

omg, that seems to be working!

Except my formula is
$$F_x = - G m_2 \frac {x} {(x^2 + y^2)^{3/2}}$$

 

[D H]

Ben posted a formula for force. Yours looks like a formula for acceleration.

 

[DaveC426913]

Ben posted a formula for force. Yours looks like a formula for acceleration.

Right. That part I remember. You end up dividing the force by the mass of the star, to arrive at the accel. which is why the m cancels out.

It is working in that I now have an arbitrary number of stars in orbit around each other (and one of those stars even has a planet! That's hard to do at this scale!)

I do see some precession - which I recall is an artifact of a rough algorithm that needs tweaking.

 

[dwahler]

It sounds like you're using Euler integration, which is simple but tends to accumulate errors pretty quickly. You can improve the precision a lot by switching to another method with better error bounds, like Runge-Kutta.

 

[DaveC426913]

It sounds like you're using Euler integration, which is simple but tends to accumulate errors pretty quickly. You can improve the precision a lot by switching to another method with better error bounds, like Runge-Kutta.

As much of a physics buff as I am, I nearly failed Calculus in HS. This is all geek to me.

 

[DaveC426913]

Definitely something still fundamentally wrong with my sim algorithm. The orbits are more like spirograph patterns - this is worse than simple precession of the orbit.

Can anyone see where I've misunderstood?

for (thisSt=0; thisSt<numStars; thisSt++){
	var accX = 0;
	var accY = 0;
	var distX = 0;
	var distY = 0;
	for (otherSt=0; otherSt<numStars; otherSt++){
		if (otherSt != thisSt){
[B]			distX = 0 + ArrSt[otherSt].dX - ArrSt[thisSt].dX;
			distY = 0 + ArrSt[otherSt].dY - ArrSt[thisSt].dY;
			var hyp = calcHypot(ArrSt[thisSt].dX,ArrSt[thisSt].dY,ArrSt[otherSt].dX,ArrSt[otherSt].dY);		
			accX += 0 + G *  ArrSt[otherSt].m * distX / Math.pow(hyp,3);
			accY += 0 + G *  ArrSt[otherSt].m * distY / Math.pow(hyp,3);[/B]
			ArrSt[thisSt].vX += accX;
			ArrSt[thisSt].vY += accY;
		}
	}
	ArrSt[thisSt].dX += ArrSt[thisSt].vX;
	ArrSt[thisSt].dY += ArrSt[thisSt].vY;
	//move symbol on stage
	ArrSt[thisSt].x = ArrSt[thisSt].dX - ArrSt[thisSt].width/2;
	ArrSt[thisSt].y = ArrSt[thisSt].dY - ArrSt[thisSt].height/2;
}

private function calcHypot(thisX,thisY,otherX,otherY){
	var dist = Math.sqrt(Math.abs(thisX-otherX) +  Math.abs(thisY-otherY));
	return dist;
}

 

[D H]

The immediate problem is that you are computing the hypotenuse incorrectly. You should be computing

$$d = \sqrt{x^2+y^2}$$

There are several other problems as well:

1. You have a fixed time step (even worse, an implicit time step) of one second (or whatever time units are in your version of G.)
2. You are intermingling the calculation of the accelerations with the state integration. This precludes you from using any but the simplest of integration techniques.
3. You are using the simplist of integration techniques.
4. You have a fixed display scale factor (even worse, an implicit scale factor) of one pixel per unit length.

 

[DaveC426913]

The immediate problem is that you are computing the hypotenuse incorrectly. You should be computing

$$d = \sqrt{x^2+y^2}$$

Doh! I was so busy suspecting the acceleration formula it never occurred to me to examine the basic hypotenuse formula.
Thanks! That's a much more well-behaved orbit - just a bit of a decay problem now.

As for 1 and 4, now that I've got the basic algorithm nailed, I can begin abstracting the variables such as t and G.

As for 2 and 3, I'm afraid that: a new language + still unsteady with OOP + only high school maths - may leave more complex algorithms a bit beyond me without some serious hand-holding.