*To*: mathgroup@smc.vnet.net*Subject*: [mg10440] RE: [mg10382] Re: Rotate3D bug solution*From*: Jean-Marie THOMAS <jmthomas@cybercable.tm.fr>*Date*: Tue, 13 Jan 1998 02:07:37 -0500

I followed with great interest all the discussions about Eulerian angles, I have not got much to say about them, (the name is the same in French litterature) but I remember having spent some time on the following code, which uses a generalization of Eulerian angles in unspecified dimension. The code uses Eulerian angles, and mostly shows the way in which these angles are defined: (* this module generates a random point on a sphere of center "center" and radius "radius". The algorithm uses generalized polar angles as random numbers and transforms the polar coordinates in cartesian co ordinates. *) randomPointOnASphere[center_,radius_]:=Module[ {dim,rp,rc,x}, dim=Length[center]; x[i_Integer,pc_]:=pc[[1]] Switch[i, 1,Product[Sin[pc[[j]] ],{j,2,dim}], _,Product[Sin[pc[[k]] ],{k,i+1,dim}] Cos[pc[[i]]] ]; rp=Flatten[ {radius,Random[Real,{0,2 Pi//N}],Table[Random[Real,{0,Pi//N}],{dim-2}]}]; rc=Table[x[i,rp],{i,dim}]; Return[rc+center] ] (* end *) Eulerian angles are here called polar angles. I'm not a mathematician, so the terminology might be wrong. Anyway this code returns random points on a sphere, random in the sense of uniform density on the surface of the sphere. Hope this helps, ----------------------------------------------- Jean-Marie THOMAS Conseil et Audit en Ingenierie de Calcul jmthomas@cybercable.tm.fr +33 (0)3 88 32 93 64 www.cybercable.tm.fr/~jmthomas ======================= -----Message d'origine----- De: John Sidles [SMTP:sidles@u.washington.edu] Date: lundi 12 janvier 1998 10:10 A: mathgroup@smc.vnet.net Objet: [mg10382] Re: Rotate3D bug solution In article <68csvb$aa4@smc.vnet.net>, Selwyn Hollis <shollis@peachnet.campus.mci.net> wrote: >Mark Evans wrote: > >> Paul Abbott wrote: >> > >> > The use of Eulerian angles for specifying rotations in 3D IS >> standard >> > (at least in maths and physics -- especially in quantum mechanics, >> > crystallography, and angular momemntum theory). Note that the >> eulerian >> > angle parametrization avoids the singularities that arise in other >> > parametrizations. >> > >> >> In the same sense, you could say that sea shells are standard legal >> tender if you live in a certain part of the world. >> >> Paul is right that there is nothing technically wrong with this kind >> of >> rotation. My point was that Mathematica packages should be written >> for >> a wider audience. It seems intuitive that the most common >> understanding of a rotation matrix is one that rotates sequentially >> about each of the three coordinate axes. The fact that Mathematica >> does not offer this rotation by default is a slip-up in my mind. > >Right on, Mark! > >I've never heard of Eulerian angles before encountering them in >Mathematica. Maybe they're the usual tricks-of-the-trade to a few >quantum physicists, but you'll be hard-pressed to find a reference to >them in any but the most esoteric mathematics literature. > Well, the latitude and longitude coordinates on the globe are (essentially) Eulerian coordinates -- so they're not *too* esoteric. Here's a very important and useful, yet simple, theorem which anyone working with rotations needs to know -- it explains why seemingly inequivalent conventions are actually precisely equivalent. Let $R(v)$ be a function which computes the three-by-three matrix associated with a three-vector $v$, where the direction of $v$ gives the axis of rotation, and the magnitude of $v$ gives the angle of rotation. (Exercise: program this useful utility function in Mathematica, together with its inverse function v(R) -- answer given at end!). Let $U$ be an arbitrary rotation matrix, with $U^{t}$ the matrix transpose of $U$. Since $U$ is a rotation, $U^{t}$ is the matrix inverse of $U$, i.e., $U^{t} U = I$. Then for any $U$ and $v$, here's the key theorem! U R(v) U^{t} = R(U v) To see how this clarifies the literature on rotations, suppose Textbook A defines Euler matrices in terms of rotation angles ${\theta, \phi, \psi}$ about fixed unit axes ${\hat{n}_1,\hat{n}_2,\hat{n}_3}$ as follows R(\theta, \phi, \psi) \edef R(\theta \hat{n}_1) R(\phi \hat{n}_2) R(\psi \hat{n}_3) \edef R_1 R_2 R_3 Now we use the above theorem to rewrite this in terms of *moving* axes. R(\theta,\phi,\psi) = R(\psi R1 R2 \hat{n_3}) R(\phi R1 \hat{n_2}) R(\theta \hat{n_1}) Cool! Its the same three angles, but now applied in the *opposite* order, and about moving instead of fixed axes. Yet the final matrix is the same. And it is perfectly reasonable for Textbook B to adopt this moving-axis convention to define Euler angles. Given all this ambiguity, I no longer use Euler angles when doing calculations involving rotations. It is just too easy to confuse the various signs and conventions! A much safer strategy, which I recommend, is to simply code the function R(v) and its inverse v(R) as utility routines in whatever language you prefer. You will find that these two functions, plus ordinary matrix multiplication, suffice for *any* calculation involving rotations. No more Euler angle torture! No more sine and cosine functions with obscure arguments! Here's some Mathematica code for R(v): rotationMatrix[x_] := Block[ {angle,xHat,pPar,pPerp,pOrthog}, angle = Sqrt[x.x]//N; If[angle<10^-6,Return[id]]; xHat = x/angle; pPar = Outer[Times,xHat,xHat]; pPerp =DiagonalMatrix[{1.0,1.0,1.0}] - pPar; pOrthog = { {0.0,xHat[[3]],-xHat[[2]]}, {-xHat[[3]],0.0,xHat[[1]]}, {xHat[[2]],-xHat[[1]],0.0} }; pPar + Cos[angle]* pPerp + Sin[angle]*pOrthog ] It is left as an exercise to (a) figure out how the above works, and (b) code the inverse function! Coding the inverse function is quite a nontrivial exercise. Hint: (a) determine the cosine() from the trace of R, then determine the sine() and the the axis of rotation from the antisymmetric part of R. This works for all rotations *except* for angles near Pi, for which the axis of rotation should be set to the (unique) eigenvector of the symmetric part of R that has unit eigenvalue (this is because sin(Pi) = 0). A final refinement: the above algorithm assumes that R is orthogonal. But what if R is just *close* to orthogonal, but has some accumulated numerical imprecision? Even worse, you can bet that some future user of your v(R) routine will hand it an R matrix that is grossly non-orthonormal! The right thing to do, therefore, is to condition R before calculating v, by (a) calculating the singular value decomposition for R, then (b) adjusting all the singular values to have unit magnitude. This will yield a cleaned-up exactly orthogonal R which (formally) is the orthogonal matrix that is closest to the input matrix in the least-mean squares sense. Better issue an error message if the singular values are far from unity -- because this indicates abuse of the inverse routine! The Mathematica routines SingularValue[] and Eigensystem[] are well suited to the above tasks -- it's tougher in "C". I wrote this up at length because I have a student who is generating animations in Mathematica and POV-Ray -- he might as well learn to do things the easy way! Happy rotating ... JAS