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);
    }
}