Roughly how physics constraints work

    /*
     * The idea of a constraint is entirely based around a cost function. A cost function takes in both the
     * position (linear position) and orientation (angular position) of both entityA and entityB, and outputs
     * a scalar, which is 0 if the configuration of the entities is valid (in the feasible region), and nonzero
     * otherwise. We denote the cost function as C(xA, rA, xB, rB). For this comment, pA and pB will denote
     * the anchor points (points where the joint is applied) of each entity in world coordinates. We will also
     * lay out the components of each entity's position and orientation in a flattened 12D vector, the same
     * for velocity.
     *
     * The cost function for 3D constraints is hard to visualize, but is sometimes easy in 2D. Imagine a fixed
     * distance constraint between two points A and B on the xy-plane. Let C(xA, rA, xB, rB) = ||pB - pA|| - d
     * where d is the target distance. If you plot pA and pB and the output of the cost function on a third
     * axis, you get an upside-down cone shape, and the cost function is 0 in the feasible region. Take the
     * time to visualize this. It should also be clear that the gradient of C points perpendicular to the
     * feasible region, so in order to get back into the feasible region you must apply an impulse in the
     * direction of this gradient, and how big an impulse (may be negative) depends on the cost. This is true
     * for the cost functions of all constraints, in 3D too although you would need at least 4D to visualize.
     *
     * But we need to correct the velocities of the objects as well as the position, so the first step is
     * differentiating the cost function with respect to time to get a velocity cost. By the chain rule:
     * dC/dt = dC/dx * dx/dt. dx/dt is just velocity. dC/dx gets you a 12D row vector of all the partial
     * derivatives of the cost function with respect to each component of x. This vector is called the
     * jacobian. We will also change dC/dt to cDot because physics. So we get cDot = J * v. The jacobian,
     * being the derivative of C with respect to position, is also the direction we must apply the impulse.
     *
     * The main difficulty with making a new type of constraint is that you have to hand-derive this jacobian.
     * This step is specific to each constraint and each subclass has a comment like this one with the
     * derivation. To derive the jacobian, you don't directly compute dC/dx, instead you directly compute
     * dC/dt, then simply rearrange into the form cDot = J * v ("isolate the velocities"), then just read the
     * jacobian by inspection.
     *
     * This is all you need to know to implement a constraint, so engineers stop reading. However we haven't
     * done the maths yet to complete how the constraint actually works.
     *
     * Since the jacobian is the direction we apply the impulse, the force, f (a vector including linear forces
     * and angular torques), is parallel to this vector, therefore a multiple of it:
     * f = J(T) * lambda
     * where (T) is the matrix transpose operation and lambda is a scalar.
     *
     * We also know that for a discrete physics time step:
     * vFinal = vInitial + a * dt
     * where a is an acceleration vector.
     *
     * By Newton's second law, for a constant mass matrix M:
     * totalF = M * a
     * a = M^-1 * totalF
     * a = M^-1 * (fExt + fC)
     * where totalF is the total force applied, fExt is the external forces and fC is the forces applied due to
     * constraints. What is going to happen to these forces is that they are going to be added to momentum every
     * time step, essentially numerical integration. Since integration is additive, we can integrate fExt and fC
     * separately, and since we only care about fC, we get:
     * aC = M^-1 * fC
     * and since we're lazy we're going to be naughty and redefine a:
     * a = M^-1 * fC
     *
     * The mass matrix is a 12x12 matrix which is mostly zeros. You get 4 3x3 matrices along the leading diagonal,
     * containing mass (linear mass) and inertia (angular mass). Linear mass is just scalar so we insert a scalar
     * multiple of the identity matrix. Inertia may not be so simple, it may be a tensor.
     *
     * Substituting, we get
     * vFinal = vInitial + M^-1 * fC * dt
     * Substituting again:
     * vFinal = vInitial + M^-1 * J(T) * lambda * dt
     * And again:
     * cDotFinal = J * (vInitial + M^-1 * J(T) * lambda * dt)
     * Rearrange:
     * cDotFinal = J * vInitial   +   J * M^-1 * J(T) * lambda * dt
     * J * M^-1 * J(T) * lambda * dt = cDotFinal - J * vInitial
     *
     * We also know that:
     * f * dt = p
     * where p is the impulse vector.
     * So to describe impulse rather than force, we're going to be slightly naughty again and redefine lambda
     * to describe that multiple for impulse rather than for force:
     * J * M^-1 * J(T) * lambda = cDotFinal - J * vInitial
     *
     * Since we aim for C = 0, we also aim for cDotFinal = 0, so:
     * J * M^-1 * J(T) * lambda = -J * vInitial
     * Rearrange:
     * lambda = (J * M^-1 * J(T))^-1 * -J * vInitial
     *
     * This J * M^-1 * J(T) term is so important for the calculations that we give it its own letter, K:
     * lambda = K^-1 * -J * vInitial
     * K^-1 is also called the "constraint mass"; K and therefore K^-1 are also scalars.
     *
     * We have lambda! We're almost there.
     * J * vInitial = cDotInitial, so:
     * lambda = K^-1 * -cDotInitial
     * lambda = - K^-1 * cDotInitial
     *
     * And remember the definition of lambda: the scalar multiple of the jacobian, which was the direction of the
     * impulse, to get the needed impulse to put cDotFinal to 0. So:
     *
     * impulse = J(T) * -(K^-1 * cDotInitial)
     *
     * And this is the formula you find in solveVelocityConstraints().
     */