function R = mfs_rigidmotion(nodes, mxdofpnt, refpnt)

# usage: R = mfs_rigidmotion(nodes, mxdofpnt, refpnt)
#
# Input  nodes       Structure with nodal point data
#        mxdofpnt    Maximum number of dofs per node
#        refpnt(:)   Coordinates of a reference point
#
# Output R(:, :)   Rigid motion matrix:
#                  2d: R(:, 3): tx, ty, rz
#                  3d: R(:, 6): tx, ty, tz, rx, ry, rz
#
# The function generates a set of unit rigid body motions from the
# coordinates.
#
# ------------------------------------------------------------------------

# Check arguments

  if (nargin != 3 || nargout != 1)
     print_usage();
  end

# Number of coordinates

  ncoor = nodes.ncoor;
  nrot  = mxdofpnt - ncoor;

# Initiallize rigid body modes

  ndofg = nodes.nofnod * mxdofpnt;

  if (ncoor == 2)
     nrig = 3;
  else
     nrig = 6;
  end

  R = zeros(ndofg, nrig);

# Rigid body translations

  k = 0;
  for n = 1 : nodes.nofnod
      for l = 1 : ncoor
          R(++k, l) = 1;
      end
      k += nrot;
  end

# Rigid body rotations

  k = 0;
  if (ncoor == 2)        % 2 dimensions
     for n = 1 : nodes.nofnod
         R(++k, 3) = -nodes.coor(n, 2) + refpnt(2);
         R(++k, 3) =  nodes.coor(n, 1) - refpnt(1);
         if (nrot)
            R(++k, 3) = 1;
         end
     end
  else                   % 3 dimensions
     for n = 1 : nodes.nofnod
         x = nodes.coor(n, :) - refpnt;
         R(++k, 5) =  x(3);
         R(  k, 6) = -x(2);
         R(++k, 4) = -x(3);
         R(  k, 6) =  x(1);
         R(++k, 4) =  x(2);
         R(  k, 5) = -x(1);
         if (nrot)
            for l = 1 : 3
                R(++k, 3 + l) = 1;
            end
         end
     end
  end

end
