Euler Angles

Generalized Euler Angles & Orthogonal Matrices for the HP-41

Overview

1°) Orthogonal Matrix >>> Generalized Euler Angles
2°) Generalized Euler Angles >>> Rotation Matrix
3°) Random Orthogonal Matrices

1°) Generalized Euler Angles

M is an orthogonal Matrix if its transpose is equal to its inverse: M M^T = M^T M = I where I is the Identity matrix

M may be written M = ( G₁₂ G₁₃ ..... G_1n ) ( G₂₃ G₂₄ ..... G_2n ) ........ G_n-1,n x U

where U is a diagonal matrix whose diagonal elements are +/-1
and G_i,j are Givens rotations matrices of an angle µ_i,j in the plane ( Ox_i , Ox_j )

                                      I      0        0       0        0
                                     0   cos µ_i,j   0   -sin µ_i,j    0
                                     0      0        I        0        0              1 <= i < j <= n
                                     0   sin µ_i,j   0     cos µ_i,j   0
                                     0      0       0        0        I

-The n(n-1)/2 angles µ_i,j are generalized Euler angles ( other definitions are possible ).

-"GEUA" finds the successive µ_i,j as follows:
M is multiplied by a Givens matrix G^T₁₂ to zero the element a₁₂ of the resulting matrix.
-The process is repeated until all the off-diagonal elements are 0.

-Since each matrix multiplication only modifies 2 rows - namely row i & row j - there is no need to call a subroutine to perform these products !

-Finally, the angles are in registers Rn²+1 thru R(3n²-n)/2 and the diagonal matrix U has replaced M in registers R01 thru Rn²

-If all the diagonal elements are +1 , M is a rotation matrix ( det M = 1 )
-If there is an odd number of (-1) , M is orthogonal but it's not a rotation matrix ( product of rotation and reflexion )
-If the off-diagonal elements are significantly different from zero ( except roundoff-errors as usual ) M was not orthogonal !

Data Registers: • R00 = n > 1 ( Registers R00 thru Rn^2 are to be initialized before executing "GEUA" )

                                • R01 = a₁₁   • Rn+1 = a₁₂ ..................... • Rn²-n+1 = a_1n
                                • R02 = a₂₁   • Rn+2 = a₂₂ ..................... • Rn²-n+2 = a_2n
                                    ........................................................................................

• Rnn = a_n1 • R2n = a_n2 ..................... • Rn² = a_nn

>>>> When the program stops, Rn²+1 = µ₁₂ , ................. , R(3n²-n)/2 = µ_n-1,n

Flags: /
Subroutines: /

01 LBL "GEUA"
02 RCL 00
03 1
04 ST+ Y
05 STO N
06 STO O
07 %
08 RCL 00
09 X^2
10 ST+ O
11 LASTX
12 -
13 DSE X
14 +
15 E3
16 /
17 1
18 +

19 STO M
20 LBL 01
21 RCL M
22 INT
23 RCL 00
24 MOD
25 LASTX
26 E3
27 /
28 +
29 ST+ N
30 LBL 02
31 RCL IND N
32 RCL IND M
33 R-P
34 X<>Y
35 STO IND O
36 1

37 ST+ O
38 P-R
39 STO P
40 X<>Y
41 STO Q
42 RCL 00
43 LASTX
44 %
45 X<>Y
46 X^2
47 +
48 PI
49 INT
50 10^X
51 /
52 RCL M
53 INT
54 +

55 RCL N
56 LBL 03
57 RCL IND X
58 RCL P
59 *
60 X<> IND Y
61 RCL Q
62 *
63 X<> IND Z
64 ST* L
65 LASTX
66 ST- IND Z
67 CLX
68 RCL P
69 *
70 ST+ IND Z
71 CLX
72 RCL 00

73 +
74 ISG Y
75 GTO 03
76 ISG N
77 GTO 02
78 ISG M
79 GTO 01
80 R^
81 INT
82 2
83 +
84 RCL O
85 DSE X
86 E3
87 /
88 +
89 CLA
90 END

( 135 bytes / SIZE (3n^2-n+2)/2 )

STACK	INPUT	OUTPUT
X	/	bbb.eee

Where bbb.eee is the control number of the n(n-1)/2 angles: bbb = n²+1 , eee = ( 3n² - n )/2

Example: n = 4 STO 00 and the orthogonal matrix is

               0.169141110   -0.444368761    -0.782087292    0.402823977        R01   R05   R09   R13
     M =   0.147032124     0.348746435   -0.549862104   -0.744586560   =   R02   R06   R10   R14
               0.527979893     0.706557116   -0.045724360    0.468960080        R03   R07   R11   R15
             -0.819152044     0.426250361   -0.289655687    0.251793847         R04   R08   R12   R16

XEQ "GEUA" >>>> 17.022 ---Execution time = 26s---

   and     R17 = µ₁₂ = 41°         R20 = µ₂₃ = 34°       R22 = µ₃₄ = -49°           ( in DEG mode )
             R18 = µ₁₃ = 67°         R21 = µ₂₄ = 48°
             R19 = µ₁₄ = -55°

Notes:

-For n = 10 , execution time is about 5mn20s
-Since there are only 319 registers, n cannot exceed 14.
-However, after zeroing a subdiagonal element, it is no more used in the rest of the calculations.
-So we can use the corresponding registers to store the Euler angles, thus overcoming this limitation:
-The following variant can deal with 17x17 orthogonal matrices !

01 LBL "GEUA"
02 RCL 00
03 1
04 ST+ Y
05 STO N
06 %
07 RCL 00
08 X^2
09 LASTX
10 -
11 DSE X
12 +
13 E3
14 /
15 1
16 +

17 STO M
18 LBL 01
19 RCL M
20 INT
21 RCL 00
22 MOD
23 LASTX
24 E3
25 /
26 +
27 ST+ N
28 LBL 02
29 RCL IND N
30 RCL IND M
31 R-P
32 X<>Y

33 STO O
34 1
35 P-R
36 STO P
37 X<>Y
38 STO Q
39 RCL 00
40 LASTX
41 %
42 X<>Y
43 X^2
44 +
45 PI
46 INT
47 10^X
48 /

49 RCL M
50 INT
51 +
52 RCL N
53 LBL 03
54 RCL IND X
55 RCL P
56 *
57 X<> IND Y
58 RCL Q
59 *
60 X<> IND Z
61 ST* L
62 LASTX
63 ST- IND Z
64 CLX

65 RCL P
66 *
67 ST+ IND Z
68 CLX
69 RCL 00
70 +
71 ISG Y
72 GTO 03
73 X<> O
74 STO IND N
75 ISG N
76 GTO 02
77 ISG M
78 GTO 01
79 CLA
80 END

( 121 bytes / SIZE n^2+1 )

STACK	INPUT	OUTPUT
X	/	µ_n-1,n

Example: With the same example, in DEG mode

XEQ "GEUA" >>>> -49 ---Execution time = 25s---

and we find in registers R01 to R16:

   R01   R05   R09   R13         1      0      0     0
   R02   R06   R10   R14        41     1      0     0
   R03   R07   R11   R15   =   67    34     1     0        the Euler angles are written in red.
   R04   R08   R12   R16       -55    48 -49    1

Notes:

-As usual, there are small roudoff-errors.
-Like the 1st version, "GEUA" works in all angular modes.
-In the 3-dimensional cases, there are many conventions for the Euler angles,
and those given by "GEUA" are different from the most commonly used.

2°) Rotation Matrix from Generalized Euler Angles

-Given n(n-1)/2 Euler angles, this program creates the corresponding rotation matrix.

Data Registers: • R00 = n > 1 ( Registers R00 and Rn^2+1 thru R(3n^2-n+2)/2 are to be initialized before executing "ROTM" )

• Rn²+1 = µ₁₂ , ................. , • R(3n²-n)/2 = µ_n-1,n

                                                                           R01 = a₁₁     Rn+1 = a₁₂ .....................    Rn²-n+1 = a_1n
                                                                           R02 = a₂₁     Rn+2 = a₂₂ .....................    Rn²-n+2 = a_2n
   >>>> When the program stops,                       ........................................................................................

Rnn = a_n1 R2n = a_n2 ..................... Rn² = a_nn

Flags: /
Subroutines: /

01 LBL "ROTM"
02 RCL 00
03 X^2
04 STO O
05   E3
06 /
07 ISG X
08 CLRGX
09 RCL 00
10 1
11 +
12   E5
13 /
14 +
15 SIGN
16 ST+ O
17 LBL 00
18 STO IND L
19 ISG L

20 GTO 00
21 RCL 00
22 1
23 %
24 X<>Y
25 X^2
26 LASTX
27 -
28 +
29 E3
30 /
31 1
32 +
33 STO M
34 LBL 01
35 RCL M
36 RCL 00
37 .1
38 %

39 +
40 +
41 STO N
42 LBL 02
43 RCL IND O
44 1
45 ST+ O
46 P-R
47 STO P
48 X<>Y
49 STO Q
50 RCL M
51 INT
52 LASTX
53 RCL 00
54 /
55 INT
56 RCL 00
57 ST* Y

58 +
59 PI
60 INT
61 10^X
62 /
63 +
64 RCL N
65 LBL 03
66 RCL IND X
67 RCL P
68 *
69 X<> IND Y
70 RCL Q
71 *
72 X<> IND Z
73 ST* L
74 LASTX
75 ST- IND Z
76 CLX

77 RCL P
78 *
79 ST+ IND Z
80 CLX
81 SIGN
82 +
83 ISG Y
84 GTO 03
85 ISG N
86 GTO 02
87 ISG M
88 GTO 01
89 FRC
90 1
91 +
92 CLA
93 END

( 138 bytes / SIZE (3n^2-n+2)/2 )

STACK	INPUT	OUTPUT
X	/	bbb.eeeii

Where bbb.eeeii is the control number of the rotation matrix: bbb = 1 , eee = n^2 , ii = nn

Example: n = 4 STO 00 and the generalized Euler angles are:

             R17 = µ₁₂ = 41°         R20 = µ₂₃ = 34°       R22 = µ₃₄ = -49°
             R18 = µ₁₃ = 67°         R21 = µ₂₄ = 48°
             R19 = µ₁₄ = -55°

XEQ "ROTM" >>>> 1.01604 ---Execution time = 27s---

and we get

Notes:

-For n = 10 , execution time is about 6mn08s
-Since there are only 319 registers, n cannot exceed 14.

3°) Random Orthogonal Matrices

-"ROM" chooses at random a diagonal matrix D whose diagonal elements are +/-1 and n(n-1)/2 angles µ_i,j
-It then performs the matrix product of D and G_i,j as above.

Data Registers: R00 = n > 1

Rnn = a_n1 R2n = a_n2 ..................... Rn² = a_nn
Flags: /
Subroutines: /

01 LBL "ROM"
02 STO O
03 X<>Y
04 STO 00
05 X^2
06   E3
07 /
08 ISG X
09 CLRGX
10 RCL 00
11 1
12 +
13   E5
14 /
15 +
16 ENTER^
17 LBL 00
18 CLX
19 RCL O
20 R-D

21 FRC
22 STO O
23 ST+ X
24 1
25 -
26 SIGN
27 STO IND Y
28 ISG Y
29 GTO 00
30 RCL 00
31 DSE X
32 E3
33 /
34 ISG X
35 STO M
36 LBL 01
37 RCL M
38 1.001
39 +
40 STO N

41 LBL 02
42 RCL O
43 R-D
44 FRC
45 STO O
46 1
47 CHS
48 ACOS
49 ST+ X
50 *
51 1
52 P-R
53 STO P
54 X<>Y
55 STO Q
56 RCL 00
57 LASTX
58 %
59 X<>Y
60 X^2

61 +
62 PI
63 INT
64 10^X
65 /
66 RCL M
67 INT
68 +
69 RCL N
70 LBL 03
71 RCL IND Y
72 RCL P
73 *
74 X<> IND Z
75 RCL Q
76 *
77 X<> IND Y
78 ST* L
79 LASTX
80 ST- IND T

81 CLX
82 RCL P
83 *
84 ST+ IND Y
85 CLX
86 RCL 00
87 +
88 ISG Y
89 GTO 03
90 ISG N
91 GTO 02
92 ISG M
93 GTO 01
94 SIGN
95 X<>Y
96 FRC
97 +
98 CLA
99 END

( 150 bytes / SIZE n^2+1 )

STACK	INPUTS	OUTPUTS
Y	n > 1	/
X	r	bbb.eeeii

Where bbb.eeeii is the control number of the orthogonal matrix: bbb = 1 , eee = n^2 , ii = nn
and r is a random seed

Example:

3 ENTER^
1 XEQ "ROM" >>> 1.00903 ---Execution time = 14s---

and we find in registers R01 thru R09:

     R01    R04    R07          0.836705341     0.426213209    -0.343898929
     R02    R05    R08 =      0.547653282    -0.651439385    0.525073910
     R03    R06    R09         -0.000235871    -0.627669522   -0.778479875

-If you execute then "GEUA" you get 10.012

and the generalized Euler angles are R10 = 33°20616782 , R11 = -0°013514413 , R12 = -141°1216061

-Since R01 = R05 = R09 = 1 , M is a rotation matrix ( det M = 1 )

Remarks:

-With n = 3 and r = 2 , you'll get an orthogonal matrix that is not a rotation matrix ( det M = -1 )
-If n = 10, execution time is about 6mn37s.
-Of course, this program may be simplified if you have an M-code function RAND and it will run faster.
-Since "ROM" only uses n^2+1 registers, you can create a random orthogonal matrix of order 17.

-The method that is used here to choose the angles is not perfect !
-A better one is described in the reference below...

Reference:

[1] T.W. Anderson, I. Olkin , L.G. Underhill - "Generation of Random Orthogonal Matrices" - TR 216

hp41programs

Generalized Euler Angles & Orthogonal Matrices for the HP-41