Function to Convert between DCM, Euler angles, Quaternions, and Euler vectors

Function which will convert any rotation data of 4 types: DCM, Euler Angles, Quaternions, and Euler

Description

SpinCalc is a consolidated matlab function that will convert any rotation data between the 4 types included. Will also convert between 2 different Euler angle set types.

Multiple orientations can be input. For N orientations:
DCM ===> 3x3xN multidimensional array
EA### ===> Nx3 matrix
Euler Vectors ===> Nx4 matrix
Quaternions ===> Nx4 matrix

Input includes error check flag that will warn when Euler angles approach singularity or when appropriate values deviate from unity. Fatal errors issued for improper DCM's etc.

*****NOTE TO USERS*****
I have gotten many questions regarding translation to Euler angle sets. When converting data to Euler angles, you MUST make sure the orientation you are translating is not near a singularity. Singular Euler sets are orientations which cannot be uniquely converted to 3 variables in that particular rotation order. The singular sets are as follows:

Type 1 Rotations: 123 - 132 - 213 - 231 - 321 - 312
Singular if second rotation angle is -90 or 90 degrees.

Type 2 Rotations: 121 - 131 - 212 - 232 - 313 - 323
Singular if second rotation angle is 0 or 180 degrees.

SpinCalc should now detect when input DCM, EV, or Q correspond to a singular Euler set output. It will prohibit output in such an event.

Naturally when converting from these singular Euler angle sets to other data types, you will receive a correct answer. Unfortunately you cannot convert that output back to the correct Euler angle set. This is why singular Euler input is prohibited.

• MATLAB release: MATLAB 7.5 (R2007b)
• Other Files: license.txt, ; SpinCalc.m

[James Tursa]

I spot tested all combinations of input/output type against independent code I happened to have and SpinCalc produced the same results, so I would conclude based on this spot checking that the code generally works as advertised. Thus I gave it a fairly high rating. I did have a few problems, but I think these can easily be improved:

------------------------------------------------------------------------------------
1) Sometimes the results were not invertible. For example:

>> q = [1 2 3 4]
q =
     1 2 3 4
>> q = q / norm(q)
q =
   0.18257418583506 0.36514837167011 0.54772255750517 0.73029674334022
>> SpinCalc('QtoEA213',q,eps,1)
ans =
  1.0e+002 *
   0.47726310993906 3.52337744339234 0.70346175941947
>> SpinCalc('EA213toQ',ans,eps,1)
??? Error using ==> SpinCalc
Error: Second input Euler angle(s) outside -90 to 90 degree range

In this case, and in other similar cases, the checks for valid Euler Angle input are too tight. Having the second input outside -90 to +90 may not make sense for some physical applications, but it is not an invalid input and should not produce an error.

------------------------------------------------------------------------------------
2) Calculating the euler rotation angle for small angle rotation quaternions is not accurate. The code in question is this line:

MU = 2*acos(Q(:,4));

For example:

>> q = [1 2 3 4e12]
q =
  1.0e+012 *
   0.00000000000100 0.00000000000200 0.00000000000300 4.00000000000000
>> q = q / norm(q)
q =
   0.00000000000025 0.00000000000050 0.00000000000075 1.00000000000000
>> SpinCalc('QtoEV',q,eps,1)
ans =
     1 0 0 0
>> SpinCalc('EVtoQ',ans,eps,1)
ans =
     0 0 0 1

You can see in this case that the small angle was completely lost in the q(4) component, so using only that component to calculate the rotation angle is not a good method for small angle cases. One fix is to use the atan2 function instead. For example, suppose the following code was used instead of the code above:

QV = sqrt(sum(Q(:,1:3).*Q(:,1:3),2));
MU = 2*atan2(QV,Q(:,4));

Now do the same calculations as above:

>> SpinCalc('QtoEV',q,eps,1)
ans =
   0.26726124191242 0.53452248382485 0.80178372573727 0.00000000010719
>> SpinCalc('EVtoQ',ans,eps,1)
ans =
   0.00000000000025 0.00000000000050 0.00000000000075 1.00000000000000

You can see that the small angle is calculated properly. Caveat: The atan2 code above is just shown as an example for improvement. It may not be the best overall way to do this.

------------------------------------------------------------------------------------
3) Some of the algorithms used are not robust for inputs that deviate from unity, etc. For example, suppose we start with a proper direction cosine matrix and perturb it a bit:

>> q = [1 2 3 4]
q =
     1 2 3 4
>> q = q / norm(q)
q =
   0.18257418583506 0.36514837167011 0.54772255750517 0.73029674334022
>> dcm=SpinCalc('QtoDCM',q,eps,1)
dcm =
   0.13333333333333 0.93333333333333 -0.33333333333333
  -0.66666666666667 0.33333333333333 0.66666666666667
   0.73333333333333 0.13333333333333 0.66666666666667
>> dcm = dcm + 1e-7*(rand(3,3)-0.5)
dcm =
   0.13333337834626 0.93333333193158 -0.33333333768657
  -0.66666669355282 0.33333337246323 0.66666661851703
   0.73333334401759 0.13333335954302 0.66666669880738
>> SpinCalc('DCMtoQ',dcm,eps,1)
Warning: Input DCM(s) matrix not orthogonal to precision tolerance.
(etc.)
ans =
   0.18257415540401 0.36514836686611 0.54772255130155 0.73029676324370
>> norm(ans)
ans =
   1.00000000382748

The algorithm used for the conversion picks up on the fact that the dcm was not quite proper (good), but does not correct the output. So maybe normalizing the returned quaternion might be better in this case. One way to accomplish this efficiently would be to use Markley's algorithm that was recently published in J. Guidance Vol 31 No 2 pp 440-442. This algorithm produces a normalized quaternion result using only one square root.

------------------------------------------------------------------------------------
4) Some of the code could be written more efficiently. For example in this section:

m1=INPUT(:,1); m2=INPUT(:,2); m3=INPUT(:,3); MU=INPUT(:,4)*pi/180;
if sqrt(m1.^2+m2.^2+m3.^2)-ones(N,1)>tol*ones(N,1),
    error('Input euler vector(s) components do not constitute a unit vector')
end
if MU<zeros(N,1) | MU>2*pi*ones(N,1),
    error('Input euler rotation angle(s) not between 0 and 360 degrees')
end
Q=[m1.*sin(MU/2),m2.*sin(MU/2),m3.*sin(MU/2),cos(MU/2)];

Why pull out the columns individually in m1, m2, and m3? Why not just use INPUT(:,1:3) in the follow-on calculations? And why calculate sin(MU/2) multiple times? Could recode this to calculate sin(MU/2) only once. For instance:

Q=[INPUT(:,1:3).*sin(MU/2),cos(MU/2)];

------------------------------------------------------------------------------------
5) Multiple warning messages appear for systematic problems when one message would suffice. For example:

>> q = [1 2 3 4]
q =
     1 2 3 4
>> q = q / norm(q)
q =
   0.18257418583506 0.36514837167011 0.54772255750517 0.73029674334022
>> q = q + 1.e-7*(rand(1,4)-.5)
q =
   0.18257418030539 0.36514838321335 0.54772258669887 0.73029678552152
>> qq = [q;q;q;q;q;q;q;q;q;q]
qq =
   0.18257418030539 0.36514838321335 0.54772258669887 0.73029678552152
   0.18257418030539 0.36514838321335 0.54772258669887 0.73029678552152
(etc.)
>> SpinCalc('QtoDCM',qq,eps,1)
ans(:,:,1) =
   0.13333335251383 0.93333342235746 -0.33333337639558
  -0.66666675533740 0.33333337341210 0.66666670795756
   0.73333338560076 0.13333335997254 0.66666675384566
ans(:,:,2) =
   0.13333335251383 0.93333342235746 -0.33333337639558
  -0.66666675533740 0.33333337341210 0.66666670795756
   0.73333338560076 0.13333335997254 0.66666675384566
(etc.)
>> SpinCalc('DCMtoQ',ans,eps,1)
Warning: Input DCM(s) matrix not orthogonal to precision tolerance.Warning: Input DCM(s) matrix not orthogonal to precision tolerance.
(etc.)
ans =
   0.18257418458450 0.36514839177157 0.54772259953620 0.73029676840507
   0.18257418458450 0.36514839177157 0.54772259953620 0.73029676840507
(etc.)

So now I have 10 error message boxes to clear. If the array was 1000 then I would have 1000 error message boxes to clear. Only one message per function call is needed for these types of cases.

My comments are admittedly a bit picky, but keep in mind that I think these can easily be addressed to improve the code and I gave it a fairly high rating.

James Tursa

[John Fuller]

Thanks for the comments. I am editing the code to accommodate notes 2-5. I'd like to leave the second euler angle constraint between -90 and 90 since reaching these bounds for Type 1 rotations causes singularities. Allowing angles past those bounds, while correct, makes detecting singularities more obscure.

[Stanley]

Works well. Agreed with others that some bounds on inputs, while perhaps obscuring internal code, make implementation into other code more cumbersome.
for instance, in this snippet, ea(1) must be in 0->360.

SpinCalc('EA321toDCM',ea,eps,0)

Alowing ea(1) = -10 degrees is easier.
Perhaps a switch to turn of bounds checking?

[John Fuller]

I am uploading an updated version with the boundaries revised on the first and third Euler angles. They must now be between -360 and 360. The bounds on the second Euler angles are going to stay since they prevent ambiguous Euler sets.

Thanks for the comments.

[Jørgen]

When using the DCM to EA323 I get some unexpected results.

This is my code:

>> R = [0 -1 0; 1 0 0; 0 0 1]
R =

     0 -1 0
     1 0 0
     0 0 1

>> SpinCalc2('DCMtoEA323',R,eps,1)

ans =

     0 0 0

------

The matrix R represents a simple 90 degrees rotation about the z-axis, but the Euler Angles are all zero.

[John Fuller]

Jørgen,
The 323 rotation you specified is a singular Euler set, because the middle angle is 0. Any type 2 rotation where the middle angle is 0 or 180 cannot be uniquely resolved by trying to translate Q, DCM, or EV back to Euler angles. If you plan on translating amongst orientations that are singular when expressed as Euler angles, I would advise you use strictly DCM, Q, or EV because they can uniquely define all orientations.

[Jørgen]

Okay, thank you for your feedback.

[Joel]

 

[James Tursa]

After comparing this submission with others in the FEX, I feel the need to bump my rating up.

[bryan]

I found this code to be very useful and it seems to work exactly as described. Thank you for your excellent contribution.