function [resp, ES, WR, rc] = mfs_dfreq(cmp, f, nb, lc)

# usage  [resp, ES, WR, rc] = mfs_dfreq(cmp, f, nb, lc)
#
# Input  cmp     Structure with component
#        f(:)    List of excitation frequencies
#        nb      Number of excitation frequencies per halfpower bandwidth
#        lc      Loadcase
#
# Output resp    Structure with frequency response
#        ES(:)   Strain energies
#        WR(:)   Work of residual
#        rc      Return code: 0 means no errors
#
# ------------------------------------------------------------------------

# Initialize

  resp = struct();
  ES   = 0;
  WR   = 0;
  rc   = 1;

# Check availability of data needed

  if (! isfield(cmp, "load"))
     printf("mfs_freqresp: Loads not found\n");
     return;
  end
  if (lc > cmp.load.nofldc || lc < 1)
     printf("mfs_freqresp: Loadcase %2d does not exist\n", lc);
     return;
  end
  if (nb > 0 && ! isfield(cmp, "modes"))
     warning("mfs_freqresp: Normal modes not found\n");
     nb = -1;
  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_dfreq: 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

# Partition stiffness and mass matrices

  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);

# Damping matrix

  dampdef = 0;

  if (isfield(cmp, "damping"))
     switch cmp.damping.type
     case "Rayleigh"
        Dll = cmp.damping.data(1) * Kll + cmp.damping.data(2) * Mll;
        if (enfmot)
           Dlp = cmp.damping.data(1) * Klp + cmp.damping.data(2) * Mlp;
        end
        dampdef = 1;
     otherwise
        warning("mfs_freqresp: Damping type %s ignored in direct response\n", ...
               cmp.damping.type);
     end
  end

  if (! dampdef) 
     warning("mfs_freqresp: No damping defined\n");
     Dll = sparse(ndofl, ndofl);
     if (enfmot)
        Dlp = sparse(ndofl, ndofp);
     end
  end

# Excitation frequencies in halfpower bandwith

  if (nb > 0)
     D = cmp.damping.data(1) * cmp.modes.omega;
     D = 0.5 * (D + cmp.damping.data(2) ./ cmp.modes.omega);
     [freq, nfreq] = mfs_addfreq(f, cmp.modes.freq, D, nb);
  else
     freq  = f;
     nfreq = length(freq);
  end

# Excitation

  if (force)
     rhs = mfs_matpartr(cmp.load.f(:, lc), cmp.dofs);
  else
     if (load == 2)
        Mup = Mlp * cmp.load.u(dofp, lc);
        Dup = Dlp * cmp.load.u(dofp, lc);
        Kup = Klp * cmp.load.u(dofp, lc);
     elseif (load == 3)
        Mvp = Mlp * cmp.load.v(dofp, lc);
        Dvp = Dlp * cmp.load.v(dofp, lc);
        Kvp = Klp * cmp.load.v(dofp, lc);
     else
        Map = Mlp * cmp.load.a(dofp, lc);
        Dap = Dlp * cmp.load.a(dofp, lc);
        Kap = Klp * cmp.load.a(dofp, lc);
     end
  end

# Response

  U  = zeros(ndofg, nfreq);
  ES = zeros(1, nfreq);
  WR = zeros(1, nfreq);

  for k = 1 : nfreq
      w  = 2 * i * pi * freq(k);
      wq = w^2;
      Kdyn = Kll + w * Dll + wq * Mll;
      if (enfmot)
         if (load == 2)
            rhs = -(wq * Mup + w * Dup + Kup);
            U(dofp, k) = cmp.load.u(dofp, lc);
         elseif (load == 3)
            rhs = -(w * Mvp + Dvp + Kvp / w);
            U(dofp, k) = cmp.load.v(dofp, lc) / w;
         else
            rhs = -(Map + Dap / w + Kap / wq);
            U(dofp, k) = cmp.load.a(dofp, lc) / wq;
         end
      end
      U(dofl, k) = Kdyn \ rhs;
      if (ndofd)
         U(dofd, k) = cmp.dofs.C(dofd, :) * U(:, k);
      end
      E          = U(:, k)' * cmp.stiff.K * U(:, k);
      ES(k)      = 0.25 * real(E);
      R          = Kdyn * U(dofl, k) - rhs;
      WR(k)      = 0.25 * real(U(dofl, k)' * R);
  end

# Build frequency response structure

  resp = struct("method", 1, "type", load, "nfreq", nfreq,
                "freq", freq, "es", ES, "wr", WR, 
                "nback", nfreq, "freqback", freq, "U", U);

end
