Barnes-Hut  算法

参考 http://arborjs.org/docs/barnes-hut

该算法对区域进行4分割。直到区域中只包含1个或者0个元素。

如下图

通过分割构造出如下树。

 

 

 

递归构造树的算法

1 bool Tree::buildTree(NbodyNode *tree, complex
     
       start, complex
      
        end, vector
       
         &plants) { 2     if (tree == nullptr) 3         return false; 4     complex
        
          check = end - start; 5     if (check.real() <= 0 || check.imag() <= 0) { 6         printf("check failed\n"); 7         return false; 8     } 9     if (plants.size() == 1) {10         tree->body() = plants.front();11         tree->isPlant() = true;12         this->plants.push_back(&tree->body());13         return true;14     }15     vector
         
           wrapers[4];16     int centerX = (start.real() + end.real()) / 2;17     int centerY = (start.imag() + end.imag()) / 2;18     complex
          
            center = complex
           
            (centerX, centerY);19 complex
            
              sub = complex
             
              ();20 for (vector
              
               ::iterator i = plants.begin(); i != plants.end(); i++) {21 sub = i->location() - center;22 if (sub.real() <= 0 && sub.imag() <= 0) {23 wrapers[0].push_back(*i);24 } else if (sub.real() < 0 && sub.imag() > 0) {25 wrapers[2].push_back(*i);26 } else if (sub.real() > 0 && sub.imag() < 0) {27 wrapers[1].push_back(*i);28 } else if (sub.real() >= 0 && sub.imag() >= 0) {29 wrapers[3].push_back(*i);30 }31 }32 int width = tree->width() / 4;33 tree->body() = Plant();34 bool ret = true;35 if (wrapers[0].size() > 0) {36 tree->leftTop() = new NbodyNode(width);37 ret = ret && buildTree(tree->leftTop(), start, center, wrapers[0]);38 tree->body() = tree->body() + tree->leftTop()->body();39 }40 if (wrapers[1].size() > 0) {41 tree->rightTop() = new NbodyNode(width);42 ret = ret && buildTree(tree->rightTop(), complex
               
                (start.real() + centerX, start.imag()),43 complex
                
                 (end.real(), centerY), wrapers[1]);44 tree->body() = tree->body() + tree->rightTop()->body();45 }46 if (wrapers[2].size() > 0) {47 tree->leftButtom() = new NbodyNode(width);48 ret = ret && buildTree(tree->leftButtom(), complex
                 
                  (start.real(), centerY),49 complex
                  
                   (centerX, end.imag()), wrapers[2]);50 tree->body() = tree->body() + tree->leftButtom()->body();51 }52 if (wrapers[3].size() > 0) {53 tree->rightButtom() = new NbodyNode(width);54 ret = ret && buildTree(tree->rightButtom(), center, end, wrapers[3]);55 tree->body() = tree->body() + tree->rightButtom()->body();56 }57 return ret;58 }
                  
                 
                
               
              
             
            
           
          
         
        
       
      
     

树中每一个非NULL节点保存该区域中星体的等效值。

若是星体，保存本身。若不是，保存该区域中的等效星体。

即

星体1 质量M1 位置（x1,y1）星体2 质量M2 位置（x2,y2）

等效星体 质量M = M1+M2 位置（x = (x1*M1+x*M2)/M, y = (y1*M1+y2*M2)/M）;

如下图

s 为该区域的宽度

d 为A星体到蓝色区域等效星体的距离

若 d/s < θ

则该区域可以被等效，否则计算该区域的子区域。

若区域本身是一个星体，则直接计算该星体对A的万有引力。不用计算 d/s