function [resps, respa, rc] = mfs_dfreq(cmp, f, lc, opts)

# usage: [resps, respa, rc] = mfs_dfreq(cmp, f, lc, opts)
#
# Input  cmp      Structure with aeroelastic component
#        f(:)     List with excitation frequencies
#        lc       Load case number
#
# Output resps   Structure with solid frequency response
#        respa   Structure with aerodynamic frequency response
#        rc      Return code: 0 means no errors
#
# ------------------------------------------------------------------------

# Copyright(c) 2024 by Johannes Wandinger

# Check arguments

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

# Initialize

  resps = struct(); respa = struct();
  rc    = 1;

# Check availability of data needed

  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
  if (! strcmp(cmp.solid.damping.type, "Rayleigh"))
     printf("*E* mfs_freqresp: only Rayleigh damping supported in ")
     printf("direct frequency response analysis\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 matrices of solid structure

  dofs  = cmp.solid.dofs;
  ndofg = dofs.ndofg;
  ndofl = dofs.ndofl;
  ndofd = dofs.ndofd;

  Kll = mfs_matpart(cmp.solid.stiff.K, dofs, [1, 0]);
  Mll = mfs_matpart(cmp.solid.mass.M, dofs, [1, 0]);
  Shl = mfs_matpartc(cmp.splines.Shg, dofs);

  aK  = cmp.solid.damping.data(1);
  aM  = cmp.solid.damping.data(2);
  Dll = aK * Kll + aM * Mll;

# Get information on symmetry

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

# Compute spline matrices

  D1h = cmp.splines.D1h;
  D2h = 2 * cmp.splines.D2h;
  Svh = full(cmp.splines.Svh);
  nh  = columns(D1h);

# Build load integrator matrices

  cref = cmp.aero.cref;
  [S1h, S2h] = mfs_vlmv2x(cmp.aero.panels, cmp.aero.pancols, Svh, cref);

# Get load data

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

  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

# 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;
  D1h      = CU \ (CL \ D1h); 
  D2h      = CU \ (CL \ D2h); 
  unit     = eye(rows(S));
  if (mnvr)
     W1C   = CU \ (CL \ W1C); 
     W2C   = CU \ (CL \ W2C); 
  endif

# Initialize frequency loop

  nfreq  = length(f);
  nofpan = cmp.aero.panels.nofpan;

  kred = zeros(1, nfreq);
  Ul   = zeros(ndofl, nfreq);
  GL   = zeros(nofpan, nfreq);

  Ihh  = eye(nh);

# 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))
        Dh = D1h + ik * D2h;
        Sh = S1h + ik * S2h;
     else
        Dh = D1h;
        Sh = S1h;
     endif

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

     G = [Dh, 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

     Qhh = Sh * G(:, 1 : nh);
     qh  = Sh * G(:, end);

     Kdyn     = -w^2 * Mll + i * w * Dll + Kll;
     KinvSlh  = Kdyn \ Shl';
     Qlh      = KinvSlh * Qhh;
     rhs      = qdyn * (KinvSlh * qh);

     Rhh       = Ihh - (qdyn * Shl) * Qlh;
     Ul(:, n)  = rhs + Qlh * (Rhh \ (Shl * (qdyn * rhs)));

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

  endfor

  printf("\n");

  U = zeros(ndofg, nfreq);
  U(dofs.dofl, :) = Ul;
  if (ndofd)
     U(dofs.dofd, :)  = dofs.C(dofs.dofd, :) * U;
  endif

# Build results

  resps = struct("method", 1, "type", 11, "nfreq", nfreq, 
                 "kred", kred, "freq", f, "nback", nfreq, 
                 "freqback", f, "U", U);
  Ua    = cmp.splines.Snh * (cmp.splines.Shg * U);
  respa = struct("method", 1, "type", 11, "nfreq", nfreq, 
                 "kred", kred, "freq", f, "nback", nfreq, 
                 "freqback", f, "U", Ua, "G", GL);

endfunction
