function [cs, mxdofpnt, rc] = mfs_cs(msg, nc, lcinp, nodes)

# usage: [cs, mxdofpnt, rc] = mfs_cs(msg, nc, lcinp, nodes)
#
# Input  msg      File handle of message file
#        nc       Number of constraint
#        lcinp    Structure with definition of linear constraint
#        nodes    Structure with nodal point data 
#
# Output cs       Structure with constraint matrix
#        mxdofpnt Maximum number of degrees of freedom per nodal
#                 point
#        rc       Return code: 0 means no errors
#                              1 means errors
#
# The function computes the constraint matrix of a linear constraint
# of type "rigbdy".
#
# Fields of structure lcinp:
#   dofs(:)     Degree of freedom identifiers
#   noda        External nodal point identifier of autonomous node
#   nodd(:)     External nodal point identifiers of dependent nodes
#
# Fields of structure cs:
#   dofa(2, :)  Autonomous degrees of freedom:
#                 internal nodal point identifier, dof identifier
#   dofd(2, :)  Dependent degrees of freedom:
#                 internal nodal point identifier, dof identifier
#    C(:, :)   Constraint matrix of one constraint
#              (rows = dep. dofs, cols = autonomous dofs)
#
# The constraint matrix reads
#
#     ux     1   0   0     0   r(3) -r(2)
#     uy     0   1   0  -r(3)   0    r(1)
#     uz     0   0   1   r(2) -r(1)   0
#     phix   0   0   0    1     0     0
#     phiy   0   0   0    0     1     0
#     phiz   0   0   0    0     0     1
#
# -------------------------------------------------------------------

  cx = [1, 2, 3, 1, 2];   % Indices to compute cross product

  legal_fields = {"dofs", "noda", "nodd"};
  nlegal = length(legal_fields);

# Initialize

  mxdofpnt = 6;
  cs       = struct("dofa", [], "dofd", [], "C", []);
  rc       = 0;

# Check for unknown fields

  fields = fieldnames(lcinp);
  defleg = ismember(fields, legal_fields);
  if (sum(defleg) != nlegal)
     rc = 1;
     for n = 1 : nlegal
         if (! defleg(n))
            fprintf(msg, "*E* constraints.rigbdy(%3.0d): ", nc);
            fprintf(msg, "unknown field %s\n", fields{n});
         end
     end
  end

# Check degree of freedom identifiers

  if (! isfield(lcinp, "dofs") || isempty(lcinp.dofs))
     rc = 1;
     fprintf(msg, "*E* constraints.rigbdy(%3d): ", nc);
     fprintf(msg, "dofs undefined\n");
  elseif ((max(lcinp.dofs)) > 6 || min(lcinp.dofs) < 1)
     fprintf(msg, "*E* constraints.rigbdy(%3.0d): ", nc);
     fprintf(msg, "illegal dof identifier found\n");
  end

# Check identifier of autonomous node

  if (! isfield(lcinp, "noda") || isempty(lcinp.noda))
     rc = 1;
     fprintf(msg, "*E* constraints.rigbdy(%3d): ", nc);
     fprintf(msg, "noda undefined\n");
  else
     if (length(lcinp.noda) > 1)
        rc = 1;
        fprintf(msg, "*E* constraints.rigbdy(%3d).noda: ", nc);
        fprintf(msg, "more than one node defined\n");
     else
        ixa = lookup(nodes.ids, lcinp.noda, "m");
        if (! ixa)
           rc = 1;
           fprintf(msg, "*E* constraints.rigbdy(%3d).noda: ", nc);
           fprintf(msg, "node %5.0d does not exist\n", lcinp.noda); 
        end
     end
  end

# Check identifiers of dependent nodes

  if (! isfield(lcinp, "nodd") || isempty(lcinp.nodd))
     rc = 1;
     fprintf(msg, "*E* constraints.rigbdy(%3d): ", nc);
     fprintf(msg, "nodd undefined\n");
  else
     nnodd = length(lcinp.nodd);
     ixd   = lookup(nodes.ids, lcinp.nodd, "m");
     bad   = find(ixd == 0);
     nbad  = length(bad);
     if (nbad)
        rc = 1;
        for n = 1 : nbad
            fprintf(msg, "*E* constraints.rigbdy(%3d).nodd(%2d): ",
                    nc, bad(n));
            fprintf(msg, "node %5.0d does not exist\n", 
                    lcinp.nodd(bad(n)));
        end
     end
  end

  if (rc) return; end

# Build structure with constraint matrix

  dofs  = sort(lcinp.dofs);
  ndofs = length(lcinp.dofs);
  ndoft = ndofs * nnodd;

  dofa = zeros(2, 6);
  dofa(1, :) = ixa;
  dofa(2, :) = 1 : 6;
  coora      = nodes.coor(ixa, :);

  dofd = zeros(2, ndoft);
  C    = sparse(ndoft, 6);

  j = 0;
  for n = 1 : nnodd
      inode = ixd(n);
      r     = nodes.coor(inode, :) - coora;
      for m = 1 : ndofs
          jdofd        = dofs(m);
          dofd(1, ++j) = inode;
          dofd(2, j)   = jdofd;
          C(j, jdofd)  = 1;
          if (jdofd < 4)
             c1 = cx(jdofd + 1); c2 = cx(jdofd + 2);
             C(j, c1 + 3) =  r(c2);
             C(j, c2 + 3) = -r(c1);
          end
      end
  end

  cs = struct("dofa", dofa, "dofd", dofd, "C", C);

end
