# Copyright (c) Erich Kaltofen 1997.

with(linalg):
# rotation parallel to the x-y plane, around the z-axis	
Txy := matrix(3,3,[ cos(theta), -sin(theta), 0,
                    sin(theta),  cos(theta), 0,
                             0,           0, 1 ]);
rotxy := proc(vec, angle)
   # vec to be rotated by angle
   evalm(subs(theta = angle, evalm(Txy)) &* vec);
end;


# rotation parallel to the x-z plane
Txz := matrix(3,3,[ cos(eta),  0, -sin(eta),
                           0,  1,         0,
                    sin(eta),  0,  cos(eta) ]);
rotxz := proc(vec, angle)
   # vec to be rotated by angle
   evalm(subs(eta = angle, evalm(Txz)) &* vec);
end;

# rotation parallel to the y-z plane
Tyz := matrix(3,3,[ 1,  0,       0,
                    0,  cos(xi), -sin(xi),
                    0,  sin(xi),  cos(xi) ]);
rotyz := proc(vec, angle)
   # vec to be rotated by angle
   evalm(subs(xi = angle, evalm(Tyz)) &* vec);
end;

# check if rotated vector changes its length
v := vector([1,2,3]);
radsim(norm(v, 2));

# radsimp is radical simplification
radsimp(norm(rotxz(v, Pi/8), 2));

# make four animations showing the different rotations

# in order to view an animation, you must set the Plot Display to
#     a separate Window via the Options menu and then Play the plot
#     by clicking the right mouse button on the plot and selecting
#     the Animation submenu

with(plots):
k := 8: # number of frames
L := []:
for i from 0 to k do
plt[i] := polygonplot3d(
              [rotxy(vector([1,0,0]), i*2*Pi/k),
               rotxy(vector([0,2,0]), i*2*Pi/k),
               rotxy(vector([0,0,3]), i*2*Pi/k)
              ],
              axes = normal,
              style = patch,
              orientation = [-25,75]):
L := [op(L), plt[i]]:
od:
display(L, insequence=true);

k := 8: # number of frames
L := []:
for i from 0 to k do
plt[i] := polygonplot3d(
              [rotxz(vector([1,0,0]), i*2*Pi/k),
               rotxz(vector([0,2,0]), i*2*Pi/k),
               rotxz(vector([0,0,3]), i*2*Pi/k)
              ],
              axes = normal,
              style = patch,
              orientation = [-25,75]):
L := [op(L), plt[i]]:
od:
display(L, insequence=true);

k := 8: # number of frames
L := []:
for i from 0 to k do
plt[i] := polygonplot3d(
              [rotyz(vector([1,0,0]), i*2*Pi/k),
               rotyz(vector([0,2,0]), i*2*Pi/k),
               rotyz(vector([0,0,3]), i*2*Pi/k)
              ],
              axes = normal,
              style = patch,
              orientation = [-25,75]):
L := [op(L), plt[i]]:
od:
display(L, insequence=true);

# composed linear transformation
T := evalm(Txy &* Txz &* Tyz);
rot := proc(vec, angle1, angle2, angle3)
   # vec to be rotated by angle
   evalm(subs({theta = angle1, eta = angle2, xi = angle3}, evalm(T))
         &* vec);
end;

k := 12: # number of frames
L := []:
for i from 0 to k do
plt[i] := polygonplot3d(
              [rot(vector([1,0,0]), i*2*Pi/k, i*2*Pi/k, i*2*Pi/k),
               rot(vector([0,2,0]), i*2*Pi/k, i*2*Pi/k, i*2*Pi/k),
               rot(vector([0,0,3]), i*2*Pi/k, i*2*Pi/k, i*2*Pi/k)
              ],
              axes = normal,
              style = patch,
              orientation = [-25,75]):
L := [op(L), plt[i]]:
od:
display(L, insequence=true);