    package flare.physics
    {
        /**
         * Force simulating an N-Body force of charged particles with pairwise
         * interaction, such as gravity or electrical charge. This class uses a
         * quad-tree structure to aggregate charge values and optimize computation.
         * The force function is a standard inverse-square law (though in this case
         * approximated due to optimization): <code>F = G * m1 * m2 / d^2</code>,
         * where G is a constant (e.g., gravitational constant), m1 and m2 are the
         * masses (charge) of the particles, and d is the distance between them.
         *
         * <p>The algorithm used is that of J. Barnes and P. Hut, in their research
         * paper <i>A Hierarchical  O(n log n) force calculation algorithm</i>, Nature,
         * v.324, December 1986. For more details on the algorithm, see one of
         * the following links:
         * <ul>
         *   <li><a href="http://www.cs.berkeley.edu/~demmel/cs267/lecture26/lecture26.html">James Demmel's UC Berkeley lecture notes</a>
         *   <li><a href="http://www.physics.gmu.edu/~large/lr_forces/desc/bh/bhdesc.html">Description of the Barnes-Hut algorithm</a>
         *   <li><a href="http://www.ifa.hawaii.edu/~barnes/treecode/treeguide.html">Joshua Barnes' implementation</a>
         * </ul></p>
         */
        public class NBodyForce implements IForce
        {
            private var _g:Number;     // gravitational constant
            private var _t:Number;     // barnes-hut theta
            private var _max:Number;   // max effective distance
            private var _min:Number;   // min effective distance
            private var _eps:Number;   // epsilon for determining 'same' location

            private var _x1:Number, _y1:Number, _x2:Number, _y2:Number;
            private var _root:QuadTreeNode;

            /** The gravitational constant to use.
             *  Negative values produce a repulsive force. */
            public function get gravitation():Number { return _g; }
            public function set gravitation(g:Number):void { _g = g; }

            /** The maximum distance over which forces are exerted.
             *  Any greater distances will be ignored. */
            public function get maxDistance():Number { return _max; }
            public function set maxDistance(d:Number):void { _max = d; }

            /** The minumum effective distance over which forces are exerted.
             *     Any lesser distances will be treated as the minimum. */
            public function get minDistance():Number { return _min; }
            public function set minDistance(d:Number):void { _min = d; }

            // --------------------------------------------------------------------

            /**
             * Creates a new NBodyForce with given parameters.
             * @param g the gravitational constant to use.
             *  Negative values produce a repulsive force.
             * @param maxd a maximum distance over which the force should operate.
             *  Particles separated by more than this distance will not interact.
             * @param mind the minimum distance over which the force should operate.
             *  Particles closer than this distance will interact as if they were
             *  the minimum distance apart. This helps avoid extreme forces.
             *  Helpful when particles are very close together.
             * @param eps an epsilon values for determining a minimum distance
             *  between particles
             * @param t the theta parameter for the Barnes-Hut approximation.
             *  Determines the level of approximation (default value if 0.9).
             */
            public function NBodyForce(g:Number=-1, max:Number=200, min:Number=2,
                                       eps:Number=0.01, t:Number=0.9)
            {
                _g = g;
                _max = max;
                _min = min;
                _eps = eps;
                _t = t;
                _root = QuadTreeNode.node();
            }

            /**
             * Applies this force to a simulation.
             * @param sim the Simulation to apply the force to
             */
            public function apply(sim:Simulation):void
            {
                if (_g == 0) return;

                // clear the quadtree
                clear(_root); _root = QuadTreeNode.node();

                // get the tree bounds
                bounds(sim);

                // populate the tree
                for (var i:uint = 0; i<sim.particles.length; ++i) {
                    insert(sim.particles[i], _root, _x1, _y1, _x2, _y2);
                }

                // traverse tree to compute mass
                accumulate(_root);

                // calculate forces on each particle
                for (i=0; i<sim.particles.length; ++i) {
                    forces(sim.particles[i], _root, _x1, _y1, _x2, _y2);
                }
            }

            private function accumulate(n:QuadTreeNode):void {
                var xc:Number = 0, yc:Number = 0;
                n.mass = 0;

                // accumulate childrens' mass
                var recurse:Function = function(c:QuadTreeNode):void {
                    if (c == null) return;
                    accumulate(c);
                    n.mass += c.mass;
                    xc += c.mass * c.cx;
                    yc += c.mass * c.cy;
                }
                if (n.hasChildren) {
                    recurse(n.c1); recurse(n.c2); recurse(n.c3); recurse(n.c4);
                }

                // accumulate own mass
                if (n.p != null) {
                    n.mass += n.p.mass;
                    xc += n.p.mass * n.p.x;
                    yc += n.p.mass * n.p.y;
                }
                n.cx = xc / n.mass;
                n.cy = yc / n.mass;
            }

            private function forces(p:Particle, n:QuadTreeNode,
                x1:Number, y1:Number, x2:Number, y2:Number):void
            {
                var f:Number = 0;
                var dx:Number = n.cx - p.x;
                var dy:Number = n.cy - p.y;
                var dd:Number = Math.sqrt(dx*dx + dy*dy);
                var max:Boolean = _max > 0 && dd > _max;
                if (dd==0) { // add direction when needed
                    dx = _eps * (0.5-Math.random());
                    dy = _eps * (0.5-Math.random());
                }

                // the Barnes-Hut approximation criteria is if the ratio of the
                // size of the quadtree box to the distance between the point and
                // the box's center of mass is beneath some threshold theta.
                if ( (!n.hasChildren && n.p != p) || ((x2-x1)/dd < _t) )
                {
                    if ( max ) return;
                    // either only 1 particle or we meet criteria
                    // for Barnes-Hut approximation, so calc force
                    dd = dd<_min ? _min : dd;
                    f = _g * p.mass * n.mass / (dd*dd*dd)
                    p.fx += f*dx; p.fy += f*dy;
                }
                else if ( n.hasChildren )
                {
                    // recurse for more accurate calculation
                    var sx:Number = (x1+x2)/2
                    var sy:Number = (y1+y2)/2;

                    if (n.c1) forces(p, n.c1, x1, y1, sx, sy);
                    if (n.c2) forces(p, n.c2, sx, y1, x2, sy);
                    if (n.c3) forces(p, n.c3, x1, sy, sx, y2);
                    if (n.c4) forces(p, n.c4, sx, sy, x2, y2);

                    if ( max ) return;
                    if ( n.p != null && n.p != p ) {
                        dd = dd<_min ? _min : dd;
                        f = _g * p.mass * n.p.mass / (dd*dd*dd);
                        p.fx += f*dx; p.fy += f*dy;
                    }
                }
            }

            // -- Helpers ---------------------------------------------------------

            private function insert(p:Particle, n:QuadTreeNode,
                x1:Number, y1:Number, x2:Number, y2:Number):void
            {
                // ignore particles with NaN coordinates
                if (isNaN(p.x) || isNaN(p.y)) return;

                // try to insert particle p at node n in the quadtree
                // by construction, each leaf will contain either 1 or 0 particles
                if ( n.hasChildren ) {
                    // n contains more than 1 particle
                    insertHelper(p,n,x1,y1,x2,y2);
                } else if ( n.p != null ) {
                    // n contains 1 particle
                    if ( isSameLocation(n.p, p) ) {
                        // recurse
                        insertHelper(p,n,x1,y1,x2,y2);
                    } else {
                        // divide
                        var v:Particle = n.p; n.p = null;
                        insertHelper(v,n,x1,y1,x2,y2);
                        insertHelper(p,n,x1,y1,x2,y2);
                    }
                } else {
                    // n is empty, add p as leaf
                    n.p = p;
                }
            }

            private function insertHelper(p:Particle, n:QuadTreeNode,
                x1:Number, y1:Number, x2:Number, y2:Number):void
            {
                // determine split
                var sx:Number = (x1+x2)/2;
                var sy:Number = (y1+y2)/2;
                var c:uint = (p.x >= sx ? 1 : 0) + (p.y >= sy ? 2 : 0);

                // update bounds
                if (c==1 || c==3) x1 = sx; else x2 = sx;
                if (c>1) y1 = sy; else y2 = sy;

                // update children
                var cn:QuadTreeNode;
                if (c == 0) {
                    if (n.c1==null) n.c1 = QuadTreeNode.node();
                    cn = n.c1;
                } else if (c == 1) {
                    if (n.c2==null) n.c2 = QuadTreeNode.node();
                    cn = n.c2;
                } else if (c == 2) {
                    if (n.c3==null) n.c3 = QuadTreeNode.node();
                    cn = n.c3;
                } else {
                    if (n.c4==null) n.c4 = QuadTreeNode.node();
                    cn = n.c4;
                }
                n.hasChildren = true;
                insert(p,cn,x1,y1,x2,y2);
            }

            private function clear(n:QuadTreeNode):void
            {
                if (n.c1 != null) clear(n.c1);
                if (n.c2 != null) clear(n.c2);
                if (n.c3 != null) clear(n.c3);
                if (n.c4 != null) clear(n.c4);
                QuadTreeNode.reclaim(n);
            }

            private function bounds(sim:Simulation):void
            {
                var p:Particle, dx:Number, dy:Number;
                _x1 = _y1 = Number.MAX_VALUE;
                _x2 = _y2 = Number.MIN_VALUE;

                // get bounding box
                for (var i:uint = 0; i<sim.particles.length; ++i) {
                    p = sim.particles[i] as Particle;
                    if (p.x < _x1) _x1 = p.x;
                    if (p.y < _y1) _y1 = p.y;
                    if (p.x > _x2) _x2 = p.x;
                    if (p.y > _y2) _y2 = p.y;
                }

                // square the box
                dx = _x2 - _x1;
                dy = _y2 - _y1;
                if (dx > dy) {
                    _y2 = _y1 + dx;
                } else {
                    _x2 = _x1 + dy;
                }
            }

            private function isSameLocation(p1:Particle, p2:Particle):Boolean {
                return (Math.abs(p1.x - p2.x) < _eps &&
                        Math.abs(p1.y - p2.y) < _eps);
            }

        } // end of class NBodyForce
    }

    // -- Helper QuadTreeNode class -----------------------------------------------

    import flare.physics.Particle;

    class QuadTreeNode
    {
        public var mass:Number = 0;
        public var cx:Number = 0;
        public var cy:Number = 0;
        public var p:Particle = null;
        public var c1:QuadTreeNode = null;
        public var c2:QuadTreeNode = null;
        public var c3:QuadTreeNode = null;
        public var c4:QuadTreeNode = null;
        public var hasChildren:Boolean = false;

        // -- Factory ---------------------------------------------------------

        private static var _nodes:Array = new Array();

        public static function node():QuadTreeNode {
            var n:QuadTreeNode;
            if (_nodes.length > 0) {
                n = QuadTreeNode(_nodes.pop());
            } else {
                n = new QuadTreeNode();
            }
            return n;
        }

        public static function reclaim(n:QuadTreeNode):void {
            n.mass = n.cx = n.cy = 0;
            n.p = null;
            n.hasChildren = false;
            n.c1 = n.c2 = n.c3 = n.c4 = null;
            _nodes.push(n);
        }
    }