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 MT = MT M = I    where I  is the Identity matrix

M may be written  M = ( G12 G13 ..... G1n ) ( G23 G24 ..... G2n ) ........  Gn-1,n x U

where U is a diagonal matrix whose diagonal elements are +/-1
and   Gi,j  are Givens rotations matrices of an angle µi,j in the plane ( Oxi , Oxj )

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 GT12  to zero the element  a12 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 Rn2+1 thru R(3n2-n)/2  and the diagonal matrix U has replaced M in registers R01 thru Rn2

-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 = a11   •  Rn+1 = a12  .....................  •  Rn2-n+1 = a1n
•  R02 = a21   •  Rn+2 = a22  .....................  •  Rn2-n+2 = a2n
........................................................................................

•  Rnn = an1   •  R2n = an2    .....................   •  Rn2 = ann

>>>>  When the program stops,   Rn2+1 = µ12 , ................. , R(3n2-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 = n2+1 ,  eee = ( 3n2 - 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 = µ12 = 41°         R20 = µ23 = 34°       R22 = µ34  = -49°           ( in DEG mode )
R18 = µ13 = 67°         R21 = µ24  = 48°
R19 = µ14 = -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" )

•  Rn2+1 = µ12 , ................. , •  R(3n2-n)/2 = µn-1,n

R01 = a11     Rn+1 = a12  .....................    Rn2-n+1 = a1n
R02 = a21     Rn+2 = a22  .....................    Rn2-n+2 = a2n
>>>>  When the program stops,                       ........................................................................................

Rnn = an1     R2n = an2    .....................     Rn2 = ann

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 = µ12 = 41°         R20 = µ23 = 34°       R22 = µ34  = -49°
R18 = µ13 = 67°         R21 = µ24  = 48°
R19 = µ14 = -55°

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

and we get

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

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 Gi,j as above.

Data Registers:               R00 = n > 1

R01 = a11     Rn+1 = a12  .....................    Rn2-n+1 = a1n
R02 = a21     Rn+2 = a22  .....................    Rn2-n+2 = a2n
>>>>  When the program stops,                       ........................................................................................

Rnn = an1     R2n = an2    .....................     Rn2 = ann
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:

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