If you run the code below, you can confirm the following with respect to the functions in the Aerospace Toolbox:
The rotation sequence coming out of the quat2angle function uses individual rotation matrices that are "coordinate system rotations of the x-, y-, and z-axes in a counterclockwise direction when looking towards the origin" using the terminology of this link:
The first output is applied first, then the second output is applied second, then the third output is applied last. So the order in which the outputs appears is the order in which they operate on a vector. When multiplied together, they get the same result as two other rotation methods seen in the code below (using quatrotate and quatmultiply).
The "sense" of the quaternion is that successive rotations are multiplied on the right (confirmed with the vq calculation). (If it had been the reverse, the calculation would have been q * vi * conj(q))
I am not sure what you mean by "intrinsic body fixed" vs "extrinsic space fixed", but maybe you can decipher that meaning from the output of the code below.
disp('------------------------------')
% sample data
q = [1 -2 -3 4]; q = q / norm(q)
vi = [1 2 3]';
% quat2angle
disp('------------------------------')
disp('dcm = quat2dcm(q)')
dcm = quat2dcm(q)
[a3,a2,a1] = quat2angle(q);
c = cos(a3);
s = sin(a3);
r3 = [ c s 0;
-s c 0;
0 0 1];
c = cos(a2);
s = sin(a2);
r2 = [ c 0 -s;
0 1 0;
s 0 c];
c = cos(a1);
s = sin(a1);
r1 = [ 1 0 0;
0 c s;
0 -s c];
disp('dce = result of individual rotation matrix multiplies')
dce = r1 * r2 * r3
% compare rotation methods
vb = dcm * vi;
vq = quatmultiply(quatconj(q),[0 vi']);
vq = quatmultiply(vq,q);
vr = quatrotate(q,vi');
disp('------------------------------')
disp('vi = input inertial vector')
disp(vi')
disp('vb = dcm * vi')
disp(vb')
disp('vq = conj(q) * vi * q')
disp(vq(2:4))
disp('vr = quatrotate(q,vi'')')
disp(vr)
disp('------------------------------')
