File: barnes-hut.js

package info (click to toggle)
axiom 20170501-4
  • links: PTS
  • area: main
  • in suites: buster
  • size: 1,048,504 kB
  • sloc: lisp: 3,600; makefile: 505; cpp: 223; ansic: 138; sh: 96
file content (202 lines) | stat: -rw-r--r-- 8,109 bytes parent folder | download | duplicates (13)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
//
//  barnes-hut.js - version 0.91
//  a graph vizualization toolkit
//
//  Copyright (c) 2011 Samizdat Drafting Co.
//  Physics code derived from springy.js, copyright (c) 2010 Dennis Hotson
// 
//  Permission is hereby granted, free of charge, to any person
//  obtaining a copy of this software and associated documentation
//  files (the "Software"), to deal in the Software without
//  restriction, including without limitation the rights to use,
//  copy, modify, merge, publish, distribute, sublicense, and/or sell
//  copies of the Software, and to permit persons to whom the
//  Software is furnished to do so, subject to the following
//  conditions:
// 
//  The above copyright notice and this permission notice shall be
//  included in all copies or substantial portions of the Software.
// 
//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
//  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
//  OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
//  NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
//  HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
//  WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
//  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
//  OTHER DEALINGS IN THE SOFTWARE.
//

//
//  barnes-hut.js
//
//  implementation of the barnes-hut quadtree algorithm for n-body repulsion
//  http://www.cs.princeton.edu/courses/archive/fall03/cs126/assignments/barnes-hut.html
//
//  Created by Christian Swinehart on 2011-01-14.
//  Copyright (c) 2011 Samizdat Drafting Co. All rights reserved.
//

  var BarnesHutTree = function(){
    var _branches = []
    var _branchCtr = 0
    var _root = null
    var _theta = .5
    
    var that = {
      init:function(topleft, bottomright, theta){
        _theta = theta

        // create a fresh root node for these spatial bounds
        _branchCtr = 0
        _root = that._newBranch()
        _root.origin = topleft
        _root.size = bottomright.subtract(topleft)
      },
      
      insert:function(newParticle){
        // add a particle to the tree, starting at the current _root and working down
        var node = _root
        var queue = [newParticle]

        while (queue.length){
          var particle = queue.shift()
          var p_mass = particle._m || particle.m
          var p_quad = that._whichQuad(particle, node)

          if (node[p_quad]===undefined){
            // slot is empty, just drop this node in and update the mass/c.o.m.
            node[p_quad] = particle
            node.mass += p_mass
            if (node.p){
              node.p = node.p.add(particle.p.multiply(p_mass))
            }else{
              node.p = particle.p.multiply(p_mass)
            }
            
          }else if ('origin' in node[p_quad]){
            // slot conatins a branch node, keep iterating with the branch
            // as our new root
            node.mass += (p_mass)
            if (node.p) node.p = node.p.add(particle.p.multiply(p_mass))
            else node.p = particle.p.multiply(p_mass)
            
            node = node[p_quad]
            queue.unshift(particle)
          }else{
            // slot contains a particle, create a new branch and recurse with
            // both points in the queue now
            var branch_size = node.size.divide(2)
            var branch_origin = new Point(node.origin)
            if (p_quad[0]=='s') branch_origin.y += branch_size.y
            if (p_quad[1]=='e') branch_origin.x += branch_size.x

            // replace the previously particle-occupied quad with a new internal branch node
            var oldParticle = node[p_quad]
            node[p_quad] = that._newBranch()
            node[p_quad].origin = branch_origin
            node[p_quad].size = branch_size
            node.mass = p_mass
            node.p = particle.p.multiply(p_mass)
            node = node[p_quad]

            if (oldParticle.p.x===particle.p.x && oldParticle.p.y===particle.p.y){
              // prevent infinite bisection in the case where two particles
              // have identical coordinates by jostling one of them slightly
              var x_spread = branch_size.x*.08
              var y_spread = branch_size.y*.08
              oldParticle.p.x = Math.min(branch_origin.x+branch_size.x,  
                                         Math.max(branch_origin.x,  
                                                  oldParticle.p.x - x_spread/2 + 
                                                  Math.random()*x_spread))
              oldParticle.p.y = Math.min(branch_origin.y+branch_size.y,  
                                         Math.max(branch_origin.y,  
                                                  oldParticle.p.y - y_spread/2 + 
                                                  Math.random()*y_spread))
            }

            // keep iterating but now having to place both the current particle and the
            // one we just replaced with the branch node
            queue.push(oldParticle)
            queue.unshift(particle)
          }

        }

      },

      applyForces:function(particle, repulsion){
        // find all particles/branch nodes this particle interacts with and apply
        // the specified repulsion to the particle
        var queue = [_root]
        while (queue.length){
          node = queue.shift()
          if (node===undefined) continue
          if (particle===node) continue
          
          if ('f' in node){
            // this is a particle leafnode, so just apply the force directly
            var d = particle.p.subtract(node.p);
            var distance = Math.max(1.0, d.magnitude());
            var direction = ((d.magnitude()>0) ? d : Point.random(1)).normalize()
            particle.applyForce(direction.multiply(repulsion*(node._m||node.m))
                                      .divide(distance * distance) );
          }else{
            // it's a branch node so decide if it's cluster-y and distant enough
            // to summarize as a single point. if it's too complex, open it and deal
            // with its quadrants in turn
            var dist = particle.p.subtract(node.p.divide(node.mass)).magnitude()
            var size = Math.sqrt(node.size.x * node.size.y)
            if (size/dist > _theta){ // i.e., s/d > Θ
              // open the quad and recurse
              queue.push(node.ne)
              queue.push(node.nw)
              queue.push(node.se)
              queue.push(node.sw)
            }else{
              // treat the quad as a single body
              var d = particle.p.subtract(node.p.divide(node.mass));
              var distance = Math.max(1.0, d.magnitude());
              var direction = ((d.magnitude()>0) ? d : Point.random(1)).normalize()
              particle.applyForce(direction.multiply(repulsion*(node.mass))
                                           .divide(distance * distance) );
            }
          }
        }
      },
      
      _whichQuad:function(particle, node){
        // sort the particle into one of the quadrants of this node
        if (particle.p.exploded()) return null
        var particle_p = particle.p.subtract(node.origin)
        var halfsize = node.size.divide(2)
        if (particle_p.y < halfsize.y){
          if (particle_p.x < halfsize.x) return 'nw'
          else return 'ne'
        }else{
          if (particle_p.x < halfsize.x) return 'sw'
          else return 'se'
        }
      },
      
      _newBranch:function(){
        // to prevent a gc horrorshow, recycle the tree nodes between iterations
        if (_branches[_branchCtr]){
          var branch = _branches[_branchCtr]
          branch.ne = branch.nw = branch.se = branch.sw = undefined
          branch.mass = 0
          delete branch.p
        }else{
          branch = {origin:null, size:null, 
                    nw:undefined, ne:undefined, sw:undefined, se:undefined, mass:0}
          _branches[_branchCtr] = branch
        }

        _branchCtr++
        return branch
      }
    }
    
    return that
  }