# 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 <vector>

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<Body*> Bodies;
std::vector<Node*> 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<Body*> pBodies, float pwidth, float pheight, float px = 0, float py = 0);
void Reset();
};


Node.cpp:

#include "Node.h"
#include <memory>

#define _MAX_DEPTH 100

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

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

void Node::GenerateChildren()
{
std::vector<Body*> q1Bodies, q2Bodies, q3Bodies, q4Bodies;

for (unsigned int i = 0; i < Bodies.size(); i++)
{
if (Bodies[i]->posX < (posX + (width / 2)))     //if true, 1st or 3rd
{
if (Bodies[i]->posY < (posY + (height / 2)))    //1
{
q1Bodies.push_back(Bodies[i]);
}
else //3
{
q3Bodies.push_back(Bodies[i]);
}
}
else                                            //2 or 4
{
if (Bodies[i]->posY < (posY + (height / 2)))    //2
{
q2Bodies.push_back(Bodies[i]);
}
else //4
{
q4Bodies.push_back(Bodies[i]);
}
}
}

Node* q1 = new Node(Depth + 1);
Node* q2 = new Node(Depth + 1);
Node* q3 = new Node(Depth + 1);
Node* q4 = new Node(Depth + 1);

q1->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 + (height / 2));

Child.push_back(q1);
Child.push_back(q2);
Child.push_back(q3);
Child.push_back(q4);

HasChildren = true;
}

void Node::SetParam(std::vector<Body*> pBodies, float pwidth, float pheight, float px, float py)
{
Bodies = pBodies;
posX = px;
posY = py;
width = pwidth;
height = pheight;

float mass = 0;
double Centerx = 0;
double Centery = 0;

for (unsigned int i = 0; i < pBodies.size(); i++)
{
mass += pBodies[i]->mass;
Centerx += pBodies[i]->posX;
Centery += pBodies[i]->posY;
}

TotalMass = mass;

unsigned int size = pBodies.size();

CenterOfMassx = Centerx / size;
CenterOfMassy = Centery / size;

if (Bodies.size() > 1 && Depth < _MAX_DEPTH)
{
GenerateChildren();
}
}

void Node::Reset()
{
Bodies.clear();

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

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

Child.clear();

HasChildren = false;
}

Node::~Node()
{

}


And finally the main.cpp file:

#include "SFML\Graphics.hpp"
#include "Node.h"
#include <ctime>
#include <math.h>
#include <vector>

#define _PI 3.14159265      //Pi, used for calculations and rounded to 8 decimal places.
#define _GRAV_CONST 0.1     //the gravitational constant. This is the timestep between each frame. Lower for slower but more accurate simulations

void BodyAttraction(std::vector<Body*> pBodies, float pSoftener);                                                                       //Attracts each body to each other body in the given vector of pointers to body objects
void CalculateForceNode(Body* bi, Node* bj, float pSoftener);                                                                           //Calculate force exerted on body from node
void CalculateForce(Body* bi, Body* bj, float pSoftener);                                                                               //Calculate force exerted on eachother between two bodies
Body* CreateBody(float px, float py, float pmass, float pvx = 0, float pvy = 0);                                                        //return a pointer to new body object defined on the heap with given paramiters
void DeleteBodies(std::vector<Body*> pBodies);                                                                                          //Deletes objects pointed to by given vector
void PollEvent(sf::RenderWindow* pTarget, bool* pIsPaused, sf::View* pSimView);                                                         //Call all polled events for the sf::window
void PopulateBodyVectorDisk(std::vector<Body*> *pBodies, unsigned int pParticlesCount, float pWidth, float pHeight, float pMaxDist, float pMinDist, float pMinMass, float pMaxMass, float pGalaticCenterMass = 0);  //populate given vector with bodies with given paramiters in a disk formation
void PopulateBodyVectorUniform(std::vector<Body*> *pBodies, unsigned int pParticlesCount, float pWidth, float pHeight, float pMinMass, float pMaxMass);
void Render(sf::RenderWindow* pTarget, std::vector<Body*> pBodies, sf::Color pObjColor);                                                //Render given body objects to given screen
void SetView(sf::View* pView, sf::RenderWindow* pTarget, float pViewWidth, float pViewHeight);                                          //set the window to the simulation view
void UpdateBodies(std::vector<Body*> pBodies);                                                                                          //Calculate velocity chance from the bodies force exerted since last update, update position based on velocity, reset force to 0
void DrawNode(Node* pNode, sf::RenderWindow* pTarget);                                                                                  //Draw a node to the screen, and all of its children (recursive)
void CheckNode(Node* pNode, Body* pBody);                                                                                               //Checks if a node is sufficently far away for force calculation, if not recureses on nodes children
void OctreeBodyAttraction();                                                                                                            //Using a calculated oct-tree, calculate the force exerted on each object
void AttractToCenter(std::vector<Body*> pBodies, float width, float height, float centerMass);                                                                                                                  //Attract each particle to the center of the simulation
void ResetForces(std::vector<Body*> pBodies);
void RepellFromCenter(std::vector<Body*> pBodies, float width, float height, float centerMass);

float const DiskRadiusMax = 20000;              //Max and min distances objects will be from the galatic center
float const GalaticCenterMass = 1000000;        //The mass of the very large object simulating a black hole at the center of a galixy;
float const ObjectMassMax = 2;                  //The max and min mass of the objects in the galixy
float const ObjectMassMin = 1;
float const SimWidth = 327680;                  //Width and height of simulation, needs to be large, particles outside of this range will not be included in the octree
float const SimHeight = 327680;
float const ViewWidth = 1920;                   //Width and height of view of the simulation for the screen.
float const ViewHeight = 1080;
float const Softener = 10;                      //A softener used for the force calculations, 10 is a good amount
unsigned int const NumParticles = 10000;        //Number of particles in simtulation, currently 2^15
double const _NODE_THRESHOLD = 0.5;             //Threshold for node calculations

float zoom = 1;                                 //The current amount of zoom in or out the user has inputed in total
float MouseX = 0;
float MouseY = 0;

std::vector<Body*> Bodies;                      //Container of all Bodies in simulation
Node GlobalNode;
bool IsPaused = false;                          //Contains the state of weather the simulation is paused or not
sf::Color ObjColor(255, 255, 255, 128);         //the defult colour of the objects
sf::View SimulationView;
sf::RenderWindow window(sf::VideoMode(1920, 1080), "N-Body simulation");

int main()
{
SetView(&SimulationView, &window, ViewWidth, ViewHeight);

while (window.isOpen())
{
PollEvent(&window, &IsPaused, &SimulationView); //These will always be done

if (!IsPaused)  //These will not if the simulation is paused
{
AttractToCenter(Bodies, SimWidth, SimHeight, GalaticCenterMass);
UpdateBodies(Bodies);
ResetForces(Bodies);
GlobalNode.Reset();
GlobalNode.SetParam(Bodies, SimWidth, SimHeight);
OctreeBodyAttraction();
}

Render(&window, Bodies, ObjColor);
}

DeleteBodies(Bodies);
}

void AttractToCenter(std::vector<Body*> pBodies, float width, float height, float centerMass)
{
Body* Temp = CreateBody(width / 2, height / 2, centerMass); //Create a body at the center of the simulation

for (unsigned int i = 0; i < pBodies.size(); i++)
{
CalculateForce(pBodies[i], Temp, Softener);
}

delete Temp;
}

void RepellFromCenter(std::vector<Body*> pBodies, float width, float height, float centerMass)
{
Body* Temp = CreateBody(width / 2, height / 2, centerMass); //Create a body at the center of the simulation

for (unsigned int i = 0; i < pBodies.size(); i++)
{
float vectorx = Temp->posX - pBodies[i]->posX;
float vectory = Temp->posY - pBodies[i]->posY;

float distSqr = vectorx * vectorx + vectory * vectory;

double Dist = (sqrt(distSqr));

double force = (pBodies[i]->mass * Dist * _GRAV_CONST * 0.0001);

pBodies[i]->forceX -= vectorx * force;
pBodies[i]->forceY -= vectory * force;
}

delete Temp;
}

void ResetForces(std::vector<Body*> pBodies)
{
for (unsigned int i = 0; i < pBodies.size(); i++)
{
pBodies[i]->forceX = 0;
pBodies[i]->forceY = 0;
}
}

void BodyAttraction(std::vector<Body*> pBodies, float pSoftener)
{
for (unsigned int i = 0; i < pBodies.size(); i++)
{
for (unsigned int j = 0; j < pBodies.size(); j++)
{
CalculateForce(pBodies.at(i), pBodies.at(j), pSoftener); //for each body in pBodies: each other body in pBodies: Calculate attractive force exerted on the first body from the second one
}
}
}

void OctreeBodyAttraction()
{
for (unsigned int i = 0; i < Bodies.size(); i++)
{
CheckNode(&GlobalNode, Bodies[i]);
}
}

inline void CheckNode(Node* pNode, Body* pBody)
{
if (pNode->Bodies.size() != 0)
{
float diffX = (pNode->CenterOfMassx - pBody->posX);
float diffY = (pNode->CenterOfMassy - pBody->posY);

float distance = sqrt((diffX) * (diffX) + (diffY) * (diffY));   //Distance from the node to the object

if ((pNode->width / distance) < (_NODE_THRESHOLD) || (pNode->HasChildren == false))     //if sufficently far away or has no children (external node) group node and calculate force
{
CalculateForceNode(pBody, pNode, Softener);
pNode->IsUsed = true;
}
else                                                                                    //if not, repeat function with children
{
CheckNode(pNode->Child[0], pBody);
CheckNode(pNode->Child[1], pBody);
CheckNode(pNode->Child[2], pBody);
CheckNode(pNode->Child[3], pBody);
}
}
}

inline void CalculateForceNode(Body* bi, Node* bj, float pSoftener)  //bi is being attracted to bj. 15 flops of calculation
{
//vector from the body to the center of mass
float vectorx = bj->CenterOfMassx - bi->posX;
float vectory = bj->CenterOfMassy - bi->posY;

//c^2 = a^2 + b^2 + softener^2
float distSqr = vectorx * vectorx + vectory * vectory + pSoftener * pSoftener;

// ivnDistCube = 1/distSqr^(3/2)
float distSixth = distSqr * distSqr * distSqr;
double invDistCube = 1.0f / (sqrt(distSixth));

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

bi->forceX += vectorx * force;
bi->forceY += vectory * force;
}

inline void CalculateForce(Body* bi, Body* bj, float pSoftener)  //bi is being attracted to bj. 15 flops of calculation
{
//std::vector from i to j
float vectorx = bj->posX - bi->posX;
float vectory = bj->posY - bi->posY;

//c^2 = a^2 + b^2 + softener^2
float distSqr = vectorx * vectorx + vectory * vectory + pSoftener * pSoftener;

// ivnDistCube = 1/distSqr^(3/2)
float distSixth = distSqr * distSqr * distSqr;
double invDistCube = 1.0f / (sqrt(distSixth));

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

bi->forceX += vectorx * force;
bi->forceY += vectory * force;
}

Body* CreateBody(float px, float py, float pmass, float pvx, float pvy)
{
Body* Temp = new Body;

Temp->posX = px;
Temp->posY = py;
Temp->mass = pmass;
Temp->velX = pvx;
Temp->velY = pvy;
Temp->forceX = Temp->forceY = 0;

return Temp;
}

void DeleteBodies(std::vector<Body*> pBodies)
{
for (unsigned int i = 0; i < pBodies.size(); i++)
{
delete pBodies[i];
}

pBodies.clear();
}

void PollEvent(sf::RenderWindow* pTarget, bool* pIsPaused, sf::View* pSimView)
{
sf::Event event;

while (pTarget->pollEvent(event))
{
if (event.type == sf::Event::Closed)
pTarget->close();
if (event.type == sf::Event::KeyPressed)
{
if (event.key.code == sf::Keyboard::Space)
*pIsPaused = !*pIsPaused;                   //toggle what is pointed to by IsPaused
}
if (event.type == sf::Event::MouseWheelScrolled)
{
zoom *= 1 + (static_cast <float> (-event.mouseWheelScroll.delta) / 10); //for each notch down, -10%, for each notch up, +10%
pSimView->zoom(1 + (static_cast <float> (-event.mouseWheelScroll.delta) / 10));
}
}

if (sf::Mouse::getPosition().x > (1920 - 20))
SimulationView.move(2 * zoom, 0);
if (sf::Mouse::getPosition().x < (0 + 20))
SimulationView.move(-2 * zoom, 0);
if (sf::Mouse::getPosition().y > (1080 - 20))
SimulationView.move(0, 2 * zoom);
if (sf::Mouse::getPosition().y < (0 + 20))
SimulationView.move(0, -2 * zoom);

pTarget->setView(*pSimView);
}

void PopulateBodyVectorDisk(std::vector<Body*> *pBodies, unsigned int pParticlesCount, float pWidth, float pHeight, float pMaxDist, float pMinDist, float pMinMass, float pMaxMass, float pGalaticCenterMass)
{
srand(static_cast<unsigned int> (time(0)));

for (unsigned int i = 0; i < pParticlesCount; i++)
{
float angle = static_cast <float> (rand()) / (static_cast <float> (RAND_MAX / (2 * static_cast <float> (_PI))));    //sets angle to random float range (0, 2 pi)

float distanceCoefficent = static_cast <float> (rand()) / static_cast <float> (RAND_MAX);
float distance = pMinDist + ((pMaxDist - pMinDist) * (distanceCoefficent * distanceCoefficent));                    //Distance point will be from the galatic center, between MinDiskRadius and MaxDiskRadius

float positionx = cos(angle) * distance + (pWidth / 2);                                                             //set positionx and positiony to be the point you get when you go in the direction of 'angle' till you have traveled 'distance'
float positiony = sin(angle) * distance + (pHeight / 2);

float orbitalVelocity = sqrt((pGalaticCenterMass * static_cast <float> (_GRAV_CONST)) / distance);                  //Calculate the orbital velocity required to orbit the galatic centre

float velocityx = (sin(angle) * orbitalVelocity);
float velocityy = (-cos(angle) * orbitalVelocity);

float mass = pMinMass + static_cast <float> (rand() % static_cast <int> (pMaxMass - pMinMass));                     //random mass (int) in range (MinObjectMass, MaxObjectMass)

pBodies->push_back(CreateBody(positionx, positiony, mass, velocityx, velocityy));
}
}

void PopulateBodyVectorUniform(std::vector<Body*> *pBodies, unsigned int pParticlesCount, float pWidth, float pHeight, float pMinMass, float pMaxMass)
{
srand(static_cast<unsigned int> (time(0)));

for (unsigned int i = 0; i < pParticlesCount; i++)
{
float positionx = static_cast <float> (rand()) / static_cast <float> (RAND_MAX) * pWidth + (SimWidth / 2 - pWidth / 2);
float positiony = static_cast <float> (rand()) / static_cast <float> (RAND_MAX) * pHeight + (SimHeight / 2 - pHeight / 2);

float mass = pMinMass + static_cast <float> (rand() % static_cast <int> (pMaxMass - pMinMass));                     //random mass (int) in range (MinObjectMass, MaxObjectMass)

pBodies->push_back(CreateBody(positionx, positiony, mass));
}
}

void Render(sf::RenderWindow* pTarget, std::vector<Body*> pBodies, sf::Color pObjColor)
{
pTarget->clear();

sf::RectangleShape Temp;
//Temp.setFillColor(pObjColor);

for (unsigned int i = 0; i < pBodies.size(); i++)
{
if (zoom > 1)
Temp.setSize(sf::Vector2f(pBodies.at(i)->mass * zoom, pBodies.at(i)->mass * zoom));
else
Temp.setSize(sf::Vector2f(pBodies.at(i)->mass, pBodies.at(i)->mass));

float ForceCoiffecent = sqrt(pBodies.at(i)->forceX * pBodies.at(i)->forceX + pBodies.at(i)->forceY * pBodies.at(i)->forceY) * (40000 * _GRAV_CONST);

if (ForceCoiffecent > 1)
ForceCoiffecent = 1;

float Red, Green, Blue;

Blue = 1 - (ForceCoiffecent);

if (ForceCoiffecent < 0.2)
Red = (ForceCoiffecent) * 5;
else
Red = 1;

if (ForceCoiffecent < 0.5)
Green = (ForceCoiffecent) * 2;
else
Green = 1;

Temp.setFillColor(sf::Color(Red * 255, Green * 255, Blue * 255, 128));
Temp.setPosition(pBodies.at(i)->posX, pBodies.at(i)->posY);
pTarget->draw(Temp);
}

//DrawNode(&GlobalNode, pTarget);

pTarget->display();
}

void DrawNode(Node* pNode, sf::RenderWindow* pTarget)
{
sf::RectangleShape Temp;
Temp.setFillColor(sf::Color(0, 0, 0, 0));
Temp.setOutlineThickness(zoom);
Temp.setOutlineColor(sf::Color(0, 255, 0, 16));
Temp.setPosition(pNode->posX, pNode->posY);
Temp.setSize(sf::Vector2f(pNode->width, pNode->height));

pTarget->draw(Temp);
if (pNode->HasChildren) //recercivly draw all children
{
DrawNode(pNode->Child[0], pTarget);
DrawNode(pNode->Child[1], pTarget);
DrawNode(pNode->Child[2], pTarget);
DrawNode(pNode->Child[3], pTarget);
}
}

void SetView(sf::View* pView, sf::RenderWindow* pTarget, float pViewWidth, float pViewHeight)
{
pView->reset(sf::FloatRect(SimWidth / 2 - pViewWidth / 2, SimHeight / 2 - pViewHeight / 2, pViewWidth, pViewHeight));
pView->setViewport(sf::FloatRect(0.f, 0.f, 1.f, 1.f));
pTarget->setView(*pView);
}

void UpdateBodies(std::vector<Body*> pBodies)
{
for (unsigned int i = 0; i < pBodies.size(); i++)
{
if ((pBodies[i]->posX > SimWidth && pBodies[i]->velX > 0) || (pBodies[i]->posX < 0 && pBodies[i]->velX < 0))
{
pBodies[i]->velX = -pBodies[i]->velX;
}

if ((pBodies[i]->posY > SimHeight && pBodies[i]->velY > 0) || (pBodies[i]->posY < 0 && pBodies[i]->velY < 0))
{
pBodies[i]->velY = -pBodies[i]->velY;
}

pBodies.at(i)->velX += pBodies.at(i)->forceX / pBodies.at(i)->mass;     //f = ma => force = mass * accelleration. Therefor
pBodies.at(i)->velY += pBodies.at(i)->forceY / pBodies.at(i)->mass;     //a = f/m => accelleration = force / mass

pBodies.at(i)->posX += pBodies.at(i)->velX;
pBodies.at(i)->posY += pBodies.at(i)->velY;
}
}

• Welcome to Code Review! I'm sure there are things that can be improved, I hope you get some good answers! – Phrancis Jul 6 '15 at 3:24
• Thanks, I'm a little scared of a down vote backlash, Ive seen questions over on stack overflow that just dump all there code and ask for help get down voted to oblivion, But i'm guessing a question of this fashion is appropriate here? – Kieren Pearson Jul 6 '15 at 3:26
• Please comment on anything you can see wrong with this. - Perfect CR question, as long as everything works. Also, you explained what it does, so you're good. – user34073 Jul 6 '15 at 3:28
• @KierenPearson Don't be scared, we don't bite! You might improve your question a bit by telling us a bit more about how it works, but either way, reviewing large amounts of code is what we do here :) – Phrancis Jul 6 '15 at 3:29
• The amount of code is not excessive. Give it a day or so, then stop by chat if you don't get an answer. – 200_success Jul 6 '15 at 6:14

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.

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.

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 < Child.size(); i++)
{
Child[i]->Reset();
}

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


to this:

for (unsigned int i = 0; i < Child.size(); i++)
{
Child[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

1. I don't like the name, because it makes me think that it is a pointer but it actually isn't.
2. 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) < (_NODE_THRESHOLD) || (pNode->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 < _NODE_THRESHOLD || !pNode->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 CalculateForceNode() and just pass in node->centerOfMassBody to CalculateForce() instead. If you did that, you could also eliminate node->TotalMass, node->CenterOfMassx, and node->CenterOfMassy since all that information would be in node->centerOfMassBody already.

## Miscellaneous

I only skimmed the rest of the code. So here are some other small things:

1. Misspellings such as Coiffecent vs Coefficient, Repell vs Repel.
2. Extra parentheses in expressions (which possibly aren't that bad, it's a style thing).
3. You defined your own pi constant, but you could just use M_PI.

1. As far as non-square nodes are concerned, I don't see any reason why rectangular nodes can't be used as long as you use the diagonal as the "width" for determining thresholds. But I'm not familiar with the Barnes Hut algorithm so there may be some theoretical reason why square nodes are better.

2. For particles ejected from the Galaxy, I suppose the bounding boxes would increase in size, but I doubt that would matter much. You could always exclude outliers, as you said.

3. I believe that SetParams() can be called with size zero vectors. Imagine if GenerateChildren() were called with 2 bodies. It creates 4 new nodes, two of which are empty and two of which have 1 body. All four of these nodes have SetParams() called on them, so two of these will have size zero vectors.

4. I know that a/b is not the same as a^2/b^2, but I was comparing to the threshold squared also. For positive a and b, this relation holds: if a/b < 0.5, then a^2/b^2 < 0.25.

5. I thought that M_PI was defined in math.h, but apparently that is not universally true. I'm used to Gnu, where M_PI is defined.

• Thank you so much for the time you put into this. Your idea to use acceleration not force is very good, I don't know how you saw that. For the bounding box section, I cant split it by center of mass, because each node must be square, and such it must be split down the center. For your idea about the bounding box, are you saying the top level node (The one that is the entire simulation) should insted be the smallest it can be to include all the particles? One issue with that is particles can be ejected from the galixy, and then the bounding box could be larger than it started! Cont. – Kieren Pearson Jul 7 '15 at 1:36
• I could impliment a maximum size of the global box, but that size would be reached quickly, I could exclude the top 5% and the botton 5% of particle position values (Remove the outliers) but that causes when the simulation starts the particles on the outscirts of the galixy to not be attracted at all. – Kieren Pearson Jul 7 '15 at 1:39
• "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." Yes this is somthing that delicitly does not occur but I will include more saftey as of now. SetParam does not get called anywhere other than GenerateChildren(), which only gets called from SetParam() if the size of its body vector is greater than one. A delicate situation indeed ;P – Kieren Pearson Jul 7 '15 at 1:45
• if (nodeWidthSquared / distanceSquared < _NODE_THRESHOLD_SQUARED) does not work, because the ratio a/b is not equal to the ratio a^2/b^2 and the only way to use a^2/b^2 is to afterwards sqrt(a^2/b^2) which just returns a/b and your back to where you started, but you needed a extra multiplication to square the pNode->width – Kieren Pearson Jul 7 '15 at 2:07
• @KierenPearson No, don't post updated code as an edit. What you can do is make a new question after you have made all your fixups, and ask for a rereview (if you want one). If you don't want a rereview, you can maybe add a short addendum to your question describing what you fixed up, or do it in comments. – JS1 Jul 7 '15 at 3:00

This is an interesting piece of code, and a well written question, so I'm sure you'll get a number of good reviews. Here are some observations that may help you improve your code.

The Node and Body objects both have exclusively public members, which is not very good design. Every piece of code that interacts with them then reaches inside and manipulates the object members directly which is prone to error and makes maintenance difficult. So as a simple example, let's consider Body. First, you could make it into a class rather than a struct. Next, some of the functions within main really should instead be member functions of Body. To take a few obvious examples, the CalculateForce and CalculateForceNode should instead be member functions:

// updates force based on other body
void CalculateForce(const Body* bj, float pSoftener);
// updates force based on other node
void CalculateForce(const Node* bj, float pSoftener);


Note also that the other Body or Node is not adjusted, so it is declared as const.

It would also be a good idea to encapsulate Bodies in main.cpp as an object with member functions rather than passing it around to stand-alone functions.

## Use constructors

The CreateBody function in main should instead be a constructor for Body:

Body(float px, float py, float pmass, float pvx=0, float pvy=0) :
posX(px),
posY(py),
velX(pvx),
velY(pvy),
forceX(0),
forceY(0),
mass(pmass)
{}


This also allows the temporary Body created in AttractToCenter (and the unused RepellFromCenter) to be created on the stack instead of the heap, avoiding both new and delete in that function.

## Use references to pass complex objects

Most of the functions in main which take a std::vector<Body*> as the first argument should instead be declared to take a reference instead. So for example,

void AttractToCenter(std::vector<Body*> pBodies, float width, float height, float centerMass)


becomes

void AttractToCenter(std::vector<Body*> &pBodies, float width, float height, float centerMass)


The difference is just the single & character in front of pBodies but the performance impact is large because it prevents the compiler from having to actually make a copy of the vector.

With all of the above suggestion and this one implemented, that function becomes quite small and simple:

void AttractToCenter(std::vector<Body*> &pBodies, float width, float height, float centerMass)
{
Body Temp(width / 2, height / 2, centerMass); //Create a body at the center of the simulation
for (unsigned int i = 0; i < pBodies.size(); i++)
{
pBodies[i]->CalculateForce(&Temp, Softener);
}
}


It could be even simpler, if your compiler isn't too old, as shown in the next suggestion.

## Use "range-for" to simplify your code

If your compiler fully supports the C++11 standard, you can use a "range-for" to simplify your code. The loop above, for example, can be written like this:

for (auto &pB : pBodies)
{
pB->CalculateForce(&Temp, Softener);
}


## Avoid the use of global variables

Most of the global variables that are not constants should probably be restricted to a smaller scope. For instance, IsPaused can simply be moved into main -- it does not have to be a global variable. Doing so will reduce the hidden linkages between functions and help you strengthen your interface and design and reduce the difficulty of maintenance.

## Use forward slashes in include paths

Using this path for your include file only works under Windows:

#include "SFML\Graphics.hpp"


If you use a forward slash instead, it works on all platforms.

#include "SFML/Graphics.hpp"


## Prefer const values to #defines

Both _GRAV_CONST and _PI should be declared as const float values rather than as #defines. Doing so allows for better type safety and in this particular case, allows for the elimination of a messy static_cast.

## Consider using a better random number generator

If you are using a compiler that supports at least C++11, consider using a better random number generator. In particular, instead of rand, you might want to look at std::uniform_real_distribution and friends in the <random> header.

## Eliminate memory leaks

There are a number of places where new is called without any corresponding delete calls. This is the recipe for a memory leak. Either fix the leak by adding specific calls to delete, or even better, use smart pointers to eliminate new and delete entirely. One basic thing that would help a great deal is to decide which objects actually own the Body and Node items. That will help you to decide if and when they should be deleted.