This is a page about »Simulating Solar System Part 2«.

Extending to N-Body Gravitational Force

Since the Earth Moon system is working, its time to extend to the enture solar system. The biggest change here is the gravitational force calculation. It is now,

And given acceleration is \(a_j=\frac{F_{ij}}{m_j}\) the calculation can be simplified to,

Hence,

 1void calculateAcceleration(std::vector<CelestialObj> celestialBodies) {
 2  // acceleration is now `CelestialObj`'s attribute
 3  acceleration = glm::vec2{0.f, 0.f};
 4  for (auto& c : celestialBodies) {
 5    if (c.name != this->name) {
 6      glm::vec2 direction = glm::vec2{c.x, c.y} - glm::vec2{this->x, this->y};
 7      float r = glm::length(direction);  // sqrt((x^2) + (y^2))
 8      glm::vec2 unitDirection = direction / r;
 9      acceleration += c.mass / (r * r) * unitDirection;
10    }
11  }
12  acceleration *= GravitationalConstant;
13}

However the change result in an unexpected behavior - both objects remains still. It took days before the bug is identified.

C++ std::vector

The celestial bodies are created and put into std::vector as described below.

 1CelestialObj earth{"earth", 0., 0., EarthRadius, EarthMass, sf::Color(159, 193, 100), false};
 2CelestialObj moon{"moon", 3.633e8, 0., MoonRadius, MoonMass, sf::Color(201, 201, 201), true};
 3std::vector<CelestialObj> solarSystem = {earth, moon};
 4
 5int main() {
 6  // ...
 7  moon.vy = 1.076e3;
 8  earth.vy = -(moon.vy * moon.mass) / earth.mass;
 9  //...
10  for (auto& c: solarSystem) {
11    c.calculateAcceleration(solarSystem);
12  }
13  for (auto& c: solarSystem) {
14    c.updatePos(c.acceleration, false);
15  }
16  window.clear(sf::Color::Black);
17  earth.render(window);
18  moon.render(window);
19}

Coming from Python, adding earth and moon to a list is as straightforward as described above. However in c++, the solarSystem actually makes a copy of earth and moon. The assignment of *.vy applies to the original copy of earth and moon, yet the calculateAcceleration and updatePos interacts with the copies. Finally the render renders the original copies.

the realization took too long as the author focus was on the calculateAcceleration correctness and the for loops (separated calculation and update position or merged)

Yet Another Bug

One of the fix attempt described below leads to an interesting outcome where earth is not rendered but the moon and its trail are.

 1void calculateAcceleration(std::vector<CelestialObj> celestialBodies) {
 2  this->acceleration.x = 0.0f;
 3  this->acceleration.y = 0.0f;
 4  for (CelestialObj c : celestialBodies) {
 5    if (c.name != this->name) {
 6      glm::vec2 direction = glm::vec2{c.x, c.y} - glm::vec2{this->x, this->y};
 7      float r = glm::length(direction);
 8      glm::vec2 unitDirection = direction / r;
 9      this->acceleration += c.mass / (r * r) * unitDirection;
10    }
11  }
12  this->acceleration *= GravitationalConstant;
13}
14
15std::vector<CelestialObj> solarSystem = {earth, moon};
16
17int main() {
18  moon.vy = 1.076e3;
19  earth.vy = -(moon.vy * moon.mass) / earth.mass;
20  earth.calculateAcceleration(solarSystem);
21  moon.calculateAcceleration(solarSystem);
22  earth.updatePos(earth.acceleration, false);
23  moon.updatePos(moon.acceleration, false);
24  earth.render(window);
25  moon.render(window);
26}

One of the main differences is between for (auto& c: solarSystem) {} and for (CelestialObj c : celestialBodies) {}. The prior with auto refers to objects in solarSystem through pointers reference i.e. any changes to c is inplace, while the latter first creates a copy of the CelestialObj i.e. the objects in solarSystem are immutable (at the cost of performance).

When slowing down the simulation to 1fps, earth start moving to the right in a curves with gradual increase of step size. At each tick, a copy of moon is created at the position 3.633e8m away from earth. The orginal copy of earth has a positive vy explains why it moves in a U shape. And as the distance between the original earth and the moon copy gets smaller, the gravitational force and acceleration increased exponentially due to \(\frac{1}{r^2}\). If the view width permits, one should observe earth moves diagonally down till the infinite in evenly spaced distance. This is due to the upward force is cancelled by the moon’s attraction to a right-downward direction force, and as earth move away from moon the gravitational force have less influence to earth’s trajectory and from Newton’s first law, it will move in a “straight” line at “constant” speed. The moon’s trajectory is correct because a copy of earth is created at fixed x/y position.

code is found on branch. Shift the centering code setPosition for better visibility together m2vp to scale the view width/height to see the trajectory. Another thing to try is to change the initial speed of earth.vy to observe different trajectory. Since a copy of moon/earth at fixed position is created throughout the iterations, the value earth.vy matters less in terms conservation of momentum.

Correct Approach to std::vector