n-body problem code in c++ example 2 (Barnes-Hut)

} } }

void printQuadtree(std::vector<Body> &bodies, QuadtreeNode *node) { // Iterate over the bodies in the region represented by the node. std::cout << "(" << node->cx << "," << node->cy << "," << node->cz << ") "; for (Body &body : bodies) { if (body.x >= node->cx) std::cout << "x=" << body.x << ", y=" << body.y << ", z=" << body.z << std::endl; if (body.y >= node->cy) std::cout << "y=" << body.y << ", z=" << body.z << std::endl; if (body.z >= node->cz) std::cout << "z=" << body.z << std::endl; } }

The code includes the following declarations:

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

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

const double G = 6.674e-11;  // gravitational constant
const double THETA = 0.5;  // accuracy parameter for Barnes-Hut algorithm

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

struct QuadtreeNode {
  double xmin, xmax, ymin, ymax, zmin, zmax;  // bounding box
  double cx, cy, cz;  // center of mass
  double mass;  // total mass
  QuadtreeNode *child[8];  // child nodes
};

// Build the quadtree by recursively dividing the space into regions and adding bodies to the appropriate nodes.
void buildQuadtree(std::vector<Body> &bodies, QuadtreeNode *node, double xmin, double xmax, double ymin, double ymax, double zmin, double zmax) {
  // Initialize the node with the bounding box and empty child nodes.
  node->xmin = xmin;
  node->xmax = xmax;
  node->ymin = ymin;
  node->ymax = ymax;
  node->zmin = zmin;
  node->zmax = zmax;
  for (int i = 0; i < 8; ++i) {
    node->child[i] = nullptr;
  }

  // Calculate the center of mass and total mass of the bodies within the region represented by the node.
  double cx = 0.0, cy = 0.0, cz = 0.0;
  double mass = 0.0;
  for (Body &body : bodies) {
    if (body.x >= xmin && body.x < xmax && body.y >= ymin && body.y < ymax && body.z >= zmin && body.z < zmax) {
      cx += body.x * body.mass;
      cy += body.y * body.mass;
      cz += body.z * body.mass;
      mass += body.mass;
    }
  }
  node->cx = cx / mass;
  node->cy = cy / mass;
  node->cz = cz / mass;
  node->mass = mass;

  // If the region contains more than one body, divide it into 8 subregions and build child nodes for each subregion.
  if (bodies.size() > 1) {
    std::vector<Body> subbodies[8];
    for (Body &body : bodies) {
      if (body.x >= xmin && body.x < xmax && body.y >= ymin && body.y < ymax && body.z >= zmin && body.z < zmax) {
        int index = 0;
        if (body.x >= node->cx) index += 1;
        if (body.y >= node->cy) index += 2;
        if (body.z >= node->cz) index += 4;
        subbodies[index].push_back(body);
      }
    }
    double subxmin, subxmax, subymin, subymax, subzmin, subzmax;
    for (int i = 0; i < 8; ++i) {
      subxmin = i & 1 ? node->cx : xmin;
      subxmax = i & 1 ? xmax : node->cx;
      subymin = i & 2 ? node->cy : ymin;
      subymax = i & 2 ? ymax : node->cy;
      subzmin = i & 4 ? node->cz : zmin;
      subzmax = i & 4 ? zmax : node->cz;
      if (!subbodies[i].empty()) {
        node->child[i] = new QuadtreeNode;
        buildQuadtree(subbodies[i], node->child[i], subxmin, subxmax, subymin, subymax, subzmin, subzmax);
    }
  }
  }