function mem = mfs_meffmass(fid, cmp, refpnt)

# usage: mem = mfs_meffmass(fid, cmp, refpnt);
#
# Input  fid          File handle
#        cmp          Structure with component
#        refpnt(:)    Coordinates of reference point (optional)
#
# Output mem          Structure with absolute and relative modal
#                     effective masses (optional)
#
# The function computes the modal effective masses and writes them to
# the output file.
#
# Fields of structure mem:
#   refpnt(ndim)      Coordinates of reference point
#   abs(nofmod, ncol) Absolute modal effective masses
#   rel(nofmod, ncol) Relative modal effective masses
#
#   nofmod is the number of normal modes and ncol is 3 or 6.
#
# --------------------------------------------------------------------

# Copyright (c) 2022 by Johannes Wandinger

  t0 = clock();

  mem = struct();

# Check arguments

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

  if (! isstruct(cmp))
     error("mfs_meffmass: second argument must be a structure");
  end

# Check component type

  if (! strcmp(cmp.type, "solid"))
     error("mfs_meffmass: component type must be \"solid\"\n");
  end

# Check availability of data needed

  if (! isfield(cmp, "modes"))
     printf("*W* Normal modes not available\n");
     return;
  end

  mfs_paths("add", "mfs_meffmass.m", "solid", "util");

# Get reference point

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

# Get rigid body motion

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

# Compute rigid body mass matrix

  dofl = cmp.dofs.dofl;
  Mll  = mfs_matpart(cmp.mass.M, cmp.dofs, [1, 0]);
  Rl   = R(dofl, :);
  MRl  = Mll * Rl;
  Mrr  = Rl' * MRl;
  Mr   = diag(Mrr)';

# Compute modal participation factors

  Xl = cmp.modes.disp(dofl, :);
  G  = Xl' * MRl;

# Compute modal effective masses

  meff    = G .* G;
  meffrel = meff ./ Mr;
  mrelsum = cumsum(meffrel);

  if (nargout == 1)
     mem.refpnt = refpnt;
     mem.abs    = meff;
     mem.rel    = meffrel;
  end

# Output results

  fr     = cmp.modes.freq;
  nmodes = length(fr);

  fprintf(fid, "\n");
  for k = 1 : 8
      fprintf(fid, "----------");
  end
  fprintf(fid, "\n\nModal effective masses of component ""%s""\n\n", ...
          inputname(2));
  fprintf(fid, "  Coordinates of reference point: ");

  switch cmp.subtype

  case "2d"

     fprintf(fid, " %10.4f, %10.4f\n\n", refpnt(1 : 2));
     fprintf(fid, "  mode   frequency      tx          ty          rz    ");
     fprintf(fid, "  sum tx   sum ty   sum rz\n\n"); 

     for m = 1 : nmodes
         fprintf(fid, "  %4d  %7.2f Hz  ", m, fr(m));
         fprintf(fid, "%10.4e  %10.4e  %10.4e  %7.5f  %7.5f  %7.5f\n", ...
                 meffrel(m, :), mrelsum(m, :));
     end

  case "3d"

     fprintf(fid, " %10.4f, %10.4f, %10.4f\n\n", refpnt);
     fprintf(fid, "  mode   frequency      tx          ty          tz    ");
     fprintf(fid, "  sum tx   sum ty   sum tz\n\n"); 

     for m = 1 : nmodes
         fprintf(fid, "  %4d  %7.2f Hz  ", m, fr(m));
         fprintf(fid, "%10.4e  %10.4e  %10.4e  %7.5f  %7.5f  %7.5f\n", ...
                 meffrel(m, 1 : 3), mrelsum(m, 1 : 3));
     end

     fprintf(fid, "\n");

     fprintf(fid, "  mode   frequency      rx          ry          rz    ");
     fprintf(fid, "  sum rx   sum ry   sum rz\n\n"); 

     for m = 1 : nmodes
         fprintf(fid, "  %4d  %7.2f Hz  ", m, fr(m));
         fprintf(fid, "%10.4e  %10.4e  %10.4e  %7.5f  %7.5f  %7.5f\n", ...
                 meffrel(m, 4 : 6), mrelsum(m, 4 : 6));
     end

  end

  fprintf(fid, "\n");

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

end
