nbodyforce.as

用于flash/flex的 as3的 2D图形图像图表的动态生成

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