我编写了一个n体模拟器，实现了Barnes-Hut算法。请对此发表任何评论，对这里有毛病请发表评论。在给定初始速度的情况下，所有粒子都在统一的磁盘中生成，以“绕动”“加拉太中心”（模拟中心处的不可见物体，因此是质量为GalaticCenterMass的星系）。然后它们像银河系一样绕着漩涡旋转，并且由于粒子相互吸引，它们形成了奇妙的图画和时间流逝。

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()
{

}


最后是main.cpp文件：

#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 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 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()
{
    PopulateBodyVectorDisk(&Bodies, NumParticles, SimWidth, SimHeight, DiskRadiusMax, DiskRadiusMin, ObjectMassMin, ObjectMassMax, GalaticCenterMass);
    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;
    }
}


评论

欢迎来到代码审查！我相信有些地方可以改善，希望您能得到一些好的答案！
谢谢，我有点害怕投下赞成票，我已经看到堆栈溢出的问题，这些问题只是将所有代码转储并寻求帮助而被投弃席，但是我猜测这种方式的问题在这里很合适？
请评论任何您可能发现与此有关的错误。 -完善的CR问题，只要一切正常。另外，您还解释了它的作用，所以您很好。
@KierenPearson别害怕，我们不咬人！您可以通过告诉我们更多有关它的工作原理的方法来改善您的问题，但是无论哪种方式，我们在这里要做的都是查看大量代码：)
代码量不要过多。给它大约一天的时间，如果没有答案，请先聊天。

#1 楼

代码很大，所以我对其的评论可能不会那么完整。我将尝试按照列出的代码从上到下的顺序进行操作。

Force vs. Acceleration

在您的body结构中，您跟踪forceX和forceY。最好跟踪accelerationX和accelerationY，因为您可以节省一些计算。计算由于重力作用在物体上的力时，公式为：


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



然后更新位置，您可以执行以下操作：


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



如果注意到，您将身体的质量乘以力，然后将其除上。如果计算加速度，则可以跳过这两个步骤。

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


我确实注意到您用力来确定颜色。为了保持这种行为，您可以将加速度乘以质量以重新计算力（仍保存除法）。或者您可以通过加速切换到着色，这可能会很有趣。

构建四叉树

边界框

构建四叉树，您可以做的一件事是为每个Node计算实际的边界框，而不是使用传入的边界框。例如，现在，四叉树的顶级节点是屏幕的四个象限。但是，如果一个象限内的实体仅占该象限的一小部分，则您仍在使用该象限的全部大小。

计算边界框很容易，因为您已经遍历了所有主体创建节点时。通过使用实际的边界框，您可以更准确地测量“阈值”，该阈值用于确定一个象限是否与物体相距足够远而不能被视为单个物体。

当然，如果执行此操作，则不应该像现在一样使用width作为量化大小。您需要计算象限的对角线宽度，以作为您的象限大小。

子象限的分割点

当前，您使用中心生成子象限以父象限为分割点。如果您按照上述方法计算出实际的边界框，然后将边界框的中心用作分割点，它可能会比以前分割您的身体更好。使用象限的质心作为分割点。我并不完全确定这会比几何中心更好，但是您可以尝试一下并找出答案。

Node :: SetParam（）

：


unsigned int size = pBodies.size();



，但是您仍然在函数中的其他两个地方使用pBodies.size()。另外这行代码：


CenterOfMassx = Centerx / size;
对我来说看起来很糟糕，因为size通常为零。尽管您可以用双精度数除以零，但我建议您检查顶部附近的size == 0，并完全避免使用大部分代码。

Node :: Reset（）

您可以缩短此代码：


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

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



对此：

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


AttractToCenter （），RepellFromCenter（）

首先，从不使用RepellFromCenter()，因此可以将其删除。至于AttractToCenter()，我只是用银河质量创建一个额外的物体，然后让您的代码像处理其他物体一样处理它。您只需要确保在每次迭代后，如果不想让银河质量移动到银河中心，就可以将其重新定位。 />我不喜欢这个名字，因为它使我认为它是一个指针，但实际上却不是。
它只能采用全局常数Softener的值。因此，您可以将其作为参数删除，然后使用全局常量代替。

CheckNode（）

这里，您检查物体是否离节点足够远能够将节点视为单个质心。检查看起来像这样：


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

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



我不太喜欢if语句中的括号，因为它们都是多余的。它看起来可能像这样：

if (pNode->width / distance < _NODE_THRESHOLD || !pNode->HasChildren)


除此之外，我认为您可以使用距离平方而不是距离来保存对sqrt()的调用。因此，换句话说，您可以执行以下操作：

if (nodeWidth / distance < _NODE_THRESHOLD)


而不是： 。CalculateForce（）

这些功能几乎相同。我建议每个Node都应包含一个Body，其中应包含该节点的质心信息。这样，您可以消除CalculateForceNode()，而只需将node->centerOfMassBody传递给CalculateForce()。如果这样做，还可以消除node->TotalMass，node->CenterOfMassx和node->CenterOfMassy，因为所有这些信息都已经在node->centerOfMassBody中。

其他

我只略过了其余的代码。因此，这里还有一些其他小事情：


Coiffecent与Coefficient，Repell与Repel之类的拼写错误。风格的东西）。
您定义了自己的pi常量，但可以只使用M_PI。

回应评论


关于正方形节点，只要您将对角线用作确定阈值的“宽度”，就看不到为什么不能使用矩形节点的任何原因。但是我对Barnes Hut算法不熟悉，因此可能存在一些理论上的理由来说明方形节点更好。
对于从银河系中弹出的粒子，我想边界框的大小会增加，但是我怀疑这会很重要。正如您所说，您始终可以排除异常值。
我相信SetParams()可以使用零向量大小进行调用。想象一下，如果用2个主体调用GenerateChildren()。它创建4个新节点，其中两个为空，其中两个具有1个主体。所有这四个节点都调用了SetParams()，因此其中两个将具有零向量。
我知道a/b与a^2/b^2不同，但我也在与阈值平方进行比较。对于正数a和b，该关系成立：如果a/b < 0.5，则a^2/b^2 < 0.25。
我认为M_PI是在math.h中定义的，但显然并非普遍如此。我习惯了定义了M_PI的Gnu。


评论

\ $ \ begingroup \ $
非常感谢您投入的时间。您使用加速而不用力的想法非常好，我不知道您是怎么看的。对于边界框部分，我无法按质心对其进行拆分，因为每个节点都必须为正方形，并且必须沿中心进行拆分。对于边界框的想法，您是说顶层节点（即整个模拟节点）应尽可能小以包含所有粒子？一个问题是粒子可以从星系中弹出，然后边界框可能比开始时要大！续
\ $ \ endgroup \ $
–基伦·皮尔森（Kieren Pearson）
15年7月7日在1:36
\ $ \ begingroup \ $
我可以隐含全局框的最大大小，但是可以很快达到该大小，可以排除粒子位置值的前5％和botton 5％（除去异常值），但这会在模拟时引起开始在银河系外围的粒子完全不被吸引。
\ $ \ endgroup \ $
–基伦·皮尔森（Kieren Pearson）
15年7月7日在1:39
\ $ \ begingroup \ $
“对我来说不好，因为大小通常为零。尽管您可以用双精度除以零，但是我建议检查顶部附近的大小== 0，并完全避免使用大多数代码。”是的，这确实不会发生，但到目前为止，我将包括更多的安全性。除了GenerateChildren（）之外，SetParam不会被调用，GenerateChildren（）仅在其主体向量的大小大于1时才从SetParam（）进行调用。确实是一个微妙的情况； P
\ $ \ endgroup \ $
–基伦·皮尔森（Kieren Pearson）
15年7月7日在1:45
\ $ \ begingroup \ $
如果（nodeWidthSquared / distanceSquared <_NODE_THRESHOLD_SQUARED）不起作用，因为比率a / b不等于比率a ^ 2 / b ^ 2，并且唯一使用a ^ 2 / b ^ 2的方法是之后，sqrt（a ^ 2 / b ^ 2）仅返回a / b并返回到您开始的位置，但是您需要额外的乘法才能对pNode-> width求平方
\ $ \ endgroup \ $
–基伦·皮尔森（Kieren Pearson）
2015年7月7日在2:07
\ $ \ begingroup \ $
@KierenPearson不，不要发布更新的代码作为编辑。您可以做的是在完成所有修复后，提出一个新问题，然后要求重新审核（如果需要）。如果您不想重新审阅，则可以在问题中添加简短的附录，以描述您要解决的问题，或者在评论中进行。
\ $ \ endgroup \ $
– JS1
2015年7月7日3:00，

#2 楼

这是一段有趣的代码，也是一个写得很好的问题，所以我相信您会得到很多好评。以下是一些可以帮助您改进代码的观察。然后，与它们交互的每一段代码都进入内部并直接操作对象成员，这容易出错并且使维护变得困难。因此，作为一个简单示例，让我们考虑Node。首先，您可以将其设置为Body而不是Body。接下来，class中的某些功能实际上应该是struct的成员功能。举几个明显的例子，main和Body应该是成员函数：作为CalculateForce。

将CalculateForceNode封装为具有成员函数的对象也是一个好主意，而不是将其传递给独立函数。

使用构造函数

Body中的Node函数应改为const的构造函数：未使用的Bodies）要在堆栈而不是堆上创建，避免在该函数中同时使用main.cpp和CreateBody。

使用引用传递复杂的对象

相反，应将以main作为第一个参数的Body声明为采用引用。例如，

// 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);


成为

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)
{}


区别只是Body前面的单个AttractToCenter字符，但是对性能的影响很大，因为它阻止了编译器实际复制矢量。

以上所有建议和一个编译器实现后，该函数将变得非常小而简单：

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


如果您的编译器不太旧，它可能会更简单，如以下建议所示。

使用“ range-for”简化代码

如果编译器完全支持C ++ 11标准，则可以使用“ range-for”简化您的代码码。例如，上面的循环可以这样写：

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


避免使用全局变量

不是常数，应该将其限制在较小的范围内。例如，RepellFromCenter可以简单地移动到new中-不必是全局变量。这样做将减少功能之间的隐藏链接，并帮助您增强界面和设计并减少维护难度。

在包含路径中使用正斜杠

使用此路径您的包含文件仅在Windows下可运行： />优先使用delete值而不是main s

应将std::vector<Body*>和&都声明为pBodies值，而不是IsPaused。这样做可以提供更好的类型安全性，在这种特殊情况下，可以消除混乱的main。

考虑使用更好的随机数生成器

至少支持C ++ 11的编译器，请考虑使用更好的随机数生成器。特别是，您可能要查看const标头中的#define和朋友，而不是_GRAV_CONST。

消除内存泄漏

在许多地方，没有任何相应的_PI调用就调用const float。这是内存泄漏的秘诀。通过向#define添加特定的调用来修复泄漏，或者甚至更好，使用智能指针完全消除static_cast和rand。一项对很多工作有帮助的基本事情是确定哪些对象实际拥有std::uniform_real_distribution和<random>物品。这将帮助您决定是否以及何时删除它们。