# Why does my binary star simulation only work for equal masses and initial speeds?

I made a 2D gravity simulation is JavaScript using something called p5.js, you can find it here: https://editor.p5js.org/christofferaakre/sketches/ZVfm6cPR

if i place two objects with certain masses into the simulation with no initial speeds, then as expected they accelerate towards each other. Something weird happens at the end, but that is just because i haven't told the simulation to anything if they get too close, so as the distance between the two objects get very small the accelerations blow up. When I give one of the objects an initial velocity perpendiculor to the distance between them, I get an elliptical orbit as expected. However, I am having trouble getting binary stars to work. Most of the time, I get something that looks two elliptical orbits about a common centre of mass, as expected, but each consecutive orbit is shifted down by a constant amount. See the image below.

I got the above image when i tried two objects with different masses and equal but opposite velocities. At first I though maybe this could be explained by the binary star system 'moving through space', since the relative positions of the stars are indeed two ellipcital orbits, but relative to the space around them they are also moving downwards.

However, I think there is some other issue with my code, because if I only give one of the stars an initial velocity, this happens:

In this picture, only the object represented by the green trajectory had an initial velocity. Somehow the fact that it had no initial velocity means it only does half an ellipse?? Indeed, if I give if even just a small velocity in the opposite direction, it does complete an orbit:

I'm hoping someone here has a clue as to what is going on here.

• What numerical scheme are you using to solve the equations? – tpg2114 Apr 20 '20 at 13:05
• I think it's called Euler Integration. I'm just calculating the force of gravity between the two stars, then using that to update the velocities, and using the velocities to update the positions every time step – Christoffer Corfield Aakre Apr 20 '20 at 13:06
• I think you'll find that Euler integration will introduce other, similar problems down the road. It does not conserve energy and so you'll see your orbits wander instead of form closed loops (when they are supposed to form closed loops, that is). You'll want to use what's called a "symplectic integrator" or "symplectic scheme." An intro to what that means is provided on our sister site, SciComp.SE. – tpg2114 Apr 20 '20 at 13:36
• @tpg2114 thanks, but I actually knew this already. I just threw this together quickly for fun, not doing it for anything serious. Thanks though :) – Christoffer Corfield Aakre Apr 20 '20 at 13:37

It simply looks like your initial condition causes the system to move through space. You need to work in the centre of mass frame. Set $$m_1v_1=-m_2v_2$$

• Thanks, this makes sense to me. However, I don't think it explains why i only got half an orbit when i gave only one of the stars an initial velocity, or does it? – Christoffer Corfield Aakre Apr 20 '20 at 13:14
• Actually nevermind my first comment, I dont think it makes much sense haha, I will accept your answer. – Christoffer Corfield Aakre Apr 20 '20 at 13:15
• Both stars need an initial velocity so that the center of mass is stationary. – R.W. Bird Apr 20 '20 at 16:25
• In your example, set moon's vy as '-earth.vy * earth.m / 1000' where 1000 is the moon's m – Jason Goemaat Apr 20 '20 at 21:52

The fact you get good behaviour with equal masses might suggest you are implementing the two-body solution incorrectly.

You can find the common centre of mass by weight ratios but the "reduced mass" for calculating the imaginary body at this point (to make this into 2 single body equations) is not as simple as just using the reduced mass equation found on wiki. (I found this rather strange)

I use the following for the doing two body calculations (in javascript)

    var m1 = body1.mass;
var m2 = body2.mass;
var m3 = m1+m2;
var mr = Math.pow(m2,3)/Math.pow(m3,2); //reduced mass


Then the gravity acceleration is defined by the reduced mass (mr) and the distance from the centre of mass (for body1)

    var g = GRAVITATIONAL_CONSTANT * mr / Math.pow( r, 2 );


If you don't use the above two body solution your calculations are affected by the time based change in the gravity field at any given point, the reduced mass equations give you a frame and method to describe the gravitational field as fixed and with an inverse square relationship from the centre of mass. When the masses are equal, it is likely your 'reduced mass' method coincides with this (I had this experience)

Yes, this is a reduced mass equation that looks different to the text-book definition, but I found that this is what worked for me when I was doing two body modelling on the path to multiple body models. It looks strange, but sketch it on paper and derive it yourself - it does work!

I am not saying the text-books are wrong, the usual equation is symmetric for both bodies - it is just it is very unclear how it is meant to be implemented.