function [resp, ES, rc] = mfs_mfreq(cmp, f, nb, lc, meth)

# usage  [resp, ES, rc] = mfs_mfreq(cmp, f, nb, lc, meth)
#
# Input  cmp     Structure with component
#        f(:)    List of excitation frequencies
#        nb      Number of excitation frequencies per halfpower bandwidth
#        lc      Loadcase
#        meth    Method: 2 = modal reduction
#                        3 = enhanced modal reduction
#
# Output resp    Structure with frequency response
#        ES(:)   Strain energies
#        rc      Return code: 0 means no errors
#
# ------------------------------------------------------------------------

# Initialize

  resp = struct();
  rc   = 1;
  ES   = [];

# Check availability of data needed

  if (! isfield(cmp, "modes"))
     printf("*E* mfs_freqresp: Normal modes not found\n");
     return;
  end
  if (! isfield(cmp, "load"))
     printf("*E* mfs_freqresp: Loads not found\n");
     return;
  end
  if (lc > cmp.load.nofldc || lc < 1)
     printf("*E* mfs_freqresp: Loadcase %2d does not exist\n", lc);
     return;
  end
  if (! isfield(cmp, "damping"))
     printf("*E* mfs_freqresp: No damping defined\n");
     return;
  end

# Determine load type

  inflc = cmp.load.inflc(lc);
  switch inflc
  case 1
     load = 1;
  case 2
     load = 2;
  case 4
     load = 3;
  case 8
     load = 4;
  otherwise
     printf("mfs_mfreq: load types cannot be mixed in one loadcase\n");
     return;
  end

  rc = 0;

  force  = (load == 1);
  enfmot = (load > 1);

# Get information on degrees of freedom

  ndofg = cmp.dofs.ndofg;
  ndofl = cmp.dofs.ndofl;
  ndofd = cmp.dofs.ndofd;
  dofl  = cmp.dofs.dofl;
  if (enfmot)
     ndofp = cmp.dofs.ndofp;
     dofp  = cmp.dofs.dofp;
  end
  if (ndofd)
     dofd = cmp.dofs.dofd;
  end
  ndofr = cmp.stiff.ndofr;
  rbs   = (ndofr && (meth == 3 || enfmot));
  if (rbs)
     dofe = cmp.stiff.dofe;
  end

# Get modal data

  fr = cmp.modes.freq;
  nr = length(fr);
  w  = cmp.modes.omega;
  X  = cmp.modes.disp;

# Partition stiffness and mass matricess

  flags = [1, 0];
  if (enfmot)
     flags = [1, 1];
  end
  [Kll, Klp] = mfs_matpart(cmp.stiff.K, cmp.dofs, flags);
  [Mll, Mlp] = mfs_matpart(cmp.mass.M, cmp.dofs, flags);
  if (rbs)
     Kee = Kll(dofe, dofe);
     Xl  = X(dofl, :);
     Xr  = Xl(:, 1 : ndofr);
     XTM = Xr' * Mll;
  end

# Get Base Motion

  if (enfmot)
     Ub       = zeros(ndofg, 1, "double");
     if (load == 2)
        Ub(dofp) = cmp.load.u(dofp, lc);
     elseif (load == 3)
        Ub(dofp) = cmp.load.v(dofp, lc);
     else
        Ub(dofp) = cmp.load.a(dofp, lc);
     end
     if (ndofr)
        Ubl = zeros(ndofl, 1, "double");
        Ubl(dofe) = - Kee \ Klp(dofe, :) * Ub(dofp);
        Ubl -= Xr * (XTM * Ubl);
        Ub(dofl) = Ubl;
     else
        Ub(dofl) = - Kll \ Klp * Ub(dofp);
     end
     if (ndofd)
        Ub(dofd) = cmp.dofs.C(dofd, :) * Ub;
     end
     Mqb = X' * cmp.mass.M * Ub;
  end

# Damping

  switch cmp.damping.type
  case "ratios"
     nd = length(cmp.damping.data);
     if (nd >= nr)
        D = cmp.damping.data(1 : nr)';
     else
        Dl = cmp.damping.data(nd);
        D  = [cmp.damping.data'; Dl( ones(nr - nd, 1))];
     end
     if (enfmot)
        Dqb = 2 * diag(w .* D) * X(dofl, :)' * Mll * Ub(dofl);
     end
  case "Rayleigh"
     if (ndofr)
        elmod = ndofr + 1 : cmp.modes.nofmod;
        D         = zeros(cmp.modes.nofmod, 1);
        D(elmod)  = 0.5 * cmp.damping.data(1) * w(elmod);
        D(elmod) += 0.5 * cmp.damping.data(2) ./ w(elmod);
     else
        D = 0.5 * (cmp.damping.data(1) * w + cmp.damping.data(2) ./ w);
     end
     if (enfmot)
        Kqg = X' * cmp.stiff.K;
        Mqg = X' * cmp.mass.M;
        Dqb = cmp.damping.data(1) * Kqg * Ub;
        Dqb = Dqb + cmp.damping.data(2) * Mqg * Ub;
     end
  end

# Excitation frequencies in halfpower bandwith

  if (nb > 0)
     [freq, nfreq] = mfs_addfreq(f, fr, D, nb);
  else
     freq  = f;
     nfreq = length(f);
  end

# Get excitation

  if (force)
     rhs = X' * cmp.load.f(:, lc);
  else
     we = 2 * pi * freq;
     if (load == 2)
        rhs = Mqb * we.^2 - i * Dqb * we;
     elseif (load == 3)
        rhs = -i * Mqb * we - Dqb;
     else
        rhs = - Mqb + i * Dqb * (1 ./ we);
     end
  end

# Get response of modal coefficients

  Ginv = (2 * pi)^2 * (fr.^2 - freq.^2 + 2 * i * (D .* fr) * freq);
  Q    = rhs ./ Ginv;

  ES = zeros(1, nfreq);
  for k = 1 : nfreq
      E     = 0.25 * Q(:, k)' * cmp.modes.Kmodal * Q(:, k);
      ES(k) = real(E);
  end

# Static response of neglected modes

  if (meth == 3)
     if (force)
        rhs = mfs_matpartr(cmp.load.f(:, lc), cmp.dofs);
     else
        MUb = cmp.mass.M * Ub;
        rhs = -mfs_matpartr(MUb, cmp.dofs);
     end
     R  = rhs - Mll * X(dofl, :) * (X(dofl, :)' * rhs);
     Uc = zeros(ndofg, 1);
     if (rbs)
        Ucl = zeros(ndofl, 1);
        Ucl(dofe) = Kee \ R(dofe);
        Ucl -= Xr * (XTM * Ucl);
        Uc(dofl) = Ucl;
     else
        Uc(dofl) = Kll \ R;
     end
     if (ndofd)
        Uc(dofd) = cmp.dofs.C(dofd, :) * Uc;
     end
     E  = 0.25 * Uc' * cmp.stiff.K * Uc;
     if (load == 1 || load == 4)
        ES = ES + E;
     elseif (load == 2)
        ES = ES + (2 * pi * freq).^4 * E;
     elseif (load == 3)
        ES = ES + (2 * pi * freq).^2 * E;
     end
  end

# Contribution of base motion to strain energy

  if (enfmot)
     E = 0.25 * Ub' * cmp.stiff.K * Ub;
     if (load == 2)
        ES = ES + E;
     elseif (load == 3)
        ES = ES + E * (2 * pi * freq).^-2;
     elseif (load == 4)
        ES = ES + E * (2 * pi * freq).^-4;
     end
  end

# Build frequency response structure

  resp = struct("method", meth, "type", load, "nfreq", nfreq,
                "freq", freq, "es", ES, "Q", Q, "nback", 0);
  if (meth == 3)
     resp = setfield(resp, "mct", Uc);
  end
  if (enfmot)
     resp = setfield(resp, "base", Ub);
  end

end
