Barnes-Hut N-body simulator

I have written an n-body simulator, implementing the Barnes-Hut algorithm. Please comment on anything you can see wrong with this.

Wikipedia Barnes-Hut page

This is a screen shot of the simulation 20 hours in. All the particles spawn in a uniform disk, given an initial velocity in order to "orbit" the "Galactic center" (an invisible object at the center of the simulation, and therefore a galaxy, with GalaticCenterMass mass). They then swirl around and around like a galaxy does and because the particles attract each other they form wonderful pictures and time lapse animations.

Node.h:

#pragma once
#include 

struct Body
{
    float posX, posY;           //position x and y
    float velX, velY;           //velocity x and y
    double forceX, forceY;      //force acting on object since last frame
    float mass;                 //mass of object
};

class Node
{
public:
    std::vector Bodies;
    std::vector Child;
    bool HasChildren;

    float posX, posY;
    float width, height;
    float TotalMass;
    float CenterOfMassx;
    float CenterOfMassy;
    unsigned int Depth;
    bool IsUsed;            //For testing, delete this later

    Node();
    Node(unsigned int pdepth);
    ~Node();
    void GenerateChildren();
    void SetParam(std::vector pBodies, float pwidth, float pheight, float px = 0, float py = 0);
    void Reset();
};

Node.cpp:

```
#include "Node.h"
#include

#define _MAX_DEPTH 100

Node::Node(unsigned int pDepth)
{
HasChildren = false;
Depth = pDepth;
}

Node::Node()
{
HasChildren = false;
Depth = 0;
}

void Node::GenerateChildren()
{
std::vector q1Bodies, q2Bodies, q3Bodies, q4Bodies;

for (unsigned int i = 0; i posX posY posY SetParam(q1Bodies, width / 2, height / 2, posX, posY);
q2->SetParam(q2Bodies, width / 2, height / 2, posX + (width / 2), posY);
q3->SetParam(q3Bodies, width / 2, height / 2, posX, posY + (height / 2));
q4->SetParam(q4Bodies, width / 2, height / 2, posX + (width / 2), posY +

The code was big so my review of it will probably not be that complete. I'll try to go in order of top to bottom of the code that was listed.

Force vs Acceleration

In your body structure, you track forceX and forceY. It would be better to track accelerationX and accelerationY, because you will save yourself some computation. When you compute the force acting on a body due to gravity, the formula is this:

double force = (bj->mass * bi->mass * invDistCube * _GRAV_CONST);

And then when you update the position, you do this:

pBodies.at(i)->velX += pBodies.at(i)->forceX / pBodies.at(i)->mass;

If you notice, you multiplied the mass of the body into the force and then divided it out later on. If you computed the acceleration, you would skip those two steps.

double accel = (bj->mass * invDistCube * _GRAV_CONST);
 ...
 pBodies.at(i)->velX += pBodies.at(i)->accelX;

I did notice that you used the force to determine the color. To maintain that behavior, you could either multiply the acceleration by the mass to recompute the force (still saving a division). Or you could switch to coloring by acceleration, which may be interesting to see.

Building the quadtree

Bounding box

When you build the quadtree, one thing you could do is to compute the actual bounding box for each Node instead of using the bounding box that was passed in. Right now for example, the top level nodes of the quadtree are the four quadrants of your screen. But if the bodies within a quadrant only occupy a small part of that quadrant, you are still using the full size of the quadrant.

Computing the bounding box is easy, since you already iterate through all the bodies when you create a node. By using the actual bounding box, you could get a more accurate measurement of the "threshold" that you use to determine if a quadrant is far enough away from a body to be considered a single body.

Of course if you do this, you shouldn't just use width to be the quandrant size as you are doing now. You would need to compute something like the diagonal width of the quadrant to be your quadrant size.

Dividing point of child quadrants

Currently you generate your child quadrants by using the center of the parent quadrant as the dividing point. If you compute the actual bounding box as mentioned above, and then use the center of the bounding box as the dividing point, it will probably divide your bodies a bit better than before.

But another idea would be to use the center of mass of the quadrant as your dividing point. I'm not entirely sure this would work better than the geometric center, but you could try it and find out.

Node::SetParam()

You have this variable:

unsigned int size = pBodies.size();

But you still use pBodies.size() two other places in the function. Also this line of code:

CenterOfMassx = Centerx / size;

looks bad to me because size is often zero. Although you can divide by zero with doubles, I would suggest checking for size == 0 near the top and avoiding most of the code entirely.

Node::Reset()

You can shorten this code:

for (unsigned int i = 0; i Reset();
}

for (unsigned int i = 0; i < Child.size(); i++)
{
    delete Child[i];
}

to this:

for (unsigned int i = 0; i Reset();
    delete Child[i];
}

AttractToCenter(), RepellFromCenter()

First of all, RepellFromCenter() is never used, so it can be removed. As for AttractToCenter(), I would just create an extra body with the galactic mass and let your code handle it like any other body. You just have to make sure that after each iteration, you reposition the galactic mass back to the galactic center if you don't want it to move.

pSoftener

• I don't like the name, because it makes me think that it is a pointer but it actually isn't.

• It can only ever take on the value of the global constant Softener. So you could remove it as a parameter and then just use the global constant instead.

CheckNode()

Here, you check if a body is far enough away from a node to be able to consider the node as a single center of mass. The check looks like this:

float distance = sqrt((diffX) * (diffX) + (diffY) * (diffY));                       

if ((pNode->width / distance) HasChildren == false))

I don't really like the parentheses in the if statement, since all of them are redundant. It could look like this:

if (pNode->width / distance HasChildren)

Beyond that, I think you could save yourself a call to sqrt() by using the distance squared instead of the distance. So in other words, instead of:

if (nodeWidth / distance < _NODE_THRESHOLD)

you could do:

if (nodeWidthSquared / distanceSquared < _NODE_THRESHOLD_SQUARED)

CalculateForceNode() vs. CalculateForce()

These are very nearly the same function. I would suggest that each Node should contain a Body which holds the center of mass information for that node. That way, you could eliminate `CalculateFor

StackExchange Code Review Q#95932, answer score: 7