Maple Example of 3-D Rotations

Here's an example of using Maple to combine two 90° rotations, one about the X-axis, the other about the Y-axis. The code is here.

      > read `rot_maple.code`;
        [dotprod]
        [crossprod]
	Warning: new definition for   norm
        [norm]
        [scalarmul]
	proc(expr) ... end

	rotate :=
    	proc(a,t,p)
    	local a2;
        a2 := normalize(a);
        evalm(cos(t)*p+(1-cos(t))*dotprod(a2,p)*a2+sin(t)*crossprod(a2,p))
    	end

	normalize := proc(a) local l; l := evalf(norm(a,2)); evalm(a/l) end

        [ 1  0  0 ]
        [         ]
        ii := [ 0  1  0 ]
        [         ]
        [ 0  0  1 ]

	rotmat := proc(a,t)
        local a2;
        a2 := normalize(a);
        evalm(cos(t)*ii+extprod(a2,a2)*(1-cos(t))+sin(t)*crossmat(a2))
        end

	rotmat2 := proc()
        local a,a1,a2,a3,c,s;
        a[1] := a1;
        a[2] := a2;
        a[3] := a3;
        evalm(c*ii+extprod(a,a)*(1-c)+s*crossmat(a))
        end

	extprod := proc(a,b)
        local m,i,j;
        m := array(1 .. 3,1 .. 3);
        for i to 3 do  for j to 3 do  m[i,j] := a[i]*b[j] od od;
        op(m)
        end

	crossmat :=
  	proc(a)
  	local m;
      	m := array(1 .. 3,1 .. 3,[[0,-a[3],a[2]],[a[3],0,-a[1]],[-a[2],a[1],0]]);
      	op(m)
  	end
	----------------------------------------------------------------
	Find the rotation matrix for an axis of (1,0,0) and an angle of Pi/2:

      > r1:=evalf(rotmat([1,0,0],Pi/2));
        [ 1.   0   0  ]
        [             ]
        r1 := [  0   0  -1. ]
        [             ]
        [  0  1.   0  ]
	----------------------------------------------------------------
      > r2:=evalf(rotmat([0,1,0],Pi/2));
        [  0    0  1. ]
        [             ]
        r2 := [  0   1.   0 ]
        [             ]
        [ -1.   0   0 ]
	----------------------------------------------------------------
      > r12:=multiply(r1,r2);
        [  0   0  1. ]
        [            ]
        r12 := [ 1.   0   0 ]
        [            ]
        [  0  1.   0 ]
	----------------------------------------------------------------
      > r21:=multiply(r2,r1);
        [  0   1.   0  ]
        [              ]
        r21 := [  0    0  -1. ]
        [              ]
        [ -1.   0   0  ]
	----------------------------------------------------------------
      > eigenvects(r12);

  	[ - .4999999997 - .8660254037 I, 1,
      	{[  - .7382716585 + .3487429165 I, .06711560405 - .8137334705 I,
        .6711560548 + .4649905518 I ]}                               ],
      	[ - .4999999997 + .8660254037 I, 1,
        {[  - .7382716585 - .3487429165 I, .06711560405 + .8137334705 I,
        .6711560548 - .4649905518 I ]}                               ],
      	[1.000000000, 1, {[ -.5773502687, -.5773502691, -.5773502691 ]}]

	This represents a 120 degree rotation about the axis 
	[ -.5773502687, -.5773502691, -.5773502691 ] .

	----------------------------------------------------------------
      > eigenvects(r21);

   	[ - .4999999997 - .8660254037 I, 1,
       	{[ .7382716585 - .3487429165 I,  - .6711560548 - .4649905518 I,
        .06711560405 - .8137334705 I ]}                             ],
       	[ - .4999999997 + .8660254037 I, 1,
        {[ .7382716585 + .3487429165 I,  - .6711560548 + .4649905518 I,
        .06711560405 + .8137334705 I ]}                             ],
       	[1.000000000, 1, {[ .5773502687, .5773502691, -.5773502691 ]}]

	This represents a 120 degree rotation about the axis 
	[ .5773502687, .5773502691, -.5773502691 ], which is different, so the
	order does matter.