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
s=. cos -:y
f=. sin -:y
s ,&.|: f*x  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 x
}."1 x 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 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.
)