Link to previous question: Barnes-Hut N-body simulator
It still has issues I'm sure.
Node.h
#pragma once
#include <vector>
struct Body
{
double posX, posY; //position x and y
double velX, velY; //velocity x and y
double AccelX, AccelY; //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;
Node();
Node(std::vector<Body*> pBodies, float pwidth, float pheight, float px = 0, float py = 0, unsigned int pDepth = 0);
~Node();
void GenerateChildren();
void SetParam(std::vector<Body*> pBodies, float pwidth, float pheight, float px = 0, float py = 0, unsigned int pDepth = 0);
void Reset();
};
Node.cpp
#include "Node.h"
#include <memory>
#define _MAX_DEPTH 100
Node::Node(std::vector<Body*> pBodies, float pwidth, float pheight, float px, float py, unsigned int pDepth)
{
SetParam(pBodies, pwidth, pheight, px, py, pDepth);
}
void Node::SetParam(std::vector<Body*> pBodies, float pwidth, float pheight, float px, float py, unsigned int pDepth)
{
HasChildren = false;
Depth = pDepth;
unsigned int size = pBodies.size();
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 < size; i++)
{
mass += pBodies[i]->mass;
Centerx += pBodies[i]->posX;
Centery += pBodies[i]->posY;
}
TotalMass = mass;
if (size > 0)
{
CenterOfMassx = static_cast<float>(Centerx / size);
CenterOfMassy = static_cast<float>(Centery / size);
}
else
{
CenterOfMassx = 0;
CenterOfMassy = 0;
}
if (size > 1 && Depth < _MAX_DEPTH)
{
GenerateChildren();
}
}
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(q1Bodies, width / 2, height / 2, posX, posY, Depth + 1);
Node* q2 = new Node(q2Bodies, width / 2, height / 2, posX + (width / 2), posY, Depth + 1);
Node* q3 = new Node(q3Bodies, width / 2, height / 2, posX, posY + (height / 2), Depth + 1);
Node* q4 = new Node(q4Bodies, width / 2, height / 2, posX + (width / 2), posY + (height / 2), Depth + 1);
Child.push_back(q1);
Child.push_back(q2);
Child.push_back(q3);
Child.push_back(q4);
HasChildren = true;
}
void Node::Reset()
{
Bodies.clear();
for (unsigned int i = 0; i < Child.size(); i++)
{
Child[i]->Reset();
delete Child[i];
}
Child.clear();
HasChildren = false;
}
Node::~Node()
{
}
Main.cpp
#include "SFML/Graphics.hpp"
#include "Node.h"
#include <ctime>
#include <vector>
void BodyAttraction(std::vector<Body*> &pBodies); //Attracts each body to each other body in the given vector of pointers to body objects
void CalculateForceNode(Body* bi, Node* bj); //Calculate force exerted on body from node
void CalculateForce(Body* bi, Body* bj); //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, sf::View* pSimView, float &pZoom, bool* pIsPaused); //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 PopulateBodyVectorExplosion(std::vector<Body*> *pBodies, unsigned int pParticlesCount, float pWidth, float pHeight, float pMinMass, float pMaxMass, float ExplosionForce);
void Render(sf::RenderWindow* pTarget, std::vector<Body*> &pBodies, float pZoom); //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, float pZoom); //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(std::vector<Body*> &pBodies, Node *pGlobalNode); //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 ResetAccel(std::vector<Body*> &pBodies);
void RepelFromCenter(std::vector<Body*> &pBodies, float width, float height, float centerMass);
double const _PI = 3.14159265; //Pi, used for calculations and rounded to 8 decimal places.
double const NodeThresholdSqr = 0.5 * 0.5; //Threshold for node calculations
double const _GRAV_CONST = 0.1; //the gravitational constant. This is the timestep between each frame. Lower for slower but more accurate simulations
float const Softener = 10; //A softener used for the force calculations, 10 is a good amount
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;
int main()
{
unsigned int const NumParticles = 10000; //Number of particles in simtulation, currently 2^15
float const DiskRadiusMax = 20000; //Max and min distances objects will be from the galatic center
float const DiskRadiusMin = 50;
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 ViewWidth = 1920; //Width and height of view of the simulation for the screen.
float const ViewHeight = 1080;
Node GlobalNode;
std::vector<Body*> Bodies; //Container of all Bodies in simulation
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");
bool IsPaused = false; //Contains the state of weather the simulation is paused or not
float zoom = 1; //The current amount of zoom in or out the user has inputed in total
srand(static_cast<unsigned int> (time(NULL)));
PopulateBodyVectorDisk(&Bodies, NumParticles, SimWidth, SimHeight, DiskRadiusMax, DiskRadiusMin, ObjectMassMin, ObjectMassMax, GalaticCenterMass);
SetView(&SimulationView, &window, ViewWidth, ViewHeight);
while (window.isOpen())
{
PollEvent(&window, &SimulationView, zoom, &IsPaused); //These will always be done
if (!IsPaused) //These will not if the simulation is paused
{
AttractToCenter(Bodies, SimWidth, SimHeight, GalaticCenterMass);
UpdateBodies(Bodies);
ResetAccel(Bodies);
GlobalNode.Reset();
GlobalNode.SetParam(Bodies, SimWidth, SimHeight);
OctreeBodyAttraction(Bodies, &GlobalNode);
}
Render(&window, Bodies, zoom);
}
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);
}
delete Temp;
}
void RepelFromCenter(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++)
{
double vectorx = Temp->posX - pBodies[i]->posX;
double vectory = Temp->posY - pBodies[i]->posY;
double distSqr = vectorx * vectorx + vectory * vectory;
double Dist = (sqrt(distSqr));
double accel = (pBodies[i]->mass * Dist * _GRAV_CONST * 0.0001);
pBodies[i]->AccelX -= vectorx * accel;
pBodies[i]->AccelY -= vectory * accel;
}
delete Temp;
}
void ResetAccel(std::vector<Body*> &pBodies)
{
for (unsigned int i = 0; i < pBodies.size(); i++)
{
pBodies[i]->AccelX = 0;
pBodies[i]->AccelY = 0;
}
}
void BodyAttraction(std::vector<Body*> &pBodies)
{
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)); //for each body in pBodies: each other body in pBodies: Calculate attractive force exerted on the first body from the second one
}
}
}
void OctreeBodyAttraction(std::vector<Body*> &pBodies, Node *pGlobalNode)
{
for (unsigned int i = 0; i < pBodies.size(); i++)
{
CheckNode(pGlobalNode, pBodies[i]);
}
}
void CheckNode(Node* pNode, Body* pBody)
{
double diffX = (pNode->CenterOfMassx - pBody->posX);
double diffY = (pNode->CenterOfMassy - pBody->posY);
double distanceSqr = ((diffX) * (diffX) + (diffY) * (diffY)); //Distance from the node to the object
double widthSqr = pNode->width * pNode->width;
if (widthSqr / distanceSqr < NodeThresholdSqr || pNode->HasChildren == false) //if sufficently far away or has no children (external node) group node and calculate force
{
CalculateForceNode(pBody, pNode);
}
else //if not, repeat function with children
{
if (pNode->Child[0]->Bodies.size() > 0)
CheckNode(pNode->Child[0], pBody);
if (pNode->Child[1]->Bodies.size() > 0)
CheckNode(pNode->Child[1], pBody);
if (pNode->Child[2]->Bodies.size() > 0)
CheckNode(pNode->Child[2], pBody);
if (pNode->Child[3]->Bodies.size() > 0)
CheckNode(pNode->Child[3], pBody);
}
}
void CalculateForceNode(Body* bi, Node* bj) //bi is being attracted to bj. 15 flops of calculation
{
//vector from the body to the center of mass
double vectorx = bj->CenterOfMassx - bi->posX;
double vectory = bj->CenterOfMassy - bi->posY;
//c^2 = a^2 + b^2 + softener^2
double distSqr = vectorx * vectorx + vectory * vectory + Softener * Softener;
// ivnDistCube = 1/distSqr^(3/2)
double distSixth = distSqr * distSqr * distSqr;
double invDistCube = 1.0f / (sqrt(distSixth));
double Accel = (bj->TotalMass * invDistCube * _GRAV_CONST);
bi->AccelX += vectorx * Accel;
bi->AccelY += vectory * Accel;
}
void CalculateForce(Body* bi, Body* bj) //bi is being attracted to bj. 15 flops of calculation
{
//std::vector from i to j
double vectorx = bj->posX - bi->posX;
double vectory = bj->posY - bi->posY;
//c^2 = a^2 + b^2 + softener^2
double distSqr = vectorx * vectorx + vectory * vectory + Softener * Softener;
// ivnDistCube = 1/distSqr^(3/2)
double distSixth = distSqr * distSqr * distSqr;
double invDistCube = 1.0f / (sqrt(distSixth));
double Accel = (bj->mass * invDistCube * _GRAV_CONST);
bi->AccelX += vectorx * Accel;
bi->AccelY += vectory * Accel;
}
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->AccelX = Temp->AccelY = 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, sf::View* pSimView, float &pZoom, bool* pIsPaused)
{
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)
{
pZoom *= 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))
pSimView->move(2 * pZoom, 0);
if (sf::Mouse::getPosition().x < (0 + 20))
pSimView->move(-2 * pZoom, 0);
if (sf::Mouse::getPosition().y > (1080 - 20))
pSimView->move(0, 2 * pZoom);
if (sf::Mouse::getPosition().y < (0 + 20))
pSimView->move(0, -2 * pZoom);
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)
{
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 PopulateBodyVectorExplosion(std::vector<Body*> *pBodies, unsigned int pParticlesCount, float pWidth, float pHeight, float pMinMass, float pMaxMass, float ExplosionForce)
{
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 velocity = static_cast <float> (rand()) / static_cast <float> (RAND_MAX);
float positionx = (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 = (pHeight / 2);
float velocityx = (cos(angle) * ExplosionForce * (1 - (velocity * velocity)));
float velocityy = (sin(angle) * ExplosionForce * (1 - (velocity * velocity)));
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)
{
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, float pZoom)
{
pTarget->clear();
sf::RectangleShape Temp;
for (unsigned int i = 0; i < pBodies.size(); i++)
{
if (pZoom > 1)
Temp.setSize(sf::Vector2f(pBodies.at(i)->mass * pZoom, pBodies.at(i)->mass * pZoom));
else
Temp.setSize(sf::Vector2f(pBodies.at(i)->mass, pBodies.at(i)->mass));
float AccelCoefficient = static_cast<float> (sqrt(pBodies.at(i)->AccelX * pBodies.at(i)->AccelX + pBodies.at(i)->AccelY * pBodies.at(i)->AccelY) * (20000 * _GRAV_CONST));
if (AccelCoefficient > 1)
AccelCoefficient = 1;
float Red, Green, Blue;
Blue = 1 - (AccelCoefficient);
if (AccelCoefficient < 0.2)
Red = (AccelCoefficient)* 5;
else
Red = 1;
if (AccelCoefficient < 0.5)
Green = (AccelCoefficient)* 2;
else
Green = 1;
Temp.setFillColor(sf::Color(static_cast<sf::Uint8> (Red * 255), static_cast<sf::Uint8> (Green * 255), static_cast<sf::Uint8> (Blue * 255), 128));
Temp.setPosition(static_cast<float>(pBodies.at(i)->posX), (static_cast<float>(pBodies.at(i)->posY)));
pTarget->draw(Temp);
}
//DrawNode(&GlobalNode, pTarget);
pTarget->display();
}
void DrawNode(Node* pNode, sf::RenderWindow* pTarget, float pZoom)
{
sf::RectangleShape Temp;
Temp.setFillColor(sf::Color(0, 0, 0, 0));
Temp.setOutlineThickness(pZoom);
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, pZoom);
DrawNode(pNode->Child[1], pTarget, pZoom);
DrawNode(pNode->Child[2], pTarget, pZoom);
DrawNode(pNode->Child[3], pTarget, pZoom);
}
}
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)->AccelX;
pBodies.at(i)->velY += pBodies.at(i)->AccelY;
pBodies.at(i)->posX += pBodies.at(i)->velX;
pBodies.at(i)->posY += pBodies.at(i)->velY;
}
}