function [resps, respa, nfreq, rc] = mfs_mfreq(cmp, f, nb, lc, meth, opts)

# usage: [resps, respa, nfreq, rc] = mfs_mfreq(cmp, f, nb, lc, meth, opts)
#
# Input  cmp      Structure with aeroelastic component
#        f(:)     List with excitation frequencies
#        nb       Number of excitation frequencies per halfpower bandwidth
#        lc       Load case number
#        meth     Method: 2 = modal reduction
#                         4 = force summation
#
# Output resps   Structure with solid frequency response
#        respa   Structure with aerodynamic frequency response
#        nfreq   Number of excitation frequencies
#        rc      Return code: 0 means no errors
#
# ------------------------------------------------------------------------

# Copyright(c) 2024 by Johannes Wandinger

# Check arguments

  if (nargin != 6 || nargout != 4)
     print_usage();
  endif

# Initialize

  resps = struct(); respa = struct(); nfreq = 0;
  fsm   = meth == 4;
  rc    = 1;

# Check availability of data needed

  if (! isfield(cmp.solid, "modes"))
     printf("*E* mfs_freqresp: normal modes not found\n");
     return;
  endif
  if (! isfield(cmp, "load"))
     printf("*E* mfs_freqresp: loads not found\n");
     return;
  endif
  if (lc > cmp.load.nofldc || lc < 1)
     printf("*E* mfs_freqresp: loadcase %2d does not exist\n", lc);
     return;
  endif
  if (! cmp.aero.cref)
     printf("*E* mfs_freqresp: reference chord length undefined\n");
     return;
  endif
  if (! isfield(cmp.solid, "damping"))
     printf("*E* mfs_freqresp: no damping defined\n");
     return;
  endif

# Check load data

  gust = bitget(cmp.load.inflc(lc), 1);
  mnvr = bitget(cmp.load.inflc(lc), 2);

  if (! (gust | mnvr))
     printf("*E* mfs_freqresp: no loads defined for loadcase %d\n", lc);
     return;
  endif

  rc = 0;

# Get modal matrices of solid structure

  ndofr  = cmp.solid.stiff.ndofr;

  nofmod = cmp.solid.modes.nofmod;
  w      = cmp.solid.modes.omega;

  Mp = speye(nofmod);
  Kp = cmp.solid.modes.Kmodal;
  X  = cmp.solid.modes.disp;

  if (fsm)
     MgX = cmp.solid.mass.M * X;
  endif

  switch cmp.solid.damping.type
  case "ratios"
     nd = length(cmp.solid.damping.data);
     if (nd >= nofmod)
        d = cmp.solid.damping.data(1 : nofmod)';
     else
        dl = cmp.solid.damping.data(nd);
        d  = [cmp.solid.damping.data, dl(ones(1, nofmod -nd))]';
     endif
     Dp = spdiags(2 * d .* w, 0, nofmod, nofmod);
     if (fsm)
        DgX = MgX * Dp;
     endif
  case "Rayleigh"
     aK = cmp.solid.damping.data(1);
     aM = cmp.solid.damping.data(2);
     d  = aK * w .^2 + aM;
     Dp = spdiags(d, 0, nofmod, nofmod);
     if (fsm)
        DgX = aK * cmp.solid.stiff.K * X + aM * MgX;
     endif
  endswitch

# Get information on symmetry

  symy = isfield(cmp.aero, "symy");
  if (symy)
     ys = cmp.aero.symy;
  else
     ys = 0;
  endif

# Compute spline matrices

  ShX = cmp.splines.Shg * X;
  D1X = cmp.splines.D1h * ShX;
  D2X = 2 * cmp.splines.D2h * ShX;
  Svh = full(cmp.splines.Svh);

# Build load integrator matrices

  cref = cmp.aero.cref;
  [S1h, S2h] = mfs_vlmv2x(cmp.aero.panels, cmp.aero.pancols, Svh, cref);
  S1X = ShX' * S1h; S2X = ShX' * S2h;
  if (fsm)
     S1g = cmp.splines.Shg' * S1h; S2g = cmp.splines.Shg' * S2h;
  endif

# Get load data

  if (gust)
     qdyn = cmp.load.gust(lc).qdyn;
     v    = cmp.load.gust(lc).v;
     wg   = cmp.load.gust(lc).wg;
     x0   = cmp.load.gust(lc).x0;
     WG0  = -(wg / v) * cmp.aero.panels.nvec(3, :)';
     dx   = (i / v) * (x0 - cmp.aero.panels.C(1, :))';
  elseif (mnvr)
     qdyn = cmp.load.mnvr.qdyn(lc);
     v    = cmp.load.mnvr.v(lc);
     W1C  = cmp.load.mnvr.D1K * cmp.load.mnvr.UK(:, lc);
     W2C  = cmp.load.mnvr.D2K * cmp.load.mnvr.UK(:, lc);
  endif

# Add frequencies in halfpower bandwidth

  if (nb)
     fr = cmp.solid.modes.freq;
     [f, nfreq] = mfs_addfreq(f, fr, d, nb);
  else
     nfreq  = length(f);
  endif

# Build aerodynamic matrices that do not depend on reduced frequency
# and initialize wake grid

  [C, S, wake] = mfs_vlmhpre(cmp.aero.ls, cmp.aero.panels, cmp.aero.pancols,
                             symy, ys, cmp.aero.cref, opts);
  [CL, CU] = lu(C); clear C;
  D1X      = CU \ (CL \ D1X);
  D2X      = CU \ (CL \ D2X);
  unit     = eye(rows(S));
  if (mnvr)
     W1C   = CU \ (CL \ W1C);
     W2C   = CU \ (CL \ W2C);
  endif

# Initialize frequency loop

  ndofg  = cmp.solid.dofs.ndofg;
  nofpan = cmp.aero.panels.nofpan;

  kred = zeros(1, nfreq);
  Up   = zeros(nofmod, nfreq);
  GL   = zeros(nofpan, nfreq);
  if (fsm) 
    FS = zeros(ndofg, nfreq);
  endif

# Loop over frequencies

  for n = 1 : nfreq

     printf("           processing frequency %3d of %3d\r", n, nfreq);
     fflush(stdout);

     w = 2 * pi * f(n); kred(n) = (w * cref) / (2 * v);
     ik = i * kred(n);

     if (kred(n))
        DX = D1X + ik * D2X;
        SX = S1X + ik * S2X;
        if (fsm)
           Sg = S1g + ik * S2g;
        endif
     else
        DX = D1X;
        SX = S1X;
        if (fsm)
           Sg = S1g;
        endif
     endif

     if (gust)
        WL = WG0 .* exp(w * dx);  % dx includes imaginary unit i
        WL = CU \ (CL \ WL);
     elseif (mnvr)
        WL = W1C + ik * W2C;
     endif

     G = [DX, WL];

     if (kred(n))
        F  = mfs_vlmhwake(wake, kred(n), cmp.aero.panels.nofpan,
                          cmp.aero.pancols.npcol, opts.ktol);
        F  = CU \ (CL \ F);
        R  = unit - S * F;
        dG = R \ (S * G);
        G += F * dG;
     endif

     Qp = SX * G(:, 1 : nofmod);
     Ql = SX * G(:, nofmod + 1);

     Kdyn     = -w^2 * Mp + i * w * Dp + Kp - qdyn * Qp;
     rhs      = qdyn * Ql;
     Up(:, n) = Kdyn \ rhs;

     GL(:, n) = G * [Up(:, n); 1];

     if (fsm)
        Qg = qdyn * Sg * GL(:, n);
        FS(:, n) = (w^2 * MgX - i * w * DgX) * Up(:, n) + Qg;
     endif

  endfor

  printf("\n");

# Displacements from Force Summation Method

  if (fsm)

     ndofl = cmp.solid.dofs.ndofl;
     ndofd = cmp.solid.dofs.ndofd;
     ndofr = cmp.solid.stiff.ndofr;

     if (ndofr)
        dofe  = cmp.solid.stiff.dofe;
        dofeg = cmp.solid.dofs.dofl(dofe);
     else
        dofe  = cmp.solid.dofs.dofl;
        dofeg = dofe;
     endif

     Kll = mfs_matpart(cmp.solid.stiff.K, cmp.solid.dofs, [1, 0]);
     Kee = Kll(dofe, dofe);
     FSl = mfs_matpartr(FS, cmp.solid.dofs);
     FSe = FSl(dofe, :);
     Ue  = zeros(ndofg, nfreq);
     Ue(dofeg, :)  = Kee \ FSe;

     if (ndofd)
        dofd = cmp.solid.dofs.dofd;
        Ue(dofd, :) = cmp.solid.dofs.C(dofd, :) * Ue;
     endif

     if (ndofr)
        XR  = cmp.solid.modes.disp(:, 1 : ndofr);
        MXR = cmp.solid.mass.M * XR;
        Ue -= XR * (MXR' * Ue);
     endif

  endif

# Build results

  resps = struct("method", meth, "type", 11, "nfreq", nfreq, 
                 "kred", kred, "freq", f, "Q", Up);
  respa = struct("method", meth, "type", 11, "nfreq", nfreq, 
                 "kred", kred, "freq", f, "Q", Up, "G", GL);

  if (fsm)
     resps = setfield(resps, "Ue", Ue);
     Ua    = cmp.splines.Snh * (cmp.splines.Shg * Ue);
     respa = setfield(respa, "Ue", Ua);
  endif

endfunction
