hp41programs

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:

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