# Orbit Simulator pseudocode

1. Jun 27, 2008

### 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?

Last edited: Jun 27, 2008
2. Jun 27, 2008

### 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 wont be a perfect circle (it will be a bit elliptic), but it will not become a spiral...

3. Jun 27, 2008

### DaveC426913

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.

4. Jun 27, 2008

### D H

Staff Emeritus
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:
PHP:

# 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];
}
}
}

5. Jun 27, 2008

### DaveC426913

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

6. Jun 27, 2008

### D H

Staff Emeritus
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 d1 and d2 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.

7. Jun 27, 2008

### DaveC426913

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.

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.

8. Jun 27, 2008

### D H

Staff Emeritus
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!

Last edited: Jun 27, 2008
9. Jun 27, 2008

### 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. =/

10. Jul 3, 2008

### DaveC426913

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

11. Jul 3, 2008

### DaveC426913

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).

12. Jul 3, 2008

### JaWiB

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

13. Jul 3, 2008

### DaveC426913

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.

14. Jul 3, 2008

### JaWiB

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

15. Jul 3, 2008

### 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.

16. Jul 3, 2008

### DaveC426913

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.

cuz that was a zillion years ago.

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

Though I'm open to any methods that work.

17. Jul 3, 2008

### D H

Staff Emeritus
Sure you have, in post #4:

18. Jul 3, 2008

### DaveC426913

omg, that seems to be working!

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

19. Jul 3, 2008

### D H

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

20. Jul 3, 2008

### DaveC426913

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.