function [cmp, rc] = mfs_freevibx(cmp, nmodes, options)

# usage: [cmp, rc] = mfs_freevibx(cmp, nmodes, options)
#
# Input  cmp     Structure with component data
#        nmodes  Number of modes to compute
#        options Structure with options
#
# Output cmp     Structure with component data 
#                (displacement matrix added)
#        rc      Return code: 0 means no errors
#
# The function computes the free vibration modes of a solid component.
# The normal modes a normalized with respect to the mass matrix.
#
# Options:
#
#   rbc          Rigid body criterion: z = 10^rbc * eps
#                Default: 6
#   rdofs(:, 2)  Definition of rigid body dofs
#   disp         passed to eigs
#   maxit        passed to eigs
#   p            passed to eigs
#   tol          passed to eigs
#
# ---------------------------------------------------------------------

# Copyright (c) 2024 by Johannes Wandinger

  legopts = {"rbc", "rdofs", "disp", "maxit", "p", "tol"};

# Check arguments

  if (nargin != 3 || nargout != 2)
     print_usage();
  endif

# Initialize

  rc = 0;

# Check existence of stiffness and mass matrix

  if (! isfield(cmp, "stiff"))
     printf("*E* mfs_freevib: stiffness matrix missing\n");
     rc = 1;
  endif
  if (! isfield(cmp, "mass"))
     printf("*E* mfs_freevib: mass matrix missing\n");
     rc = 1;
  endif

# Process options

  defopts = fieldnames(options);
  nopts   = length(defopts);
  lopts   = ismember(defopts, legopts);

  if (sum(lopts) != nopts)
     for n = 1 : nopts
         if (! lopts(n))
            printf("*E* freevib: unknown option \"%s\"\n",
                   defopts{n});
         endif
     endfor
     rc = 1;
  endif

  if (rc) return; endif

  rc = 1;

  rbc = options.rbc;

  if (isfield(options, "rdofs"))
     rdofs = options.rdofs;
     rddef = 1;
  else
     rddef = 0;
  endif

# Get information on degrees of freedom

  ndofg = cmp.dofs.ndofg;
  ndofl = cmp.dofs.ndofl;
  dofl  = cmp.dofs.dofl;
  if ((ndofd = cmp.dofs.ndofd))
     dofd = cmp.dofs.dofd;
  endif

# Extract l-dof stiffness and mass matrix

  Kll = mfs_matpart(cmp.stiff.K, cmp.dofs, [1, 0]);
  Mll = mfs_matpart(cmp.mass.M, cmp.dofs, [1, 0]);

  if (! nnz(Mll))
     printf("E* mfs_freevib: mass matrix is empty\n");
     return;
  endif

# Look for rigid body modes

  if (rddef)

     [ixdofr, rcc] = mfs_finddof(cmp.nodes.ids, cmp.dofs.mxdofpnt,
                                 rdofs(:, 1), rdofs(:, 2));
     if (rcc)
        printf("*E* mfs_freevib: Bad definition of rigid degrees of freedom\n");
        return;
     endif

     dofr = lookup(dofl, ixdofr);
     if (! isempty(find(dofr == 0)))
        printf("*E* mfs_freevib: Bad definition of rigid degrees of freedom\n");
        return;
     endif
     dofe  = lookup(dofl, setdiff(dofl, ixdofr));
     ndofr = length(dofr);
     ndofe = ndofl - ndofr;

     n  = ndofe;
     nm = nmodes - ndofr;

  else

     [L, U, p, q, R] = lu(Kll, "vector");

     UD(p) = abs(diag(U));
     clear L; clear U;

     z     = eps * 10^rbc;
     dofe  = find(UD > z);
     ndofe = length(dofe);
     ndofr = ndofl - ndofe;

     if (ndofr)
        dofr = find(UD <= z);
        printf("           %2.0d rigid body modes detected\n", ndofr);
        n  = ndofe;
        nm = nmodes - ndofr;
     else
        n  = ndofl;
        nm = nmodes;
     endif

  endif

# Matrices for systems with rigid body modes

  if (ndofr)
     Kee = Kll(dofe, dofe);
     Ker = Kll(dofe, dofr);
     Xr  = zeros(ndofl, ndofr);
     if (nnz(Ker))    % avoid bug #52365
        Xr(dofe, :) = Kee \ Ker;
     endif
     Xr(dofr, :) = -eye(ndofr);
     Mrr = Xr' * Mll * Xr;
     [Crr, p] = chol(Mrr);
     if (p)
        printf("*E* freevib: rigid body mass matrix not positive definite\n");
        return;
     endif
     Xr  = Xr / Crr;
     Mee = Mll(dofe, dofe);
     MX  = Mll * Xr;
     MXe = MX(dofe, :);
  endif

  if (nm)

# Function for power step

     if (ndofr)
        [L, p, q] = chol(Kee, "vector", "lower");
        if (p)
           printf("*E* mfs_freevib: there may be undetected rigid body modes\n");
           printf("    To detect more rigid body modes, increase value of options.rbc.\n");
           printf("    Current value is %5.2f.\n", options.rbc);
           return; 
        endif
        Mqq = Mee(q, q);
        Mqx = MXe(q, :);
        pst = @(x) powerstep(x, L, Mqq, Mqx);
     else
        [L, p, q] = chol(Kll, "vector", "lower");
        if (p)
           printf("*E* mfs_freevib: there may be undetected rigid body modes\n");
           printf("    To detect more rigid body modes, increase value of options.rbc.\n");
           printf("    Current value is %5.2f.\n", options.rbc);
           return; 
        endif
        Mqq = Mll(q, q);
        pst = @(x) L \ (Mqq * (L' \ x));
     endif

# Solve the eigenvalue problem

     opts.issym  = 1;
     opts.isreal = 1;
     [X, lambda, flag] = eigs(pst, n, nm, "la", opts);
     if (flag)
        printf("*E* mfs_freevib: eigs failed with error %d\n", flag);
        return;
     endif
     [nrow, ncol] = size(X);
     if (ncol != nm)
        printf("*E* mfs_freevib: eigs returned %d modes instead of %d\n", ...
                ncol, nm);
        return;
     endif

# Sort results

     w2 = 1./ diag(lambda);
     [w2, index] = sort(w2);
     w2 = w2(index);
     X(q, :)  = L' \ X(:, index);

# Operations for systems with rigid body modes

     if (ndofr)
        Xe = zeros(ndofl, nm);
        Xe(dofe, :) = X;
        Xe = Xe - Xr * (MXe' * X);
        X = [Xr, Xe];
        w2r = diag(Xr' * Kll * Xr);
        w2  = [w2r; w2];
     endif

  else

     X   = Xr;
     w2r = diag(Xr' * Kll * Xr);
     w2  = w2r;

  endif

# Scale eigenvectors

  s = diag(sqrt(diag(X' * Mll * X)));
  X = X / s; 

# Circular frequencies

  omega = real(sqrt(w2));

# Frequencies and g-dof modes

  f = omega / (2 * pi);
  if (ndofl < ndofg)
     disp = zeros(cmp.dofs.ndofg, nmodes, "double");
     disp(dofl, :) = X;
     if (ndofd)
        disp(dofd, :) = cmp.dofs.C(dofd, :) * disp;
     endif
  else
     disp = X;
  endif

# Reaction loads and internal constraint loads

  if ((ndofp = cmp.dofs.ndofp) || ndofd)

     W2 =  spdiags(w2, 0, nmodes, nmodes);
     FR =  cmp.stiff.K * disp;
     FR -= cmp.mass.M * (disp * W2);

     if (ndofp)
        dofp = cmp.dofs.dofp;
        if (ndofd)
           Fp = cmp.dofs.C(:, dofp)' * FR;
           FR(dofp, :) -= Fp;
        else
           Fp = FR(dofp, :);
        endif
        reac = mfs_mat2nodedata(Fp, dofp, cmp.dofs.mxdofpnt, 4);
        reac = mfs_resnodedata(reac, cmp.nodes);
     endif

     if (ndofd)
        dofda = unique([cmp.dofs.dofa; dofd]);
        Fda   = FR(dofda, :);
        icsl  = mfs_mat2nodedata(Fda, dofda, cmp.dofs.mxdofpnt, 4);
        icsl  = mfs_resnodedata(icsl, cmp.nodes);
     endif

  endif
     
# End

  if (ndofr)
     w2(1 : ndofr) = 0;
  endif
  Kmodal = spdiags(w2, 0, nmodes, nmodes);

  modes = struct("nofmod", nmodes, "Kmodal", Kmodal, 
                 "omega", omega, "freq", f, "disp", disp);
  if (ndofp)
     modes.reac = reac;
  endif

  if (ndofd)
     modes.icsl = icsl;
  endif

  cmp.modes = modes;

  cmp.stiff.ndofr = ndofr;
  cmp.stiff.ndofe = ndofe;

  if (ndofr)
     cmp.stiff.dofr  = dofr;
     cmp.stiff.dofe  = dofe;
  endif

  rc = 0;

endfunction

function y = powerstep(x, L, Mqq, Mqx) 

# This function performs the power step in case K is not positive
# definite.
#
# ---------------------------------------------------------------------

  x = L' \ x;
  y = L \ (Mqq * x - Mqx * (Mqx' * x));

endfunction
