For planar robots, we can calculate how small changes in the angle of each joint corresponds to the displacement of the end effector, by constructing a jacobian consisting of cross products between the axis of rotation of each joint and difference between the position of each joint and the end joint. I tried to use the same logic for a spherical joint with 3dof by calculating three angles for each global axis of rotation, from the displacement of the end effector at each global axis. I then tried to construct a rotational matrix to correctly orient the joint from these three Euler angles. However, since rotations of each angle stack on top of each other in rotation matrices, this didn't work. Is there a practical way to construct a rotation matrix from these angles or is there a different solution to this problem?
Write the rotation matrix from the first coordinate system to the coordinate system moved by the first motor. Four of the elements should include a cosine or sine of the first motor angle. Then write the rotation matrix from that coordinate system to the coordinate system moved by the second motor. Do this again for the third motor.
Multiply these three rotation matrixes together and you will get the rotation matrix for the 3 dof spherical joint.
Incidentally, you can easily write the individual 3x3 rotation matrices as follows:
The first column is composed of the direction cosines of the moving x axis relative to the fixed coordinate system; the second column is composed of the direction cosines of the moving y axis with respect to the fixed coordinate system; the third column is composed of the direction cosines of the moving z axis with respect to the fixed coordinate system.