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 "rigfit".
#
# Fields of structure lcinp:
#   nodd        External nodal point identifier of dependent node
#   dofa{:, 2}  Cell array defining autonomous degrees of freedom:
#               {:, 1} Array with external nodal point identifiers
#               {:, 2} Array with degree of freedom identifiers
#
# 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 is obtained from a QR-decomposition of the
# constraint matrix of a rigid body (see rigbdy).
#
# -------------------------------------------------------------------

  signs = [-1, 1]; cols = [2, 1];   % to compute cross product

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

# Initialize

  mxdofpnt = 3;
  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.rigfit(%3.0d): ", nc);
            fprintf(msg, "unknown field %s\n", fields{n});
         end
     end
  end

# Check definition of autonomous degrees of freedom

  ndofa = 0;

  if (! isfield(lcinp, "dofa") || isempty(lcinp.dofa))
     rc = 1;
     fprintf(msg, "*E* constraints.rigfit(%3d): ", nc);
     fprintf(msg, "dofa undefined\n");
  else
     [nrowa, ncol] = size(lcinp.dofa); 
     if (ncol != 2 || ndims(lcinp.dofa) != 2)
        rc = 1;
        fprintf(msg, "*E* constraints.rigfit(%3d): ", nc);
        fprintf(msg, "bad definition of dofa\n");
     else
        nodea = cell2mat(lcinp.dofa(:, 1));
        ixa   = lookup(nodes.ids, nodea, "m");
        bad   = find(ixa == 0);
        nbad  = length(bad);
        if (nbad)
           rc = 1;
           for n = 1 : nbad
               fprintf(msg, "*E* constraints.rigfit(%3d).dofa: ", nc);
               fprintf(msg, "node %5.0d does not exist\n",
                       nodea(bad(n)));
           end
        end
        dofs = lcinp.dofa(:, 2);
        for n = 1 : nrowa
            ln(n)   = length(lcinp.dofa{n, 1});
            ndofa  += length(dofs{n}) * ln(n);
            dofs{n} = sort(dofs{n});
            if (dofs{n}(end) > 3 || dofs{n}(1) < 1)
               rc = 1;
               fprintf(msg, "*E* constraints.rigfit(%3d).dofa{%2d, 2}: ",
                       nc, n);
               fprintf(msg, "illegal dof identifier found\n");
            end
        end
     end
  end

  if (ndofa > 0 && ndofa < 3)
     rc = 1;
     fprintf(msg, "*E* constraints.rigfit(%3.0d): ", nc);
     fprintf(msg, "unstable support of rigid body\n");
  end

# Check identifier of dependent node

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

  if (rc) return; end

# Build degree of freedom arrays

  dofa = zeros(2, ndofa);
  dofd = zeros(2, 3); 
  dofd(1, :) = ixd;
  dofd(2, :) = 1 : 3;

# Build inverse constraint matrix

  coord = nodes.coor(ixd, :);
  Cinv  = sparse(ndofa, 3);
  ixca  = mat2cell(ixa(:), ln);

  j = 0;
  for n = 1 : nrowa
      ixnod = ixca{n};
      r     = nodes.coor(ixnod, :) - coord;
      for k = 1 : length(ixnod)
          for jdofa = dofs{n}
              dofa(1, ++j)   = ixnod(k);
              dofa(2, j)     = jdofa;
              Cinv(j, jdofa) = 1;
              if (jdofa < 3)
                 Cinv(j, 3) = signs(jdofa) * r(k, cols(jdofa));
              end
          end
      end
  end

# Compute constraint matrix

  [Q, R] = qr(Cinv);

  d = abs(diag(R));
  if (min(d) < eps * max(d))
     rc = 1;
     fprintf(msg, "*E* constraints.rigfit(%3.0d): ", nc);
     fprintf(msg, "unstable support of rigid body\n");
     return;
  end

  C = R \ Q';
  
# Build the constraint structure

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

end
