/** 
* PHY2D physics 2D engine
* C# and AS3 code by Paul Firth http://www.wildbunny.co.uk/
* JavaScript adaptation by @ge1doot http://codepen.io/ge1doot/
*/ 

 
! function ()
{
	"use strict";
	
	/** 
	* init screen and pointer
	*/ 

	var screen = ge1doot.screen.init("screen", function ()
	{
		PHY2D.zoom = 1;
		PHY2D.deleteStatic();
		PHY2D.rectangle(screen.width / 2, screen.height + 50, screen.width, 100, 0, 0);
		PHY2D.rectangle(screen.width / 2, -3000, screen.width, 100, 0, 0);
		PHY2D.rectangle(-50, -1500 + screen.height, 100, 3000 + screen.height, 0, 0);
		PHY2D.rectangle(screen.width + 50, -1500 + screen.height, 100, 3000 + screen.height, 0, 0);
		PHY2D.zoom = screen.width / 1660;
	}, false);
	
	var ctx = screen.ctx;
	
	if (!ctx.setLineDash)
	{
		ctx.setLineDash = function () {}
	}
	
	var pointer = screen.pointer.init(
	{
		up: function ()
		{
			PHY2D.stopManipulate();
		}
	});

	/** 
	* PHY2D engine
	*/ 
		
	var PHY2D = {
		kGravity: 5,
		kTimeStep: 1 / 60,
		kFriction: 0.5,
		kMotionAABB: true,
		kImgPath: "http://www.dhteumeuleu.com/images/"
	}
	
	/** 
	* Vector 2D library
	*/ 
	
	PHY2D.Vector = function (x, y)
	{
		this.x = x || 0.0;
		this.y = y || 0.0;
	}
	
	PHY2D.Vector.prototype = {
	
	
		set: function (x, y)
		{
			this.x = x;
			this.y = y;
			return this;
		},
		
		
		dot: function (v)
		{
			return this.x * v.x + this.y * v.y;
		},
		
		
		lenSqr: function ()
		{
			return this.x * this.x + this.y * this.y;
		},
		
		
		len: function ()
		{
			return Math.sqrt(this.x * this.x + this.y * this.y);
		},
		
		
		dist: function (v)
		{
			var dx = v.x - this.x;
			var dy = v.y - this.y;
			return Math.sqrt(dx * dx + dy * dy);
		},
		
		
		transform: function (v, m)
		{
			this.x = m.cos * v.x - m.sin * v.y + m.tx;
			this.y = m.sin * v.x + m.cos * v.y + m.ty;
			return this;
		},
		
		
		rotate: function (v, m)
		{
			this.x = m.cos * v.x - m.sin * v.y;
			this.y = m.sin * v.x + m.cos * v.y;
			return this;
		},
		
		
		normal: function (a, b)
		{
			var x = a.x - b.x,
				y = a.y - b.y,
				len = Math.sqrt(x * x + y * y);
			this.x = -y / len;
			this.y =  x / len;
			return this;
		},
		
		
		project: function (a, b, n)
		{
			var x = a.x - b.x,
				y = a.y - b.y,
				len = Math.sqrt(x * x + y * y);
			return (-y / len) * n.x + (x / len) * n.y;
		},
		
		
		addScale: function (v1, v2, s)
		{
			this.x = v1.x + (v2.x * s);
			this.y = v1.y + (v2.y * s);
			return this;
		},
		
		
		subScale: function (v1, v2, s)
		{
			this.x = v1.x - (v2.x * s);
			this.y = v1.y - (v2.y * s);
			return this;
		},
		
		
		add: function (v1, v2)
		{
			this.x = v1.x + v2.x;
			this.y = v1.y + v2.y;
			return this;
		},
		
		
		sub: function (v1, v2)
		{
			this.x = v1.x - v2.x;
			this.y = v1.y - v2.y;
			return this;
		},
		
		
		scale: function (v1, s)
		{
			this.x = v1.x * s;
			this.y = v1.y * s;
			return this;
		},
		
		
		perpSelf: function ()
		{
			var x = this.x;
			this.x = -this.y;
			this.y = x;
			return this;
		},
		
		
		perp: function (v1)
		{
			this.x = -v1.y;
			this.y =  v1.x;
			return this;
		},
		
		
		inv: function (v1)
		{
			this.x = -v1.x;
			this.y = -v1.y;
			return this;
		},
		
		
		copy: function (v1)
		{
			this.x = v1.x;
			this.y = v1.y;
			return this;
		},
		
		
		unit: function ()
		{
			var invLen = 1.0 / Math.sqrt(this.x * this.x + this.y * this.y);
			this.x *= invLen;
			this.y *= invLen;
			return this;
		},
		
		
		unitLen: function (len)
		{
			var invLen = 1.0 / len;
			this.x *= invLen;
			this.y *= invLen;
			return this;
		},
		
		
		clamp: function (v, min, max)
		{
			if (v > max) v = max;
			else if (v < min) v = min;
			return v;
		},
		
		
		rotateIntoSpaceOf: function (a, m)
		{
			var dx = -a.x,
				dy = -a.y;
			this.x = dx *  m.cos + dy * m.sin;
			this.y = dx * -m.sin + dy * m.cos;
			return this;
		},
		
		
		projectPointOntoEdge: function (p, e0, e1)
		{
			var ex = e1.x - e0.x;
			var ey = e1.y - e0.y;
			var t = (ex * (p.x - e0.x) + ey * (p.y - e0.y)) / (ex * ex + ey * ey);
			if (t < 0) t = 0; else if (t > 1) t = 1;
			this.x = e0.x + (ex * t);
			this.y = e0.y + (ey * t);
			return this;
		},
		
		
		/** 
		* pseudo SIMD array vectors
		*/ 
		
		array: function (n, values)
		{
			var array = new Array(n);
			array.min = new PHY2D.Vector();
			array.max = new PHY2D.Vector();
			for (var i = 0; i < n; i++)
			{
				array[i] = new PHY2D.Vector(
					values ? values[i * 2 + 0] : 0.0,
					values ? values[i * 2 + 1] : 0.0
				);
			}
			
			
			array.transform = function (v, m)
			{
				for (var i = 0, len = this.length; i < len; i++)
				{
					var vi = v[i],
						elem = this[i];
					var x = m.cos * vi.x - m.sin * vi.y + m.tx;
					var y = m.sin * vi.x + m.cos * vi.y + m.ty;
					if (x < this.min.x) this.min.x = x;
					if (y < this.min.y) this.min.y = y;
					if (x > this.max.x) this.max.x = x;
					if (y > this.max.y) this.max.y = y;
					elem.x = x;
					elem.y = y;
				}
				return this;
			}
			
			
			array.rotate = function (v, m)
			{
				for (var i = 0, len = this.length; i < len; i++)
				{
					var vi = v[i],
						elem = this[i];
					elem.x = m.cos * vi.x - m.sin * vi.y;
					elem.y = m.sin * vi.x + m.cos * vi.y;
				}
				return this;
			}
			
			
			array.resetMinmax = function ()
			{
				this.min.x = 100000.0;
				this.min.y = 100000.0;
				this.max.x = -100000.0;
				this.max.y = -100000.0;
			}
			
			
			array.normal = function (points)
			{
				for (var i = 0; i < this.length; i++)
				{
					this[i].normal(
						points[(i + 1) % this.length],
						points[i]
					);
				}
				return this;
			}
			
			
			return array;
		}
	}
	

	/** 
	* Matrix container
	*/ 
	
	PHY2D.Matrix = function ()
	{
		this.cos = 0.0;
		this.sin = 0.0;
		this.tx  = 0.0;
		this.ty  = 0.0;
	}
	
	
	PHY2D.Matrix.prototype.set = function (x, y, a)
	{
		this.cos = Math.cos(a);
		this.sin = Math.sin(a);
		this.tx  = x;
		this.ty  = y;
		return this;
	}
	
	/** 
	* temporary working vectors
	*/ 
	
	PHY2D.v0 = new PHY2D.Vector();
	PHY2D.v1 = new PHY2D.Vector();
	PHY2D.v2 = new PHY2D.Vector();
	PHY2D.v3 = new PHY2D.Vector();
	PHY2D.v4 = new PHY2D.Vector();
	PHY2D.v5 = new PHY2D.Vector();
	PHY2D.z0 = new PHY2D.Vector();
	
	/** 
	* objects and contacts list arrays
	*/ 
	
	PHY2D.objects = [];
	PHY2D.manipulate = null;
	PHY2D.joints = [];
	PHY2D.contacts = [];
	PHY2D.contacts.index = 0;
	PHY2D.contacts.create = function (A, B, pa, pb, dist, normal)
	{
		if (!this[this.index]) this[this.index] = new PHY2D.Contact();
		this[this.index++].set(A, B, pa, pb, dist, normal);
	}
	
	
	/** 
	* AABB container
	*/ 

	
	PHY2D.AABB = function ()
	{
		this.x = 0.0;
		this.y = 0.0;
		this.w = 0.0;
		this.h = 0.0;
	}
	
	/** 
	* Rigidbody base PROTOTYPE definition
	*/ 
	
	PHY2D.RigidBody = {};
	
	
	/** 
	* Init new Rigidbody object
	*/ 
	
	PHY2D.RigidBody.initBase = function (x, y, w, h, mass, angle, img)
	{
		this.img = null;
		if (img)
		{
			var image = new Image();
			image.src = PHY2D.kImgPath + img.src;
			this.img = {
				w: (img.w ? img.w * PHY2D.zoom : w),
				h: (img.h ? img.h * PHY2D.zoom : h),
				x: (img.x ? img.x * PHY2D.zoom : 0),
				y: (img.y ? img.y * PHY2D.zoom : 0),
				elem: image
			}
		}
		this.pos       = new PHY2D.Vector(x, y);
		this.vel       = new PHY2D.Vector();
		this.ang       = angle;
		this.angVel    = 0.0;
		this.invMass   = mass ? 1 / mass : 0;
		this.matrix    = new PHY2D.Matrix().set(x, y, angle);
		this.nextFrame = new PHY2D.Matrix();
		this.aabb      = new PHY2D.AABB();
		this.c1        = new PHY2D.Vector();
		this.c0        = new PHY2D.Vector();
		this.normal    = new PHY2D.Vector();
		this.static    = false;
	}
	
	/** 
	* draw rigid body image
	*/ 
	
	PHY2D.RigidBody.draw = function ()
	{
		var img = this.img;
		if (img)
		{
			ctx.save();
			ctx.translate(this.pos.x, this.pos.y);
			ctx.rotate(this.ang);
			ctx.drawImage(img.elem, -img.w * 0.5 + img.x, -img.h * 0.5 + img.y, img.w, img.h);
			ctx.restore();
		}
	}
	
	/** 
	* Rigidbody integration
	*/ 

	PHY2D.RigidBody.integrate = function ()
	{
		// gravity
		if (this.invMass > 0) this.vel.y += PHY2D.kGravity;
		// update position
		var vx = this.vel.x * PHY2D.kTimeStep;
		var vy = this.vel.y * PHY2D.kTimeStep;
		var va = this.angVel * PHY2D.kTimeStep;
		this.pos.x += vx;
		this.pos.y += vy;
		this.ang   += va;
		// update matrix
		this.matrix.set(this.pos.x, this.pos.y, this.ang);
		if (PHY2D.kMotionAABB)
		{
			this.nextFrame.set(this.pos.x + vx, this.pos.y + vy, this.ang + va);
		}
		// compute motion AABB
		if (!this.static) this.motionAABB();
		else
		{
			if (this.invMass === 0)
			{
				this.static = true;
				this.motionAABB();
			}
		}
		// manipulate
		if (pointer.active && !PHY2D.manipulate && this.invMass)
		{
			if (this.isPointInPoly(pointer.pos.x, pointer.pos.y))
			{
				PHY2D.manipulate = new PHY2D.Manipulate(this);
				PHY2D.manipulate.solve();
			}
		}
	}	
	

	/** 
	* Polygon PROTOTYPE definition
	* extends RigidBody base prototype
	*/ 
	
	PHY2D.Polygon = Object.create(PHY2D.RigidBody);
	
	/** 
	* Init new Polygon object
	*/ 
	
	PHY2D.Polygon.init = function (x, y, w, h, vertices, mass, angle, img)
	{
		this.type = "poly";
		this.initBase(x, y, w, h, mass, angle, img);
		this.length = (vertices.length / 2) | 0;
		// vertices
		this.localSpacePoints  = new PHY2D.Vector().array(this.length, vertices);
		this.localSpaceNormals = new PHY2D.Vector().array(this.length).normal(this.localSpacePoints);
		this.worldSpaceNormals = new PHY2D.Vector().array(this.length);
		this.worldSpacePoints  = new PHY2D.Vector().array(this.length);
		// calculate inverse inertia tensor
		this.invI = (this.invMass > 0) ? 1 / ((1 / this.invMass) * (w * w + h * h) / 3) : 0
		return this;
	}

	/** 
	* AABB motion
	*/ 
	
	PHY2D.Polygon.motionAABB = function ()
	{
		this.worldSpacePoints.resetMinmax();
		PHY2D.kMotionAABB && this.worldSpacePoints.transform(this.localSpacePoints, this.nextFrame);
		this.worldSpacePoints.transform(this.localSpacePoints, this.matrix);
		this.worldSpaceNormals.rotate(this.localSpaceNormals, this.matrix);
		var min = this.worldSpacePoints.min;
		var max = this.worldSpacePoints.max;
		this.aabb.x = (min.x + max.x) * 0.5;
		this.aabb.y = (min.y + max.y) * 0.5;
		this.aabb.w = (max.x - min.x) * 0.5;
		this.aabb.h = (max.y - min.y) * 0.5;
	}
	
	/** 
	* is point in polygon ?
	*/ 
	
	PHY2D.Polygon.isPointInPoly = function (x, y)
	{
		var c = false;
		for (var p = this.worldSpacePoints, i = -1, l = this.length, j = l - 1; ++i < l; j = i)
		{
			((p[i].y <= y && y < p[j].y) || (p[j].y <= y && y < p[i].y))
			&& (x <= (p[j].x - p[i].x) * (y - p[i].y) / (p[j].y - p[i].y) + p[i].x)
			&& (c = !c);
		}
		return c;
	}

	/** 
	* Generate contact points
	*/ 
	
	PHY2D.Polygon.generateContacts = function (B)
	{
		if (B.type === "poly")   return this.polyVsPoly(B);
		if (B.type === "circle") return this.polyVsCircle(B);
	}
	
	/** 
	* Polygon vs Polygon collision
	* http://www.wildbunny.co.uk/blog/2011/04/20/collision-detection-for-dummies/
	* Minkowski Difference
	*/ 

	PHY2D.Polygon.polyVsPoly = function (that)
	{
		var face, vertex, vertexRect, faceRect, fp, va, vb, vc, wsN, wdV0, wdV1, wsV0, wsV1;
		var mp = PHY2D.mostPenetrating, ms = PHY2D.mostSeparated;
		ms.set(100000, -1, -1, 0);
		mp.set(-100000, -1, -1, 0);
		// face of A, vertices of B
		this.featurePairJudgement(that, 2);
		// faces of B, vertices of A
		that.featurePairJudgement(this, 1);
		if (ms.dist > 0 && ms.fpc !== 0)
		{
			// objects are separated
			face = ms.edge;
			vertex = ms.closestI;
			fp = ms.fpc;
		}
		else if (mp.dist <= 0)
		{
			// objects are penetrating
			face = mp.edge;
			vertex = mp.closestI;
			fp = mp.fpc;
		}
		if (fp === 1)
		{
			vertexRect = this;
			faceRect = that;
		}
		else
		{
			vertexRect = that;
			faceRect = this;
		}
		wsN = faceRect.worldSpaceNormals[face];
		va = vertexRect.worldSpacePoints[(vertex - 1 + vertexRect.length) % vertexRect.length];
		vb = vertexRect.worldSpacePoints[vertex];
		vc = vertexRect.worldSpacePoints[(vertex + 1) % vertexRect.length];
		if (PHY2D.v0.project(vb, va, wsN) < PHY2D.v1.project(vc, vb, wsN))
		{
			wdV0 = va;
			wdV1 = vb;
		}
		else
		{
			wdV0 = vb;
			wdV1 = vc;
		}
		// world space edge
		wsV0 = faceRect.worldSpacePoints[face];
		wsV1 = faceRect.worldSpacePoints[(face + 1) % faceRect.length];
		// form contact
		if (fp === 1)
		{
			this.c0.projectPointOntoEdge(wsV0, wdV0, wdV1);
			this.c1.projectPointOntoEdge(wsV1, wdV0, wdV1);
			that.c0.projectPointOntoEdge(wdV1, wsV0, wsV1);
			that.c1.projectPointOntoEdge(wdV0, wsV0, wsV1);				
			this.normal.inv(wsN);
		}
		else
		{
			this.c0.projectPointOntoEdge(wdV1, wsV0, wsV1);
			this.c1.projectPointOntoEdge(wdV0, wsV0, wsV1);
			that.c0.projectPointOntoEdge(wsV0, wdV0, wdV1);
			that.c1.projectPointOntoEdge(wsV1, wdV0, wdV1);
			this.normal.copy(wsN);
		}
		var d0 = PHY2D.v1.sub(that.c0, this.c0).dot(this.normal);
		var d1 = PHY2D.v1.sub(that.c1, this.c1).dot(this.normal);
		// always return closest contact first
		if ( d0 < d1 )
		{
			PHY2D.contacts.create(this, that, this.c0, that.c0, d0, this.normal);
			PHY2D.contacts.create(this, that, this.c1, that.c1, d1, this.normal);
		}
		else
		{
			PHY2D.contacts.create(this, that, this.c1, that.c1, d1, this.normal);
			PHY2D.contacts.create(this, that, this.c0, that.c0, d0, this.normal);
		}
	}
	
	/** 
	* Polygon vs circle
	*/ 
	
	PHY2D.Polygon.polyVsCircle = function (that)
	{
		var mp = PHY2D.mostPenetrating, ms = PHY2D.mostSeparated;
		var v1 = PHY2D.v1;
		ms.set(100000, -1, -1, 0);
		mp.set(-100000, -1, -1, 0);
		var wsN = this.worldSpaceNormals;
		var wsP = this.worldSpacePoints;
		// track distance from circle centre to polygon
		for (var i = 0; i < this.length; i++)
		{
			// world space edge
			var wsV0 = wsP[i];
			var wsV1 = wsP[(i + 1) % this.length];
			// project onto poly edge and get distance
			var dist = v1.projectPointOntoEdge(that.pos, wsV0, wsV1).dist(that.pos);
			// track negative
			if (dist > mp.dist)
			{
				mp.dist = dist;
				mp.face = i;
				this.c0.copy(v1);
			}
			// track positive
			if (dist > 0 && dist < ms.dist)
			{
				ms.dist = dist;
				ms.face = i;
				this.c1.copy(v1);
			}
		}
		// account for circle radius - the above check was for sphere centre
		var penetration = mp.dist - that.radius;
		if (penetration <= 0)
		{
			// penetration
			if (mp.dist <= 0)
			{
				// normal is always a face normal
				this.normal.copy(this.worldSpaceNormals[mp.face]);
			}
			else
			{
				// normal can still be vertex normal
				this.normal.sub(that.pos, this.c0).unit();
			}
		}
		else
		{
			// separation
			this.normal.sub(that.pos, this.c1).unit();
			// normal can also be vertex normal
			penetration = ms.dist - that.radius;
			this.c0.copy(this.c1);
			
		}
		this.c1.subScale(that.pos, this.normal, that.radius);
		PHY2D.contacts.create(this, that, this.c0, this.c1, penetration, this.normal);
	}
	
	/** 
	* Feature Pair Judgement
	*/ 

	PHY2D.Polygon.featurePairJudgement = function (that, fpc)
	{
		var wsN, closestI, closest, dist;
		var mp = PHY2D.mostPenetrating, ms = PHY2D.mostSeparated;
		for (var edge = 0; edge < this.length; edge++)
		{
			// get support vertices
			wsN = this.worldSpaceNormals[edge];
			// rotate into RigidBody space
			PHY2D.v5.rotateIntoSpaceOf(wsN, that.matrix);
			var closestI = -1,
				closestD = -100000;
			// Get the vertex most in the direction of the given vector
			for (var i = 0; i < that.length; i++)
			{
				var d = PHY2D.v5.dot(that.localSpacePoints[i]);
				if (d > closestD)
				{
					closestD = d;
					closestI = i;
				}
			}
			var closest = that.worldSpacePoints[closestI];
			PHY2D.v0.sub(closest, this.worldSpacePoints[edge]);
			// distance from origin to face	
			var dist = PHY2D.v0.dot(wsN);
			if (dist > 0)
			{
				// recompute distance to clamped edge
				PHY2D.v1.sub(closest, this.worldSpacePoints[(edge + 1) % this.length]);
				// project onto minkowski edge
				dist = this.c0.projectPointOntoEdge(PHY2D.z0, PHY2D.v0, PHY2D.v1).lenSqr();
				// track separation
				if (dist < ms.dist)
				{
					ms.set(dist, closestI, edge, fpc);
				}
			}
			else
			{
				// track penetration
				if (dist > mp.dist)
				{
					mp.set(dist, closestI, edge, fpc);
				}
			}
		}
	}
	
	/** 
	* Circle PROTOTYPE definition
	* extends RigidBody base prototype
	*/ 
	
	PHY2D.Circle = Object.create(PHY2D.RigidBody);
	
	/** 
	* Init new circle
	*/ 
	
	PHY2D.Circle.init = function (x, y, w, mass, img)
	{
		this.type = "circle";
		this.initBase(x, y, w, w, mass, 0, img);
		this.radius = w / 2;
		this.invI = (1 / 10) * mass * this.radius * this.radius;
		return this;
	}
	
	/** 
	* Is point in circle ?
	*/ 
	
	PHY2D.Circle.isPointInPoly = function (x, y)
	{
		var dx = (x - this.pos.x);
		var dy = (y - this.pos.y);
		if (dx * dx + dy * dy < this.radius * this.radius) return true;
		return false;
	}
	
	/** 
	* AABB motion -- circle version
	*/ 
	
	PHY2D.Circle.motionAABB = function ()
	{
		var x = this.pos.x + (this.vel.x * PHY2D.kTimeStep);
		var y = this.pos.y + (this.vel.y * PHY2D.kTimeStep);
		this.aabb.x = (this.pos.x + x) * 0.5;
		this.aabb.y = (this.pos.y + y) * 0.5;
		this.aabb.w = this.radius + Math.abs(this.pos.x - x) * 0.5;
		this.aabb.h = this.radius + Math.abs(this.pos.y - y) * 0.5;
	}
	
	/** 
	* generate contact point (circle vs circle)
	*/ 
	
	PHY2D.Circle.generateContacts = function (that)
	{
		if (that.type === "circle")
		{
			this.normal.sub(that.pos, this.pos);
			var distCentre = this.normal.len();
			var penetration = distCentre - (this.radius + that.radius);
			if ( distCentre === 0 )
			{
				// default for circles exactly on top of each other
				this.normal.set(1, 0);
				penetration = -(this.radius + that.radius);
			}
			else
			{
				this.normal.unitLen(distCentre);
			}
			this.c0.addScale(this.pos, this.normal, this.radius);
			this.c1.subScale(that.pos, this.normal, that.radius);
			// generate point on edge
			PHY2D.contacts.create(this, that, this.c0, this.c1, penetration, this.normal);
		}
		else
		{
			// jump to poly vs circle
			that.polyVsCircle(this);
		}
	}
	
	/** 
	* Feature Pair Container
	*/ 

	PHY2D.FeaturePair = function ()
	{
		this.dist     = 0;
		this.closestI = 0;
		this.edge     = 0;
		this.fpc      = 0;
	}
	
	PHY2D.FeaturePair.prototype.set = function (dist, closestI, edge, fpc)
	{
		this.dist = dist;
		this.closestI = closestI;
		this.edge = edge;
		this.fpc = fpc;
	}
	
	PHY2D.mostSeparated   = new PHY2D.FeaturePair();
	PHY2D.mostPenetrating = new PHY2D.FeaturePair();

	/** 
	* Contacts Constructor
	*/ 

	PHY2D.Contact = function ()
	{
		this.a           = null;
		this.b           = null;
		this.normal      = new PHY2D.Vector();
		this.tangent     = new PHY2D.Vector();
		this.ra          = new PHY2D.Vector();
		this.rb          = new PHY2D.Vector();
		this.impulse     = new PHY2D.Vector();
		this.dv          = new PHY2D.Vector();
		this.v1          = new PHY2D.Vector();
		this.v2          = new PHY2D.Vector();
		this.dist        = 0;
		this.impulseN    = 0;
		this.impulseT    = 0;
		this.invDenom    = 0;
		this.invDenomTan = 0;
	}
	
	
	PHY2D.Contact.prototype = {
		
		/** 
		* set existing contact object
		*/ 
		
		set: function (A, B, pa, pb, dist, normal)
		{
			var ran, rbn;
			this.a = A;
			this.b = B;
			this.normal.copy(normal);
			this.tangent.perp(normal);
			this.dist = dist;
			this.impulseN = 0;
			this.impulseT = 0;
			// calculate radius arms
			this.ra.sub(pa, A.pos).perpSelf();
			this.rb.sub(pb, B.pos).perpSelf();
			// compute denominator in impulse equation
			ran = this.ra.dot(this.normal);
			rbn = this.rb.dot(this.normal);
			this.invDenom = 1 / (A.invMass + B.invMass + (ran * ran * A.invI) + (rbn * rbn * B.invI));
			ran = this.ra.dot(this.tangent);
			rbn = this.rb.dot(this.tangent);
			this.invDenomTan = 1 / (A.invMass + B.invMass + (ran * ran * A.invI) + (rbn * rbn * B.invI));
		},

		/** 
		* Apply impulse
		*/ 
		
		applyImpulse: function ()
		{
			// linear
			this.a.vel.addScale(this.a.vel, this.impulse, this.a.invMass);
			this.b.vel.subScale(this.b.vel, this.impulse, this.b.invMass);
			// angular
			this.a.angVel += this.impulse.dot(this.ra) * this.a.invI;
			this.b.angVel -= this.impulse.dot(this.rb) * this.b.invI;
		},

		/** 
		* Solve contact
		*/ 

		solve: function ()
		{
			var newImpulse, absMag;
			// get all of relative normal velocity
			this.dv.sub(
				this.v1.addScale(this.b.vel, this.rb, this.b.angVel),
				this.v2.addScale(this.a.vel, this.ra, this.a.angVel)
			);
			// accumulated impulse
			newImpulse = (this.dv.dot(this.normal) + this.dist / PHY2D.kTimeStep) * this.invDenom + this.impulseN;
			// push only
			if (newImpulse > 0) newImpulse = 0;
			// apply impulse
			this.impulse.scale(this.normal, newImpulse - this.impulseN)
			this.applyImpulse();
			this.impulseN = newImpulse;
			// friction
			if (this.dist <= 0)
			{
				absMag = -this.impulseN * PHY2D.kFriction;
				newImpulse = this.impulse.clamp(this.dv.dot(this.tangent) * this.invDenomTan + this.impulseT, -absMag, absMag);
				this.impulse.scale(this.tangent, newImpulse - this.impulseT);
				this.applyImpulse();
				this.impulseT = newImpulse;
			}
		}
	}
	
	/** 
	* constraint (joint) constructor
	*/ 
	
	PHY2D.Constraint = function (A, va, B, vb)
	{
		this.a        = A;
		this.b        = B;
		this.ra       = new PHY2D.Vector();
		this.rb       = new PHY2D.Vector();
		this.axis     = new PHY2D.Vector();
		this.normal   = new PHY2D.Vector();
		this.impulse  = new PHY2D.Vector();
		this.dv       = new PHY2D.Vector();
		this.va       = new PHY2D.Vector(va[0] * PHY2D.zoom, va[1] * PHY2D.zoom);
		this.vb       = new PHY2D.Vector(vb[0] * PHY2D.zoom, vb[1] * PHY2D.zoom);
		this.pa       = new PHY2D.Vector();
		this.pb       = new PHY2D.Vector();
		this.invDenom = 0;
		this.impulseN = 0;
		return this;
	}
	
	/** 
	* constraint pre-solve
	*/ 
	
	PHY2D.Constraint.prototype.presolve = function ()
	{
		// transform constraint points
		this.pa.transform(this.va, this.a.matrix);
		this.pb.transform(this.vb, this.b.matrix);
		// projection axis
		this.axis.sub(this.pb, this.pa);
		this.normal.scale(this.axis, 1 / (this.axis.len() + 0.0001));
		// calculate radius arms
		this.ra.sub(this.pa, this.a.pos).perpSelf();
		this.rb.sub(this.pb, this.b.pos).perpSelf();
		// compute denominator in impulse equation
		var ran = this.ra.dot(this.normal), rbn = this.rb.dot(this.normal);
		this.invDenom = 1 / (this.a.invMass + this.b.invMass + (ran * ran * this.a.invI) + (rbn * rbn * this.b.invI));
	}
	
	/** 
	* constraint solve
	*/ 
	
	PHY2D.Constraint.prototype.solve = function ()
	{
		// calculate relative velocity in the axis, we want to remove this
		this.dv.sub(
			PHY2D.v1.addScale(this.b.vel, this.rb, this.b.angVel),
			PHY2D.v2.addScale(this.a.vel, this.ra, this.a.angVel)
		);
		// accumulated impulse
		var dist = this.axis.dot(this.normal);
		var newImpulse = (this.dv.dot(this.normal) + dist / PHY2D.kTimeStep) * this.invDenom + this.impulseN;
		// apply impulse
		this.impulse.scale(this.normal, newImpulse - this.impulseN)
		this.applyImpulse();
		this.impulseN = newImpulse;
	}
	
	/** 
	* Apply Impulse
	*/ 
	
	PHY2D.Constraint.prototype.applyImpulse = PHY2D.Contact.prototype.applyImpulse;
	
	/** 
	* Manipulate Constructor
	*/ 
		
	PHY2D.Manipulate = function (B)
	{
		this.b        = B;
		this.rb       = new PHY2D.Vector();
		this.axis     = new PHY2D.Vector();
		this.normal   = new PHY2D.Vector();
		this.impulse  = new PHY2D.Vector();
		this.dv       = new PHY2D.Vector();
		this.vb = (B.type === "circle") ? new PHY2D.Vector() : new PHY2D.Vector().rotateIntoSpaceOf(PHY2D.v1.sub(B.pos, pointer.pos), B.matrix);
		this.pb       = new PHY2D.Vector();
		return this;
	}
	
	/** 
	* solve Manipulate Constraint
	*/ 
	
	PHY2D.Manipulate.prototype.solve = function ()
	{
		this.pb.transform(this.vb, this.b.matrix);
		this.axis.sub(this.pb, pointer.pos);
		this.normal.scale(this.axis, 1 / (this.axis.len() + 0.0001));
		this.rb.sub(this.pb, this.b.pos).perpSelf();
		this.dv.addScale(this.b.vel, this.rb, this.b.angVel);
		this.impulse.scale(this.normal, this.dv.dot(this.normal) + this.axis.dot(this.normal) * 2);
		this.b.vel.sub(this.b.vel, this.impulse);
		this.b.angVel -= this.impulse.dot(this.rb) * this.b.invI;
	}
	
	/** 
	* Rendering
	*/ 
	
	PHY2D.render = function ()
	{

		/** 
		* brute force AABB broad phase
		* SAP, grids, quadtrees don't do any better < 100 objects
		*/ 
	
		PHY2D.contacts.index = 0;
		for (var i = 0, len = PHY2D.objects.length; i < len - 1; i++)
		{
			var A = PHY2D.objects[i];
			for (var j = i + 1; j < len; j++)
			{
				var B = PHY2D.objects[j];
				if (A.invMass || B.invMass)
				{
					var a = A.aabb,
						b = B.aabb;
					if (
						Math.abs(b.x - a.x) - (a.w + b.w) < 0 &&
						Math.abs(b.y - a.y) - (a.h + b.h) < 0
					) {
					
						/** 
						* generate speculative contact
						*/ 
						
						A.generateContacts(B);
						
					}
				}
			}
		}
		
		var contacts = PHY2D.contacts, len = contacts.index;
		var joints = PHY2D.joints, cln = joints.length;

		/** 
		* pre-solve constraints
		*/ 
		
		for (var i = 0; i < cln; i++)
		{
			PHY2D.joints[i].presolve();
		}
		
		
		for (var j = 0; j < 5; j++) // numIterations
		{
			
			/** 
			* solve constraints
			*/ 
			
			for (var i = 0; i < cln; i++)
			{
				PHY2D.joints[i].solve();
			}
			
			/** 
			* solve contacts
			*/ 
			
			for (var i = 0; i < len; i++)
			{
				PHY2D.contacts[i].solve();
			}
		}
		
		/** 
		* integration
		*/ 

		for (var i = 0, len = PHY2D.objects.length; i < len; i++)
		{
			PHY2D.objects[i].integrate();
		}
		
		/** 
		* draw images
		*/ 

		for (var i = 0; i < len; i++)
		{
			var rb = PHY2D.objects[i];
			rb.draw();
		}
		
		/** 
		* Manipulate object
		*/ 
		
		if (PHY2D.manipulate)
		{
			PHY2D.manipulate.solve();
			ctx.beginPath();
			ctx.lineWidth = 2;
			ctx.setLineDash([2,2]);
			ctx.moveTo(pointer.pos.x, pointer.pos.y);
			ctx.lineTo(PHY2D.manipulate.pb.x, PHY2D.manipulate.pb.y);
			ctx.strokeStyle = "rgb(255,255,255)";
			ctx.stroke();
			ctx.closePath();
			ctx.fillStyle = "rgb(255,255,255)";
			ctx.arc(PHY2D.manipulate.pb.x, PHY2D.manipulate.pb.y, 5, 0, 2 * Math.PI);
			ctx.fill();
		}
	}
	
	/** 
	* create new rectangle
	*/ 
	
	PHY2D.rectangle = function (x, y, w, h, mass, angle, img)
	{
		x *= this.zoom;
		y *= this.zoom;
		w *= this.zoom;
		h *= this.zoom;
		var vertices = [
			w / 2, -h / 2, -w / 2, -h / 2, -w / 2, h / 2, w / 2, h / 2
		];
		var rb = Object.create(PHY2D.Polygon).init(x, y, w, h, vertices, mass, angle, img);
		PHY2D.objects.push(rb);
		return rb;
	}
	
	/** 
	* create new triangle
	*/ 
	
	PHY2D.triangle = function (x, y, w, h, mass, angle, img)
	{
		x *= this.zoom;
		y *= this.zoom;
		w *= this.zoom;
		h *= this.zoom;
		var vertices = [
			w / 2, h / 2, 0, - h / 2, -w / 2, h / 2
		];
		var rb = Object.create(PHY2D.Polygon).init(x, y, w, h, vertices, mass, angle, img);
		PHY2D.objects.push(rb);
		return rb;
	}
	
	/** 
	* create new circle
	*/ 
	
	PHY2D.circle = function (x, y, w, mass, img)
	{
		x *= this.zoom;
		y *= this.zoom;
		w *= this.zoom;
		var rb = Object.create(PHY2D.Circle).init(x, y, w, mass, img);
		PHY2D.objects.push(rb);
		return rb;
	}
	
	/** 
	* create new constraint
	*/ 
	
	PHY2D.addJoint = function (A, va, B, vb)
	{
		PHY2D.joints.push(
			new PHY2D.Constraint(A, va, B, vb)
		);
	}
	
	/** 
	* delete static objects
	* used for screen resize
	*/ 
	
	PHY2D.deleteStatic = function ()
	{
		var k = PHY2D.objects.length;
		while (k--)
		{
			var p = PHY2D.objects[k];
			if (!p.invMass)
			{
				PHY2D.objects.splice(k, 1);
			}
		}
	}
	
	/** 
	* default engine zoom
	*/ 
	
	PHY2D.zoom = 1;
	
	/** 
	* cancel manipulation
	*/ 
	
	PHY2D.stopManipulate = function () { PHY2D.manipulate = null; }
	

	////////////////////////////////////////////////////////////////

	screen.resize();
	
	for (var i = 0; i < 6; i++) PHY2D.circle(100 + i * 200, screen.height - 100, 200, 1, {src: "e6.png"});
	for (var i = 1; i < 6; i++) PHY2D.circle(100 + i * 233, screen.height - 333, 233, 1, {src: "e8.png"});
	PHY2D.rectangle(500, -400, 370, 420, 4, Math.PI / 2, {
		src: "p25.png", 
		x: 10, 
		w: 380
	});
	
	// ==== main loop ====
	
	function run()
	{
		requestAnimationFrame(run);
		ctx.clearRect(0, 0, screen.width, screen.height);
		PHY2D.render();
	}
	
	// ==== start animation ====
	run();
	
}();