I don't get it. If a and b are length-2 lists, the result has shape 4 2. Is that right?
Henry Rich > -----Original Message----- > From: [EMAIL PROTECTED] > [mailto:[EMAIL PROTECTED] On Behalf Of Kip Murray > Sent: Friday, April 18, 2008 9:24 PM > To: Programming forum > Subject: Re: [Jprogramming] Quaternions > > The quaternion > > a_0 + i a_1 + j b_0 + k b_1 > > (underscores indicate subscripts) can be represented by the matrix > > (a,b) ,: (-+b),+a NB. b is b_0 + i b_1, monadic + is conjugate > > Then the determinant gives the square of the magnitude, matrix > multiplication gives the product, matrix inverse the > reciprocal, matrix > addition the sum. > > Kip Murray > > > On Sat, 19 Apr 2008, Chris Burke wrote: > > See also: > > open '~system/packages/math/quatern.ijs' > > J. Patrick Harrington wrote: > > On Fri, 18 Apr 2008, Dominic Genest wrote: > >> Hi, > >> > >> Are there phrases for quaternions ? > >> Mainly, I need to convert from quaternions to rotation matrices. > >> > >> Thanks > >> > >> Dom > >> > > I've used quaternions on occasion, and find two sets of routines > > that differ depending on the way you store an array of 4-element > > quaternions: (4 n) or (n 4). Here are both. I think they work but > > no warranty... In the following, > > > > mm=: +/ . * NB. matrix multiply > > > > NB. Quaternion maths. NB. Four-element quaternion is > scalar followed by > > 3-vector > > NB. Should work on (n 4) arrays of quaternions. > > > > L=: +/"1 &. (*:"_) NB. Gives the lengths of vector arrays > > L2=: [: +/"1 *:"_ NB. squares of lengths > > dot=: [: +/"1 * NB. dot products of arrays of vectors > > pm1=: 1&|."1 NB. 1st cyclic permutation > > pm2=: 2&|."1 NB. 2nd cyclic permutation > > pm12=: 4 : '(pm1 x)*pm2 y' > > cross=: pm12 - pm12~ NB. cross products of vectors > > NB. cross=: pm12~ - pm12 NB. for left-handed coordinate system > > > > NB. Multiply quaternions: q1 QM q2 --> q3 > > QM=: 4 : 0 > > s1=. 0{"1 x > > v1=. }."1 x > > s2=. 0{"1 y > > v2=. }."1 y > > s=. (s1*s2) - v1 dot v2 > > v=. (s1*v2)+(s2*v1)+(v1 cross v2) > > s ,&.|: v NB. catenate under transpose > > ) > > > > NB. Quaternion for rotation of points by theta about unit vector u > > NB. (Rotates *points*; use -theta to rotate coordinate axes) > > NB. Usage: u Qrot theta --> q then q Pr P --> P' > > Qrot=: 4 : 0 > > u=. x > > s=. cos -:y > > f=. sin -:y > > s ,&.|: (f*u) NB. catenate under transpose > > ) > > > > NB. The conjuate of the quaternions: Qb q --> q_bar > > Qb=: 3 : '(1 _1 _1 _1)*"1 y' > > > > NB. The inverse of quaternions: Qi q --> q^(-1) > > Qi=: Qb % L2 > > > > NB. Rotate point using quaternion q > > NB. Usage: q Pr (x y z) --> (x' y' z') > > Pr=: 4 : 0 > > qb=. Qb q=. x > > }."1 q QM (0 ,&.|: y) QM qb > > ) > > > > norm=: ] % L NB. normalize to unit length > > > > NB. Rotate points P (P[n,3]) about some vector u NB. by > angle theta: (0 > > 1 0; 0.25p1) Qrtn P --> P' > > NB. (theta in radians; u need not be a unit vector) > > Qrtn=: 4 : 0 > > 'u theta'=. x > > ((norm u) Qrot theta) & Pr"1 y > > ) > > > > NB. Usage: DC Trans P --> P', where DC = direction > cosines of z-axis > > NB. P point in scattering frame, P' point in fixed system > > Trans=: 4 : 0 > > NB. rotation axis: u=, (x y z) cross (0 0 1), normalized > > u=. norm (1 0 2){"1 (_1 1 0)*"1 x theta=. -acos 2{"1 x > NB. rotate point > > back to fixed axes > > q=. u Qrot theta NB. the rotation quaternion > > q Pr y NB. rotate P in scattering frame to > P' in fixed > > ) > > > > NB. Form 3x3 rotation matrix which effects the same > transformation as > > NB. the quaternion q. Thus, Quat_2_Mat q --> M ; M mm P --> P' > > Quat_2_Mat=: 3 : 0 > > 'w x yy z'=. |: y NB. seperate components of q > > n=. #x2=. *:x NB. n is number of quaternions > > y2=. *:yy > > z2=. *:z > > a=. (1- +: y2+z2),(+:(x*yy)-(w*z)),(+:(x*z)+(w*yy)) > > b=. (+:(x*yy)+(w*z)),(1- +: x2+z2),(+:(yy*z)-(w*x)) > > c=. (+:(x*z)-(w*yy)),(+:(yy*z)+(w*x)),(1- +:x2+y2) > > if. n=1 do. > > (3,3)$ a,b,c > > else. > > |: (3,3,n)$ |. a,b,c > > end. > > ) > > > > ================================================================ > > NB. Quaternion maths. NB. Four-element quaternion is > scalar followed by > > 3-vector > > NB. Should work on (4 n) arrays of quaternions. > > > > L=: +/ &. (*:"_) NB. Gives the lengths of vector arrays > > L2=: [: +/ *:"_ NB. squares of lengths > > dot=: [: +/ * NB. dot products of arrays of vectors: > [4 n] dot [4 n] > > dot2=: [: +/ ]* [: ,[ NB. 1 vector dotted into array: [4 > 1] dot2 [4 n] > > NB. pm1=: 1&|. NB. 1st cyclic permutation > > NB. pm2=: 2&|. NB. 2nd cyclic permutation > > NB. pm12=: 4 : '(pm1 x)*pm2 y' > > NB. cross=: pm12 - pm12~ NB. cross products of vectors > > NB. cross=: pm12~ - pm12 NB. for left-handed coordinate system > > cross=: ((1: |.[)*(_1: |. ]))-((_1: |.[)*(1:|.])) > > > > NB. Multiply quaternions: q1 QM q2 --> q3 > > QM=: 4 : 0 > > s1=. 0{ x > > v1=. }. x > > s2=. 0{ y > > v2=. }. y > > s=. (s1*s2) - v1 dot v2 > > s,(s1*"1 v2)+(s2*"1 v1)+(v1 cross v2) > > ) > > > > NB. Quaternion for rotation of points by theta about unit vector u > > NB. (Rotates *points*; use -theta to rotate coordinate axes) > > NB. Usage: u Qrot theta --> q then q Pr P --> P' > > NB. u may be [3 n] array and theta [n] array. > > Qrot=: 4 : 0 > > s=. cos -: y > > f=. sin -: y > > s, (f*"1 x) ) > > > > NB. The conjuate of the quaternions: Qb q --> q_bar > > Qb=: 3 : '(1 _1 _1 _1)* y' > > > > NB. The inverse of quaternions: Qi q --> q^(-1) > > Qi=: Qb %"1 L2 > > > > NB. Rotate point using quaternion q > > NB. Usage: q Pr P --> P' where q is [4 n] array of quaternions, > > NB. P is [3 n] array of point coordinates, and P' also [3 n] > > NB. Pr=: 4 : '}. x QM (0,y) QM Qb x' > > Pr=: [: }. [ QM (0"_ , ]) QM [: Qb [ > > > > norm=: ] %"1 L NB. normalize to unit length > > > > NB. Rotate points P (P[n,3]) about some vector u NB. by > angle theta: (0 > > 1 0; 0.25p1) Qrtn P --> P' > > NB. (theta in radians; u need not be a unit vector) > > Qrtn=: 4 : 0 > > 'u theta'=. x > > ((norm u) Qrot theta) & Pr y > > ) > > > > NB. Usage: DC Trans P --> P', where DC = direction > cosines of z-axis > > NB. P point in scattering frame, P' point in fixed system > > Trans=: 4 : 0 > > NB. rotation axis: u=, (x y z) cross (0 0 1), normalized > > u=. norm (1 0 2){ (_1 1 0)* x theta=. -acos 2{ x NB. > rotate point > > back to fixed axes > > q=. u Qrot theta NB. the rotation quaternion > > q Pr y NB. rotate P in scattering frame > to P' in fixed > > ) > > > > NB. Form 3x3 rotation matrix which effects the same > transformation as > > NB. the quaternion q. Thus, Quat_2_Mat q --> M ; M mm P --> P' > > Quat_2_Mat=: 3 : 0 > > 'w x yy z'=. |: y NB. seperate components of q > > n=. #x2=. *:x NB. n is number of quaternions > > y2=. *:yy > > z2=. *:z > > a=. (1- +: y2+z2), (+:(x*yy)-(w*z)), (+:(x*z)+(w*yy)) > > b=. (+:(x*yy)+(w*z)), (1 - +: x2+z2), (+:(yy*z)-(w*x)) > > c=. (+:(x*z)-(w*yy)), (+:(yy*z)+(w*x)), (1- +:x2+y2) > > if. n=1 do. > > (3,3)$ a,b,c > > else. > > |: (3,3,n)$ |. a,b,c > > end. > > ) > > > ---------------------------------------------------------------------- > For information about J forums see > http://www.jsoftware.com/forums.htm > > > Kip Murray > > [EMAIL PROTECTED] > http://www.math.uh.edu/~km > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
