Barnes-Hut 算法
参考 http://arborjs.org/docs/barnes-hut
该算法对区域进行4分割。直到区域中只包含1个或者0个元素。
如下图
通过分割构造出如下树。
递归构造树的算法
1 bool Tree::buildTree(NbodyNode *tree, complex<double> start, complex<double> end, vector<Plant> &plants) { 2 if (tree == nullptr) 3 return false; 4 complex<double> 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<Plant> wrapers[4]; 16 int centerX = (start.real() + end.real()) / 2; 17 int centerY = (start.imag() + end.imag()) / 2; 18 complex<double> center = complex<double>(centerX, centerY); 19 complex<double> sub = complex<double>(); 20 for (vector<Plant>::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<double>(start.real() + centerX, start.imag()), 43 complex<double>(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<double>(start.real(), centerY), 49 complex<double>(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