Unable to get a proper orbit on my solar system model

Question

I am making a solar system model using just p5.js and without box2D though the language/ platform would not matter for this question. Also, all the numbers and variables used to describe the problem are not 100% accurate, but the behavior is the same.

I used the Newton formula ( F = G Mm / r^2 )for finding mutual gravity between two objects, say A and B. Now to attract A towards B, I am dividing this mutual gravity by A's mass to find the centripetal acceleration on A and then multiplying this by a unit vector pointing towards B. On applying this relation to both A and B, they both experience attraction towards each other, inversely proportional to their mass.

Now if i leave both of them to interact with each other with A having mass = 1000 units and B having mass = 10 units, as expected A does not get pulled and remains stationary but B gets accelerated towards A. Now what happens is as B reaches the center and flies of in the opposite direction, it goes farther then the distance I placed it at initially. This keeps on snowballing during each acceleration cycle and at some point, goes off the screen. This seems to as a violation of the conservation of energy or some major flaw in my physics.

Alright then moving on to the second problem, we have the same objects and masses. The difference is that i give B (the lighter object) an initial velocity of some value, say x in the positive direction of x axis. Now i place B perpendicular to A's x axis and let them interact. this time, B moves in a sort of elliptical orbit with two issues. first is the A (heavier object) is not at one of the focii of the ellipse, and instead is at the exact center of it and second is that over time, the orbit itself starts rotating. I feel that this rotation is caused by the initial velocity provided and just to make it clear, the velocity is applied initially only and is not applied on each frame. The traced path of this orbit is as follows:

Also notice how the maximum extent of each orbit is a bit more then the previous. This is pretty much the previous problem mixed in with this.

The next thing which I am currently trying is, applying a constant velocity tangential to motion along with the gravital centripetal acceleration. Let me know if this would be useful or if my whole approach needs to be changed.

Also this my code for the simulation :

var constG;
var axisX;
var planets = [];

function setup() {
  createCanvas(500, 500);
  //createCanvas(displayWidth, displayHeight);
  //fullscreen(true);
  constG = 0.0001;//6.67 * pow(10, -11);
  axisX = createVector(1, 0);
}

function draw() {
  background(0, 5);
  for (var planet of planets) {   
    planet.update();
    planet.display();
  }

  for (var i = 0; i < planets.length; i++){
    var selfPlanet = planets[i];
    for (var j = 0; j < planets.length; j++){
      if (j == i){
        continue;
      }
      var otherPlanet = planets[j];
      var gravitalAcc = calcGravitalAcc(selfPlanet, otherPlanet);
      selfPlanet.applyForce(gravitalAcc);
    }
  }
  
  if (planets.length > 0){
    planets[0].radius = 15;
    planets[0].mass = 100;    // this just makes the first planet heavy so that i 
    planets[0].vel.mult(0);   // can test stuff while making it the sun.
    planets[0].speed = 0;
  }
}

function mousePressed() {
    planets.push(new CelestialBody(mouseX, mouseY, 7));
}

function calcGravitalAcc(self, other){
  var tempVec = p5.Vector.sub(other.pos, self.pos);
  return tempVec.normalize().mult(constG * (other.mass)/pow(tempVec.mag(), 2))
}

and this is the Celestialbody class, just any typical class in simple physics simulations:

class CelestialBody {
  constructor(x, y, radius) {
    this.pos = createVector(x, y);
    this.radius = radius;
    this.color = color(255);
    
    this.mass = 1;
    this.speed = 1;
    this.vel = createVector(1, 0) //p5.Vector.random2D();
    this.vel.setMag(this.speed);
    this.acc = createVector(0, 0);
  }

  display() {
    fill(this.color);
    stroke(this.color);
    circle(this.pos.x, this.pos.y, this.radius * 2);
  }

  update() {
    this.pos.add(this.vel);
    this.vel.add(this.acc);
    this.acc.mult(0);
  }
  
  applyForce(vForce){
    this.acc.add(vForce);
  }
}

Answer 1

the first problem is usually caused by too big time step of the simulation in conjuction with absent collision handling. when your objects get close the forces get big and the incremental step in simulation became too big so the next iterated position is after the collision and usually further away then before it so the breaking force is smaller leading to bigger and bigger orbits in time...

To remedy that you can do:

• add collision handling
• limit max speed or force to some sane limit

I never encountered your second problem and without code I can not even guess ... other than rounding errors

Take a look at this:

• s it possible to make realistic n-body solar system simulation in matter of size and mass?

an also all sublinks in there ...

[Edit1] after I saw your code

The architecture of the code looks OK the problem is that your equations are slightly off. it should be:

vel+=acc*dt;
pos+=vel*dt;
acc=0.0;

instead of yours:

pos+=vel;
vel+=acc;
acc=0.0;

so you got wrong order and missing the *dt where dt is the time step. Because of this no matter how you change your timer interval the result is the same (just slower/faster) and also the direction of acceleration is applied one step latter than should causing the spinning of orbit (because acceleration was computed from different position than it was applied to the final position so its direction is always off).