NB. Rotate n vectors v about axes u by angles w.
NB. Makes use of this quaternion-based formula: 
NB. with  u'= sin(w/2)*u   and   t= 2(u' X v)
NB. v' = v + [t * cos(w/2)] + (u' X t)
NB. Dimensions of arrays v and u are $V --> (n 3). 
NB. The axes u must be normalized to unit length. 
NB. Usage:   v qRot (u;w) --> v'
NB. Timing: 1e6 rotations in 0.26 sec on Mac mini

qRot=: dyad define
'u w'=. y
NB. u=. norm u  NB. Add this if u's are not normalized.
u=. u*sin-:w
t=. +: u cross x
x+ (t*cos-:w)+ u cross t
)

L=: +/"1 &. (*:"_)         NB. Lengths of vectors 
norm=: ] % L               NB. normalize to unit length
NB. cross products of vectors: v1 cross v2 --> v1 X v2
cross=: (([: 1&|."1 [)* [: 2&|."1 ])- (([: 1&|."1 [)* [: 2&|."1 ])~

--------------------------------------------------------------
NB. The axes u must be normalized to unit length. 
NB. Usage:   v qRot3 (u;w) --> v'  (v' rotated v's)
NB. This version can handle cases where any inputs
NB. v,u & w may be single values or arrays.

qRot3=: dyad define
'u w'=. y
NB. u=. norm u        NB. Add this if u's are not normalized.
if. 2>#$u do. u=. ($x)$u end.  NB. for single axis, clone it
if. (2>#$u) *. (1<#w) do.      NB. for single u & v ...
   u=. u */~ sin-:w            NB. this will produce (#w) u's
else. 
   u=. u * sin-:w       NB. if there's a u for each w, "*" OK
end.
NB. for one axis and one vector v, but multiple w's, clone v: 
if. 2>#$x do. x=. ($u)$x end.  
t=. +: u cross x
x+ (t*cos-:w)+ u cross t
)