function mp = mfs_massproperties(fid, cmp, refpnt)

# usage: mp = mfs_massproperties(fid, cmp, refpnt)
#
# Input  fid         File handle of output file
#        cmp         Structure with component
#        refpnt(:)   Coordinates of reference point (optional)
#
# Output mp          Structure with mass properties (optional)
#
# The function computes the mass properties and writes them to the
# output file. If requested, the mass properties are also returned
# in the structure mp which has the following fields:
#
#   mp.refpnt  Coordinates of reference point
#   mp.mrr     Rigid body mass matrix
#   mp.cm      Coordinates of the center of mass
#   mp.m       Mass
#   mp.J       Inertia tensor
#
# ------------------------------------------------------------------------

# Copyright (c) 2016 by Johannes Wandinger

  t0 = clock();

# Check arguments

  if (nargin < 2 || nargout > 1)
     print_usage();
  end

  if (! is_valid_file_id(fid))
     error("mfs_massproperties: first argument does not refer to an open file\n");
  end
  if (! isstruct(cmp))
     fclose(fid);
     error("mfs_massproperties: second argument must be a structure\n");
  end

  ncoor = cmp.nodes.ncoor;

  if (nargin < 3)
     refpnt = zeros(1, ncoor);
  end

  if (! strcmp(cmp.type, "solid"))
     printf("*W* mfs_massproperties: component type must be \"solid\"\n");
     return;
  end

  if (! isfield(cmp, "mass"))
     printf("*W* mfs_massproperties: mass matrix mssing\n");
     return;
  end

  mfs_paths("add", "mfs_massproperties.m", "solid");

# Rigid body motion

  R = mfs_rigidmotion(cmp.nodes, cmp.dofs.mxdofpnt, refpnt);

# Rigid body mass matrix

  mrr = R' * cmp.mass.M * R;

# Mass

  m   = mrr(1);

# Inertia tensor

  if (ncoor == 2)
     J = mrr(3, 3);
  else
     J = mrr(4 : 6, 4 : 6);
  end

# Center of Mass

  if (m > 0)
     minv = 1 / m;
     if (ncoor == 2)
        cm(1) = refpnt(1) + mrr(2, 3) * minv;
        cm(2) = refpnt(2) - mrr(1, 3) * minv;
     else
        cm(1) =  refpnt(1) + mrr(2, 6) * minv;
        cm(2) =  refpnt(2) + mrr(3, 4) * minv;
        cm(3) =  refpnt(3) + mrr(1, 5) * minv;
     end
  else
     cm = refpnt;
  end

# Inertia tensor with respect to center of mass

  dist = cm - refpnt;
  rsq  = dist * dist';
  if (ncoor == 2)
     JS = J - m * rsq;
  else
     JS = J - m * (eye(3, 3) * rsq - dist' * dist);
  end

# Output results

  fprintf(fid, "\nMass properties of component ""%s""\n\n", ...
          inputname(2));

  if (ncoor == 2)
     fprintf(fid, "  Coordinates of reference point: ");
     fprintf(fid, " %10.4f, %10.4f\n\n", refpnt(1 : 2));
     fprintf(fid, "  Rigid body mass matrix:\n\n");
     for k = 1 : 3
         fprintf(fid, "     %11.4e  %11.4e  %11.4e\n", 
                 mrr(k, 1), mrr(k, 2), mrr(k, 3));
     end
     fprintf(fid, "\n  Mass = %11.4e, Moment of inertia = %11.4e\n", ...
             m, J);
     fprintf(fid, "\n  Coordinates of center of mass: ");
     fprintf(fid, "%10.4f, %10.4f\n", cm);
     fprintf(fid, "\n  Moment of inertia with respect to center of mass: ");
     fprintf(fid, "%11.4e\n\n", JS);
  else
     fprintf(fid, "  Coordinates of reference point: ");
     fprintf(fid, " %10.4f, %10.4f,  %10.4f\n\n", refpnt);
     fprintf(fid, "  Rigid body mass matrix:\n\n");
     for k = 1 : 6
         fprintf(fid, "   ")
         for l = 1 : 6
             fprintf(fid, " %11.4e", mrr(k, l));
         end
         fprintf(fid, "\n");
     end
     fprintf(fid, "\n  Mass = %11.4e\n", m);
     fprintf(fid, "\n  Inertia tensor with respect to reference point:\n");
     for k = 1 : 3
         fprintf(fid, "\n   ")
         for l = 1 : 3
             fprintf(fid, " %11.4e", J(k, l));
         end
     end
     fprintf(fid, "\n\n  Coordinates of center of mass: ");
     fprintf(fid, "%10.4f, %10.4f, %10.4f\n", cm);
     fprintf(fid, "\n  Inertia tensor with respect to center of mass:\n");
     for k = 1 : 3
         fprintf(fid, "\n   ")
         for l = 1 : 3
             fprintf(fid, " %11.4e", JS(k, l));
         end
     end
     fprintf(fid, "\n");
  end

# Structure with mass properties

  if (nargout == 1)
     mp = struct("mrr", mrr, "m", m, "J", J, "cm", cm, "JS", JS);
  end

# End

  mfs_paths("remove");
  elapsed_time = etime(clock(), t0);
  printf("%10.4f seconds needed to compute mass properties ", ...
         elapsed_time);
  printf("of component %s\n", inputname(2));

  fflush(stdout);

end
