n-body problem parallel code in c++ example 1

Piotr Borowiecki
#include <cmath>
#include <iostream>
#include <vector>
#include <omp.h>

const double G = 6.674e-11;  // gravitational constant

struct Body {
  double x, y, z;  // position
  double vx, vy, vz;  // velocity
  double mass;  // mass
};

// Calculate the gravitational force between two bodies.
void gravitationalForce(Body &body1, Body &body2, double &fx, double &fy, double &fz) {
  double dx = body2.x - body1.x;
  double dy = body2.y - body1.y;
  double dz = body2.z - body1.z;
  double distSqr = dx * dx + dy * dy + dz * dz;
  double f = G * body1.mass * body2.mass / distSqr;
  double fx = f * dx / distSqr;
  double fy = f * dy / distSqr;
  double fz = f * dz / distSqr;
}

int main() {
  std::vector<Body> bodies;
  // Initialize the array of bodies with their positions, velocities, and masses.

  double dt = 0.01;  // time step
  for (int t = 0; t < 100; ++t) {  // simulate for 100 time steps
    // Use OpenMP to parallelize the loop over all pairs of bodies.
    #pragma omp parallel for
    for (int i = 0; i < bodies.size(); ++i) {
      for (int j = i + 1; j < bodies.size(); ++j) {
        Body &body1 = bodies[i];
        Body &body2 = bodies[j];
        double fx, fy, fz;
        gravitationalForce(body1, body2, fx, fy, fz);
        body1.vx += dt * fx / body1.mass;
        body1.vy += dt * fy / body1.mass;
        body1.vz += dt * fz / body1.mass;
        body2.vx -= dt * fx / body2.mass;
        body2.vy -= dt * fy / body2.mass;
        body2.vz -= dt * fz / body2.mass;
      }
    }
    
    // Update the position of each body according to its velocity.
    for (Body &body : bodies) {
      body.x += dt * body.vx;
      body.y += dt * body.vy;
      body.z += dt * body.vz;
    }
  }

  return 0;
}