I'm a bit more ashamed than you not knowing how a physics solver actually works :)

Newton engine builds the jacobians, what i give as input for a motor is angular acceleration about an axis.

In theory we would need only need one axis and angle to get any rotation we want, but newton likes to do this for 3 orthonormal axis, so everything is defined and stable. I guess something similar applies to you.

I refer to the orientation build from this 3 axis as 'constraint space' (not sure it's the proper term), and the primary axis where the moter acts as 'error axis' (axis from Dirks 'qr' quat).

The problem is: Newton wants constraints order independent, which makes sense because we want them to be solved at the same time.

But it is impossible to decompose a rotation to multiple rotations so you can combine (multiply) them together in any order and get the initial result. Whatever order you take, you get a different result.

So my trick here limits the problem, because at least the measurement of angles is guaranteed to be order independent, and by aligning constraint space with error axis any error becomes negligible.

However, if your (or Dirks) solver solves the constraints in given order, then you don't have this problem.

The problem in your video looks like a jump over 180 degrees or similar gimbal lock issue coming from eulers. I guess not your jacobians are wrong but more likely the angle measurement.

But i also can see some unnatural twisting, also a typical euler issue. Euler angles == evil. Both Dirks or mine suggestion should fix those things.

Other than that it looks good - seems you're not far off.

I tried to work on twist / swing again, but Newton seems somehow broken at the moment due to heavy changes for new softbody feature. I need to clarify this...

@JoeJ

I changed my angle computation to figure out the error using your suggested method. It seemed to make it better -- the method in which it fails now appears to be consistent -- it will correctly rotate unless the target angle for any degree of freedom crosses the 180 degree mark. When moving the target back inside the working quadrant, it becomes stable again. Aside from that, it appears to be stable, even when it's shot with the gun.

I remember that when you do an Euler decomposition the middle angle can vary a lot. So you need to be careful with this. Personally I have never used Euler angles for the above reasons.

I realized that I might be using axis angle, not euler angles for angular velocity. I didn't realize they're not the same thing.

This is what I am doing to compute the velocity bias:

                const XMVECTOR targetAngles = XMVectorSet(xRot, yRot, 0, 0);
		const XMVECTOR targetQuaternion = XMQuaternionRotationRollPitchYawFromVector(targetAngles);
		const XMVECTOR quaternionBetweenObjects = XMQuaternionMultiply(this->object2->qOrientation, XMQuaternionInverse(this->object1->qOrientation));
		const XMVECTOR quaternionError = XMQuaternionMultiply(quaternionBetweenObjects, XMQuaternionInverse(targetQuaternion));
		
		XMVECTOR axisBetweenObjects;
		float angleBetweenObjects;
		XMQuaternionToAxisAngle(&axisBetweenObjects, &angleBetweenObjects, quaternionError);
		
		const XMVECTOR targetAngularVelocity = XMVectorMultiply(XMVector3Normalize(axisBetweenObjects), XMVectorReplicate(angleBetweenObjects));
		this->velocityBias = XMVectorMultiply(targetAngularVelocity, XMVectorReplicate(-FractalTreeResources::angleBB));

Whenever yRot in targetAngles is between the radian equivalent of -90 degrees and +90 degrees, it works correctly with any value for xRot. If xRot is 0, any value for yRot works. If xRot is non-zero, then as soon as yRot is no longer set to be in the range of -90 and +90, it explodes. Moving yRot back into the working range un-explodes it and it becomes stable again, even after being unstable for a long time. After this works, the values of targetAngles will come from my inverse kinematics solver, so it needs 3 degrees of freedom because the robot arm segments can have the ability to move on any axis the player sets.

Dirk at one point helped me go through derivation of some joints. Maybe you could open a new topic to figure out the calculus for deriving the joint (or just reuse this topic). It's not too complicated and could be a great learning experience. I would be happy to try and help with the derivation. I understand the mathematics can seem very farfetched, but the problem here is truly a math problem.

Edit: Oh and I couldn't help with quaternion based derivations, but I could help with dot product based ones. Dot product does not give a full range of motion, but can be good enough to work with. To get started, you need to remove 2 DOF (two axes) and leave the third axis unrestricted. This adds two rows to jacobian. Then we add in limit on last axis for another row (when it is active), and finally the motor as a non-holonomic constraint (constraints velocity, no position level equation). When all constraints are active (including the 3 linear constraints) in total the joint would have 3 + 2 + 2 constraints. Once a good design is found to remove these freedoms, the rest is calculus.