Herman J. Woltring

02-21-1990, 10:40 PM

Dear Biomch-L readers,

Some readers on the list have shown interest in the mathematics and programming

details underlying the current joint attitude parametrization debate. For this

purpose, I'd like to refer them to a test package PRP FORTRAN which is available

to the readership by sending the request

SEND PRP FORTRAN BIOMCH-L

to LISTSERV@HEARN.EARN (or .BITNET); the request should be sent from the address

at which the Biomch-L listserver knows the requester. A FORTRAN-77 package is

then returned which accomodates all possible Cardanic conventions and the

proposed helical convention. While the PRP package was written for use under a

VAX/VMS system, only a few I/O statements in the package's testprogramme may

have to be modified for use in different environments.

Basically, the relation between a 3*3 orthonormal matrix Rijk (with |Rijk| = +1)

and the selected Cardanic angles can be described as follows:

Let Rijk = Ri(PHIi) * Rj(PHIj) * Rk(PHIk), where PHIi, PHIj, and PHIk denote

`planar' rotations, in the right-handed sense, about selected axes of a Carte-

sian co-ordinate system (1: x-axis, 2: y-axis, 3: z-axis). Defining ijk as

cyclic (i.e., 123, 231, or 312), the elements of Ri(PHIi) are defined as

rii = 1; rij = rji = rik = rki = 0

rjj = rkk = cos(PHIi); rkj = -rjk = sin(PHIi)

whence the elements r.. of Rijk can be expressed as:

rkj + rji = {1 + sin(PHIj)} sin(PHIi + PHIk)

rjj - rki = {1 + sin(PHIj)} cos(PHIi + PHIk)

rkj - rji = {1 - sin(PHIj)} sin(PHIi - PHIk)

rjj + rki = {1 - sin(PHIj)} cos(PHIi - PHIk)

This yields two unique pairs -- as per Codman's Paradox -- for PHIi and PHIk

unless |PHIj| = PI/2 -- gimbal-lock; in the latter case, two of the four rela-

tions above are zero. Note that all angles are modulo 2*PI. Given the choosen

PHIi and PHIk, PHIj can be derived from

rik = sin(PHIj)

{rii cos(PHIk) + rkk cos(PHIi)} - {rij sin(PHIk) + rjk sin(PHIi)} = cos(PHIj)

If ijk is anti-cyclic, i.e., 321, 213 or 132, the cyclic case can be obtained

via the intermediate transformation

Rijk(PHIi,PHIj,PHIk) = Rkji(-PHIk,-PHIj,-PHIi)'

Note that for other ijk, i.e., i=k, the n e u t r a l attitude corresponds to

gimbal-lock, with either PHIi + PHIk or PHIi - PHIk undefined. Interestingly,

Leonhard Euler's original publication `De Immutatione Coordinatarum, Caput IV,

Appendix de Superficiebus' in his `Introductio in Analysin Infinitorum' (Lausan-

ne, France 1748) provides such a case.

For the helical convention, the following formulae apply:

(Rijk - Rijk')/2 = sin(theta) A(N)

(Rijk + Rijk')/2 = cos(theta) I + {1-cos(theta)} N N'

THETA = theta N, N = (n1,n2,n3)', |N| = 1.

where A(N) is a skew-symmetric matrix with elements aji = -aij = nk and akk = 0,

for cyclic ijk. This system has a unique solution for 0 < theta < PI. If theta

= 0, no solution exists for N, but THETA is zero; for theta = PI, there are two

solutions, one THETA, the other -THETA, as one would expect for any periodic

solution modulo 2*PI.

Those who are interested in written material may get in touch with me directly,

and I'll be glad to send a mathematically oriented paper presented during a

1989 CAMARC Workshop on Clinical Protocols in Ancona/Italy.

Regards -- Herman J. Woltring, tel & fax +31.40.41 37 44

Some readers on the list have shown interest in the mathematics and programming

details underlying the current joint attitude parametrization debate. For this

purpose, I'd like to refer them to a test package PRP FORTRAN which is available

to the readership by sending the request

SEND PRP FORTRAN BIOMCH-L

to LISTSERV@HEARN.EARN (or .BITNET); the request should be sent from the address

at which the Biomch-L listserver knows the requester. A FORTRAN-77 package is

then returned which accomodates all possible Cardanic conventions and the

proposed helical convention. While the PRP package was written for use under a

VAX/VMS system, only a few I/O statements in the package's testprogramme may

have to be modified for use in different environments.

Basically, the relation between a 3*3 orthonormal matrix Rijk (with |Rijk| = +1)

and the selected Cardanic angles can be described as follows:

Let Rijk = Ri(PHIi) * Rj(PHIj) * Rk(PHIk), where PHIi, PHIj, and PHIk denote

`planar' rotations, in the right-handed sense, about selected axes of a Carte-

sian co-ordinate system (1: x-axis, 2: y-axis, 3: z-axis). Defining ijk as

cyclic (i.e., 123, 231, or 312), the elements of Ri(PHIi) are defined as

rii = 1; rij = rji = rik = rki = 0

rjj = rkk = cos(PHIi); rkj = -rjk = sin(PHIi)

whence the elements r.. of Rijk can be expressed as:

rkj + rji = {1 + sin(PHIj)} sin(PHIi + PHIk)

rjj - rki = {1 + sin(PHIj)} cos(PHIi + PHIk)

rkj - rji = {1 - sin(PHIj)} sin(PHIi - PHIk)

rjj + rki = {1 - sin(PHIj)} cos(PHIi - PHIk)

This yields two unique pairs -- as per Codman's Paradox -- for PHIi and PHIk

unless |PHIj| = PI/2 -- gimbal-lock; in the latter case, two of the four rela-

tions above are zero. Note that all angles are modulo 2*PI. Given the choosen

PHIi and PHIk, PHIj can be derived from

rik = sin(PHIj)

{rii cos(PHIk) + rkk cos(PHIi)} - {rij sin(PHIk) + rjk sin(PHIi)} = cos(PHIj)

If ijk is anti-cyclic, i.e., 321, 213 or 132, the cyclic case can be obtained

via the intermediate transformation

Rijk(PHIi,PHIj,PHIk) = Rkji(-PHIk,-PHIj,-PHIi)'

Note that for other ijk, i.e., i=k, the n e u t r a l attitude corresponds to

gimbal-lock, with either PHIi + PHIk or PHIi - PHIk undefined. Interestingly,

Leonhard Euler's original publication `De Immutatione Coordinatarum, Caput IV,

Appendix de Superficiebus' in his `Introductio in Analysin Infinitorum' (Lausan-

ne, France 1748) provides such a case.

For the helical convention, the following formulae apply:

(Rijk - Rijk')/2 = sin(theta) A(N)

(Rijk + Rijk')/2 = cos(theta) I + {1-cos(theta)} N N'

THETA = theta N, N = (n1,n2,n3)', |N| = 1.

where A(N) is a skew-symmetric matrix with elements aji = -aij = nk and akk = 0,

for cyclic ijk. This system has a unique solution for 0 < theta < PI. If theta

= 0, no solution exists for N, but THETA is zero; for theta = PI, there are two

solutions, one THETA, the other -THETA, as one would expect for any periodic

solution modulo 2*PI.

Those who are interested in written material may get in touch with me directly,

and I'll be glad to send a mathematically oriented paper presented during a

1989 CAMARC Workshop on Clinical Protocols in Ancona/Italy.

Regards -- Herman J. Woltring, tel & fax +31.40.41 37 44