function RotMovieBody(L, axis1, axis2, axis3)
  % -----------------------------------------------------------------------
  % RotMovieBody - Translate the line to the origin, scale it to unit
  %                length, and rotate it about axis1 and then about axis2
  %                so as to align it with axis3.
  %
  % Args: L - a matrix of endpoints of the line to be operated upon
  %       axis1 - 'x', 'y' or 'z': axis of the first rotation
  %       axis2 - 'x', 'y' or 'z': axis of the second rotation
  %       axis3 - 'x', 'y' or 'z': axis to which L will be aligned
  %
  % Rtns: none
  % -----------------------------------------------------------------------

  % video parameters
  FileName = ['Rot3D' axis1 axis2 axis3 '.avi'];
  FrameRate = 10; % frames per second
  KeyRate = 10; % key-frames per sec (=FrameRate => no interpolation)

  RotationRate = 0.04; % rotations per second
  rate = 2 * pi * RotationRate / FrameRate; % radians per frame

  % create the avi file and set its parameters
%  video = avifile( ...
%    FileName,'fps', FrameRate, 'keyframe', KeyRate, ...
%    'quality', 100, 'compression', 'None' ...
%  );

  % define colors
  Clr1 = [.0 .0 .0]; % the line being operated on
  Clr2 = [.7 .4 .4]; % the arc traced by the endpoint in rotation 1
  Clr3 = [.4 .7 .7]; % the arc traced by the endpoint in rotation 2
  Clr4 = [.7 .7 .4]; % construction lines
  Clr5 = [.95 .90 .90]; % srface swepth out by the line in rotation 1
  Clr6 = [.90 .95 .95]; % srface swepth out by the line in rotation 2

  % calculate matrices required to normalize the line
  T = [ % translate the starting point of the line to the origin
    1 0 0 -L(1,1)
    0 1 0 -L(2,1)
    0 0 1 -L(3,1)
    0 0 0    1
  ];
  u = L(1:3,2) - L(1:3,1); % delta x, delta y and delta z of the line
  h = 1 / sqrt(dot(u, u)); % 1 / (length of line)
  S = [ % scale the line to unit length
    h 0 0 0
    0 h 0 0
    0 0 h 0
    0 0 0 1
  ];

  % calculate the rotation matrices
  [angle1, angle2] = AxesToAngles(L, axis1, axis2, axis3);
  frames1 = abs(round(angle1 / rate));
  frames2 = abs(round(angle2 / rate));
  R1 = ConstructRotation(axis1, angle1 / frames1);
  R2 = ConstructRotation(axis2, angle2 / frames2);

  % calculate matrix to project onto the plane normal to axis1
  plane = NormalPlane(axis1);
  Pn = ConstructProjection2(plane);

  % calculate matrices to project onto axes of the plane normal to axis1
  P2 = ConstructProjection1(plane(1));
  P3 = ConstructProjection1(plane(2));

  L = S * T * L; % put the line in its original position
  Lp = Pn * L; % project onto the normal plane (normal to axis1)
  Lax1 = [P2*L(:,2) Lp(:,2) P3*L(:,2)]; % project onto axes of normal plane
  Lpn = [Lp L(:,2)]; % the projecton and a perpendicular down to it
  Lo1 = L; Lpo = Lpn; % save initial positions

  % setup and draw the initial frame
  RotMovieFrameInit(axis1, axis2, axis3, angle1, angle2);
  PlotFig3(L, [1 0], Clr1); % the line
  PlotFig3(Lpn, [1 1 0], Clr4); % projection with perpendicular
  PlotFig3(Lax1, [1 1 0], Clr4); % projection to axes of normal plane

%  hold off; video = addframe(video, getframe);
  hold off; getframe;

  % perform the rotation about axis 1
  Arc1 = L(:,2); Op1 = []; Arc2 = Lpn(:,2); Op2 = []; % initialize the arcs
  X1 = []; Y1 = []; Z1 = []; C1 = [];
  for i = 1:frames1
    V1 = L; % get frst two patch vertices
    % rotate incrementally, extending the arcs as we go
    L = R1 * L; Lpn = R1 * Lpn;
    Arc1 = [Arc1 L(:,2)]; Op1 = [Op1 1];
    Arc2 = [Arc2 Lpn(:,2)]; Op2 = [Op2 1];
    V1 = [V1 L(:,2)]; % get last patch vertex
    X1 = [X1 V1(1,:)']; Y1 = [Y1 V1(2,:)']; Z1 = [Z1 V1(3,:)'];
    C1(1,i,1:3) = Clr5;

    % setup the frame and draw the constant parts
    RotMovieFrameInit(axis1, axis2, axis3, angle1, angle2);
    PlotFig3(Lo1, [1 0], Clr1);
    PlotFig3(Lpo, [1 1 0], Clr4);
    PlotFig3(Lax1, [1 1 0], Clr4);

    % draw the evolving parts
    PlotFig3(L, [1 0], Clr1);
    PlotFig3(Lpn, [1 1 0], Clr4);
    PlotFig3(Arc1, Op1, Clr2);
    PlotFig3(Arc2, Op2, Clr4);

    set(patch(X1, Y1, Z1, C1), ... surface swept out by L
      'EdgeLighting','phong', 'EdgeColor', 'none', ...
      'FaceLighting','phong', 'FaceAlpha', .3 ...
    );

%    hold off; video = addframe(video, getframe);
    hold off; getframe;
  end

  % calculate matrix to project onto the plane normal to axis2
  plane = NormalPlane(axis2);
  P2 = ConstructProjection1(plane(1));
  P3 = ConstructProjection1(plane(2));

  Lax2 = [P2*L(:,2) L(:,2) P3*L(:,2)]; % project onto axes of normal plane
  Lo2 = L; % save the postiion of the line
 
  % setup the frame and draw the constant parts
  RotMovieFrameInit(axis1, axis2, axis3, angle1, angle2);
  PlotFig3(Lo1, [1 0], Clr1);
  PlotFig3(Lpo, [1 1 0], Clr4);
  PlotFig3(Lax1, [1 1 0], Clr4);
  PlotFig3(Arc1, Op1, Clr2);
  PlotFig3(Arc2, Op2, Clr4);

  set(patch(X1, Y1, Z1, C1), ... surface swept out by L
    'EdgeLighting','phong', 'EdgeColor', 'none', ...
    'FaceLighting','phong', 'FaceAlpha', .3 ...
  );

  % draw the new parts
  PlotFig3(L, [1 0], Clr1);
  PlotFig3(Lax2, [1 1 0], Clr4);

%  hold off; video = addframe(video, getframe);
  hold off; getframe;

  % perform the rotation about axis 2
  Arc3 = L(:,2); Op3 = []; % initialize the arc
  for i = 1:frames2
    V1 = L; % get first two patch vertices
    % rotate incrementally, extending the arc as we go
    L = R2 * L; Arc3 = [Arc3 L(:,2)]; Op3 = [Op3 1];
    V1 = [V1 L(:,2)]; % get last patch vertex
    X1 = [X1 V1(1,:)']; Y1 = [Y1 V1(2,:)']; Z1 = [Z1 V1(3,:)'];
    C1(1,i+frames1,1:3) = Clr6;

    % setup the frame and draw the constant parts
    RotMovieFrameInit(axis1, axis2, axis3, angle1, angle2);
    PlotFig3(Lo1, [1 0], Clr1);
    PlotFig3(Lpo, [1 1 0], Clr4);
    PlotFig3(Lax1, [1 1 0], Clr4);
    PlotFig3(Arc1, Op1, Clr2);
    PlotFig3(Arc2, Op2, Clr4);
    PlotFig3(Lo2, [1 0], Clr1);
    PlotFig3(Lax2, [1 1 0], Clr4);

    % draw the evolving parts
    PlotFig3(L, [1 0], Clr1);
    PlotFig3(Arc3, Op3, Clr3);

    set(patch(X1, Y1, Z1, C1), ... surface swept out by L
      'EdgeLighting','phong', 'EdgeColor', 'none', ...
      'FaceLighting','phong', 'FaceAlpha', .25 ...
    );

%    hold off; video = addframe(video, getframe);
    hold off; getframe;
  end
%  video = close(video);
end
