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
Lsq=: [: +/"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=: dyad define
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. In right handed coordinates, looking in the +z direction, a
NB. positive rotation (from x to y) is *clockwise* -- we must use
NB. -theta to rotate counter clockwise (or make the axis -z).
NB. Usage: u Qrot theta --> q  then  q Pr P --> P'
NB. u may be [3 n] array and theta [n] array.
Qrot=: dyad define
s=. cos -:y
f=. sin -:y
s ,&.|: f*x  NB. catenate under transpose
)

NB. For rotations about the y axis: QY angles --> rotation quaternions
QY=: [:>([:cos-:);0;0;~[:sin-:  NB. quaternion for rotation about y-axis
QZ=: [:>([:cos-:);0;0;[:sin-:   NB. quaternion for rotation about z-axis

NB. The conjuate of the quaternions: Qb q --> q_bar
Qb=: monad : '(1 _1 _1 _1)*"1 y'

NB. The inverse of quaternions: Qi q --> q^(-1)
Qi=: Qb % Lsq

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=: dyad : '}. x QM (0,y) QM Qb x'
Pr=: [: }. [ QM (0"_ , ]) 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=: dyad define
'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=: dyad define
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=: monad define
'w x y0 z'=. |: y  NB. seperate components of q
n=. #x2=. *:x      NB. n is number of quaternions
y2=. *:y0
z2=. *:z
a=. (1- +: y2+z2),(+:(x*y0)-(w*z)),(+:(x*z)+(w*y0))
b=. (+:(x*y0)+(w*z)),(1- +: x2+z2),(+:(y0*z)-(w*x))
c=. (+:(x*z)-(w*y0)),(+:(y0*z)+(w*x)),(1- +:x2+y2)
if. n=1 do.
  (3,3)$ a,b,c
else.
  |: (3,3,n)$ |. a,b,c
end.
)