rotation around a line

format('v',15)

//rotation axis
u = [rand(), rand(), rand()];
u = u./ sqrt(u*u')

//angle
t = 2*%pi*rand(); //picked randomly

c = cos(t);
s = sin(t);

//translation vector
T= [rand(); rand(); rand()] //picked randomly

//rotation matrix
R = [c+u(1)*u(1)*(1-c), u(1)*u(2)*(1-c)-u(3)*s, u(1)*u(3)*(1-c)+u(2)*s;
 u(1)*u(2)*(1-c)+u(3)*s, c+u(2)*u(2)*(1-c), u(2)*u(3)*(1-c)-u(1)*s;
 u(1)*u(3)*(1-c)-u(2)*s, u(2)*u(3)*(1-c)+u(1)*s, c+u(3)*u(3)*(1-c)]

// FIRST TEST
//transformation matrix: rotation around the line u*t+T
// - the upper left 3x3 matrix is actually a rotation
// - the translation part depends on the line and the rotation
M = [eye(3,3),T;0,0,0,1]*[R,[0;0;0];0,0,0,1]*[eye(3,3),-T;0,0,0,1]

// SECOND TEST
//transformation matrix: rotation and post-rotation translation
// - the upper left 3x3 matrix is actually a rotation
// - the translation part is generated randomly (no link to the rotation)
//M = [eye(3,3),T;0,0,0,1]*[R,[0;0;0];0,0,0,1]

//identity
I = eye(4,4);

//system
A = I-M

//row reduction
S = rref(A)

//line direction
nu = S(1:3,3);
nu(3)=1;
nu = nu./ sqrt(nu'*nu)

//line position
nt = S(1:3,4)