```  1 /**
2  * Constructs a new charge force, with an optional charge constant. The charge
3  * constant can be negative for repulsion (e.g., particles with electrical
4  * charge of equal sign), or positive for attraction (e.g., massive particles
5  * with mutual gravity). The default charge constant is -40.
6  *
7  * @class An n-body force, as defined by Coulomb's law or Newton's law of
8  * gravitation, inversely proportional to the square of the distance between
9  * particles. Note that the force is independent of the <i>mass</i> of the
10  * associated particles, and that the particles do not have charges of varying
11  * magnitude; instead, the attraction or repulsion of all particles is globally
12  * specified as the charge {@link #constant}.
13  *
14  * <p>This particular implementation uses the Barnes-Hut algorithm. For details,
15  * see <a
16  * href="http://www.nature.com/nature/journal/v324/n6096/abs/324446a0.html">"A
17  * hierarchical O(N log N) force-calculation algorithm"</a>, J. Barnes &
18  * P. Hut, <i>Nature</i> 1986.
19  *
20  * @name pv.Force.charge
21  * @param {number} [k] the charge constant.
22  */
23 pv.Force.charge = function(k) {
24   var min = 2, // minimum distance at which to observe forces
25       min1 = 1 / min,
26       max = 500, // maximum distance at which to observe forces
27       max1 = 1 / max,
28       theta = .9, // Barnes-Hut theta approximation constant
29       force = {};
30
31   if (!arguments.length) k = -40; // default charge constant (repulsion)
32
33   /**
34    * Sets or gets the charge constant. If an argument is specified, it is the
35    * new charge constant. The charge constant can be negative for repulsion
36    * (e.g., particles with electrical charge of equal sign), or positive for
37    * attraction (e.g., massive particles with mutual gravity). The default
38    * charge constant is -40.
39    *
40    * @function
41    * @name pv.Force.charge.prototype.constant
42    * @param {number} x the charge constant.
43    * @returns {pv.Force.charge} this.
44    */
45   force.constant = function(x) {
46     if (arguments.length) {
47       k = Number(x);
48       return force;
49     }
50     return k;
51   };
52
53   /**
54    * Sets or gets the domain; specifies the minimum and maximum domain within
55    * which charge forces are applied. A minimum distance threshold avoids
56    * applying forces that are two strong (due to granularity of the simulation's
57    * numeric integration). A maximum distance threshold improves performance by
58    * skipping force calculations for particles that are far apart.
59    *
60    * <p>The default domain is [2, 500].
61    *
62    * @function
63    * @name pv.Force.charge.prototype.domain
64    * @param {number} a
65    * @param {number} b
66    * @returns {pv.Force.charge} this.
67    */
68   force.domain = function(a, b) {
69     if (arguments.length) {
70       min = Number(a);
71       min1 = 1 / min;
72       max = Number(b);
73       max1 = 1 / max;
74       return force;
75     }
76     return [min, max];
77   };
78
79   /**
80    * Sets or gets the Barnes-Hut approximation factor. The Barnes-Hut
81    * approximation criterion is the ratio of the size of the quadtree node to
82    * the distance from the point to the node's center of mass is beneath some
83    * threshold.
84    *
85    * @function
86    * @name pv.Force.charge.prototype.theta
87    * @param {number} x the new Barnes-Hut approximation factor.
88    * @returns {pv.Force.charge} this.
89    */
90   force.theta = function(x) {
91     if (arguments.length) {
92       theta = Number(x);
93       return force;
94     }
95     return theta;
96   };
97
98   /**
99    * @ignore Recursively computes the center of charge for each node in the
100    * quadtree. This is equivalent to the center of mass, assuming that all
101    * particles have unit weight.
102    */
103   function accumulate(n) {
104     var cx = 0, cy = 0;
105     n.cn = 0;
106     function accumulateChild(c) {
107       accumulate(c);
108       n.cn += c.cn;
109       cx += c.cn * c.cx;
110       cy += c.cn * c.cy;
111     }
112     if (!n.leaf) {
113       if (n.c1) accumulateChild(n.c1);
114       if (n.c2) accumulateChild(n.c2);
115       if (n.c3) accumulateChild(n.c3);
116       if (n.c4) accumulateChild(n.c4);
117     }
118     if (n.p) {
119       n.cn += k;
120       cx += k * n.p.x;
121       cy += k * n.p.y;
122     }
123     n.cx = cx / n.cn;
124     n.cy = cy / n.cn;
125   }
126
127   /**
128    * @ignore Recursively computes forces on the given particle using the given
129    * quadtree node. The Barnes-Hut approximation criterion is the ratio of the
130    * size of the quadtree node to the distance from the point to the node's
131    * center of mass is beneath some threshold.
132    */
133   function forces(n, p, x1, y1, x2, y2) {
134     var dx = n.cx - p.x,
135         dy = n.cy - p.y,
136         dn = 1 / Math.sqrt(dx * dx + dy * dy);
137
138     /* Barnes-Hut criterion. */
139     if ((n.leaf && (n.p != p)) || ((x2 - x1) * dn < theta)) {
140       if (dn < max1) return;
141       if (dn > min1) dn = min1;
142       var kc = n.cn * dn * dn * dn,
143           fx = dx * kc,
144           fy = dy * kc;
145       p.fx += fx;
146       p.fy += fy;
147     } else if (!n.leaf) {
148       var sx = (x1 + x2) * .5, sy = (y1 + y2) * .5;
149       if (n.c1) forces(n.c1, p, x1, y1, sx, sy);
150       if (n.c2) forces(n.c2, p, sx, y1, x2, sy);
151       if (n.c3) forces(n.c3, p, x1, sy, sx, y2);
152       if (n.c4) forces(n.c4, p, sx, sy, x2, y2);
153       if (dn < max1) return;
154       if (dn > min1) dn = min1;
155       if (n.p && (n.p != p)) {
156         var kc = k * dn * dn * dn,
157             fx = dx * kc,
158             fy = dy * kc;
159         p.fx += fx;
160         p.fy += fy;
161       }
162     }
163   }
164
165   /**
166    * Applies this force to the specified particles. The force is applied between
167    * all pairs of particles within the domain, using the specified quadtree to
168    * accelerate n-body force calculation using the Barnes-Hut approximation
169    * criterion.
170    *
171    * @function
172    * @name pv.Force.charge.prototype.apply
173    * @param {pv.Particle} particles particles to which to apply this force.