//Original from tutorial
public Vector Wheel::CalculateForce(Vector relativeGroundSpeed, float dt)
{
    //calculate speed of tire patch at ground
    Vector patchSpeed = LocalForwardDir * PatchSpeed;

    //get velocity difference between ground and patch
    Vector velDifference = relativeGroundSpeed + patchSpeed;

    //project ground speed onto side axis
    float forwardMag = 0;
    Vector sideVel = velDifference.Project(mSideAxis);
    Vector forwardVel = velDifference.Project(mForwardAxis, out forwardMag);

    //calculate super fake friction forces
    //calculate response force
    Vector responseForce = -sideVel * 2.0f;
    responseForce -= forwardVel;

    //calculate torque on wheel
    mTorque += forwardMag * mRadius;

    //return force acting on body
    return responseForce;
}


// New physics model:
public void ApplyWheelConstraint(Wheel wheel, Vector worldGravity, float dt)
{
    // Gather a lot of data
    float maxKineticAcc = cGravity * wheel.mKineticFrictionCoef;
    float maxKineticDV = maxKineticAcc / mMass * dt;

    Vector offset = ToWorld(wheel.mPosition);
    Vector forwardDir = ToWorld(wheel.LocalForwardDir);
    Vector constraintDir = ToWorld(wheel.LocalConstraintDir);

    Vector groundVel = PointVel(offset);
    Vector constraintVel = groundVel.Project(constraintDir);
    Vector forwardVel = groundVel.Project(forwardDir);
    float  forwardVelMag = forwardVel.Length;
    Vector patchVel = forwardDir * wheel.PatchSpeed;
    Vector patchVelDiff = groundVel + patchVel;
    float  patchVelDiffMag = patchVelDiff.Length;

    Vector torqueDV = new Vector( );

    if (wheel.mStaticFriction)
        torqueDV = forwardDir * wheel.ComputeTorqueDV(mMass, dt);
    else
    {
        float forwardVelDiff = patchVelDiff.Dot(forwardDir);
        torqueDV = forwardDir * forwardVelDiff;
    }

    // Accumulate accelerations
    Vector desiredDeltaV = constraintVel + torqueDV;

    // Compute the corrective impulse and its resultant acceleration
    Vector deltaVDir = desiredDeltaV;
    float deltaVMag = 0.0f;
    deltaVDir.Normalize(out deltaVMag);
    Vector roadImpulse = ComputeImpulseToChangePointVel(dt, offset, deltaVDir, deltaVMag);
    Vector roadImpulseAcc = roadImpulse / mMass / dt;

    // Determine whether to switch to kinetic or static friction
    Vector frictionAcc = roadImpulseAcc;
    float accMag = frictionAcc.Length;

    if (patchVelDiffMag < 0.25f)
        wheel.mStaticFriction = true;
    else if (wheel.mStaticFriction)
    {
        // Static friction
        float maxAcc = cGravity * wheel.Friction;
        float totalPatchAcc = accMag;
        if (totalPatchAcc > maxAcc)
            wheel.mStaticFriction = false;
    }

    if (!wheel.mStaticFriction)
    {
        // Kinetic friction, clamp acceleration
        accMag = Math.Min(accMag, maxKineticAcc);
        frictionAcc = frictionAcc.Normalize() * accMag;
        roadImpulse = frictionAcc * mMass * dt;
    }

    // Gather some debugging stats
    wheel.mStats.mWorldLinearAcceleration.mValue = frictionAcc;

    // Apply result
    Vector totalImpulse = roadImpulse;
    AddImpulse(totalImpulse, offset);
    Vector newGroundV = PointVel(offset);

    // Integrate wheel
    if (wheel.mStaticFriction)
    {
        wheel.mAngVel = newGroundV.Dot(forwardDir) / wheel.mRadius;
    }
    else
    {
        float frictionForce = frictionAcc.Dot(forwardDir) * mMass;
        float frictionTorque = frictionForce * wheel.mRadius;
        wheel.mAngVel += frictionTorque / wheel.mInertia * dt;
    }

    wheel.mAngle += wheel.mAngVel * dt;
}