function [angle1, angle2] = AxesToAngles(L, axis1, axis2, axis3)
  % -----------------------------------------------------------------------
  % AxesToAngles - Calculate the angles by which a line must be rotated
  %                about axis1 and axis2 to align it with axis3.
  %
  % Args: L - a matrix of points on the line to be rotated
  %       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 align the line with
  %
  % Rtns: angle1 - the angle of rotation about axis1
  %       angle2 - the angle of rotation about axis2
  % -----------------------------------------------------------------------

  d = L(1:3, 2) - L(1:3, 1);    % calculate deltax, deltay and deltaz
  x = d(1); y = d(2); z = d(3); % extract components of the vector

  switch axis1
  case 'x' % rotate about x first

    h = sqrt(y*y + z*z);
    switch axis2
    case 'y'
      angle1 = atan2(y, z); % by (90-a) about x => zx plane
      switch axis3
      case 'x'
        angle2 = atan2(h, x); % by (b) about y => x axis
      case 'z'
        angle2 = atan2(-x, h); % by -(90-b) about y => z axis
      otherwise
        error('axis3 = %s is invalid when [axis1 axis2] = %s%s', ...
          axis3, axis1, axis2 ...
        );
      end
    case 'z'
      angle1 = atan2(-z, y); % by (-a) about x => xy plane
      switch axis3
      case 'y'
        angle2 = atan2(x, h); % by (b) about z => y axis
      case 'x'
        angle2 = atan2(-h, x); % by (90-b) about z => x axis
      otherwise
        error('axis3 = %s is invalid when [axis1 axis2] = %s%s', ...
          axis3, axis1, axis2 ...
        );
      end
    otherwise
      error('axis2 = %s is invalid when axis1 = %s', axis2, axis1);
    end

  case 'y' % rotate about y first

    h = sqrt(x*x + z*z);
    switch axis2
    case 'z'
      angle1 = atan2(z, x); % by (a) about y => xy plane
      switch axis3
      case 'y'
        angle2 = atan2(h, y); % by (-b) about z => y axis
      case 'x'
        angle2 = atan2(-y, h); % by (90-b) about z => x axis
      otherwise
        error('axis3 = %s is invalid when [axis1 axis2] = %s%s', ...
          axis3, axis1, axis2 ...
        );
      end
    case 'x'
      angle1 = atan2(-x, z); % by -(90-a) about y => yz plane
      switch axis3
      case 'z'
        angle2 = atan2(y, h); % by (90-b) about x => z axis
      case 'y'
        angle2 = atan2(-h, y); % by (-b) about x => y axis
      otherwise
        error('axis3 = %s is invalid when [axis1 axis2] = %s%s', ...
          axis3, axis1, axis2 ...
        );
      end
    otherwise
      error('axis2 = %s is invalid when axis1 = %s', axis2, axis1);
    end

  case 'z' % rotate about z first

    h = sqrt(x*x + y*y);
    switch axis2
    case 'x'
      angle1 =atan2(x, y); % by (a) about z => yz plane
      switch axis3
      case 'z'
        angle2 = atan2(h, z); % by (90-b) about x => z axis
      case 'y'
        angle2 = atan2(-z, h); % by (-b) about x => y axis
      otherwise
        error('axis3 = %s is invalid when [axis1 axis2] = %s%s', ...
          axis3, axis1, axis2 ...
        );
      end
    case 'y'
      angle1 = atan2(-y, x); % by (90-a) about z => zx plane
      switch axis3
      case 'x'
        angle2 = atan2(z, h); % by (b) about y => x axis
      case 'z'
        angle2 = atan2(-h, z); % by -(90-b) about y => z axis
      otherwise
        error('axis3 = %s is invalid when [axis1 axis2] = %s%s', ...
          axis3, axis1, axis2 ...
        );
      end
    otherwise
      error('axis2 = %s is invalid when axis1 = %s', axis2, axis1);
    end
  otherwise
    error('axis1 = %s is invalid', axis1);
  end
end
