function cmp = mfs_flutter_k(cmp, kred, rho, opts)

# usage: cmp = mfs_flutter_k(cmp, kred, rho, opts)
#
# Input  cmp       Stucture with aeroelastic component
#        kred(:)   List of reduced frequencies
#        rho       Mass density of air
#        opts      Structure with options
# Output cmp       Structure with aeroelastic component
#
# The function performs a flutter analysis using the k-method.
#
# --------------------------------------------------------------------

# Copyright(c) 2024 by Johannes Wandinger

# Check arguments

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

# Initialize

  nk   = length(kred);
  np   = cmp.aero.panels.nofpan;
  r    = 0.5 * rho;
  cref = cmp.aero.cref;

  printf("           Flutter analysis using the k-method\n");

# Make sure reduced frequencies are ascending

  kred = sort(kred);

# Get information on symmetry

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

# Compute modal spline matrices

  ShX = cmp.splines.Shg * cmp.solid.modes.disp;
  D1X = cmp.splines.D1h * ShX;
  D2X = 2 * cmp.splines.D2h * ShX;
  SvX = cmp.splines.Svh * ShX;

# Build modal load integrator matrices

  [S1X, S2X] = mfs_vlmv2x(cmp.aero.panels, cmp.aero.pancols, SvX, cref);

# Identify modes that do not contribute

  MTM = (S1X + i * S2X) * (D1X + i * D2X);
  mtm = norm(MTM, 2, "columns") / norm(MTM, 2);
  mtol = opts.mtol * eps;
  msel = find(mtm > mtol);

  if (opts.msg)
     printf("mfs_flutter_k: mtol = %10.5e\n", mtol)
     disp("mtm = "); disp(mtm);
     disp("msel = "); disp(msel);
  endif

  D1X = D1X(:, msel); D2X = D2X(:, msel);
  S1X = S1X(msel, :); S2X = S2X(msel, :);

# Get modal structural matrices

  nx  = length(msel);
  KXX = cmp.solid.modes.Kmodal(msel, msel);
  MXX = (4 / cref^2) * eye(nx, nx);

# Initialize results

  XX = zeros(nx, nx, nk);
  p  = zeros(nx, nk);

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

# Loop over reduced frequencies

  for n = 1 : nk

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

     if (kred(n))
        ik = i * kred(n);
        G  = D1X + ik * D2X;
        SX = S1X + ik * S2X;
     else
        G  = D1X; 
        SX = S1X;
     endif

     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

     QXX = SX * G;
     BXX = kred(n)^2 * MXX + r * QXX;

     [x, my] = eig(KXX, BXX);

     p(:, n)     = diag(my);
     XX(:, :, n) = x / diag(norm(x, 2, "columns"));

  endfor

  printf("\n");

# Sort the results

  for n = 2 : nk

      mac = abs(XX(:, :, n-1)' * XX(:, :, n));
      [mm, ix] = max(mac);
      chk = find(diff(sort(ix)) == 0);

      if (isempty(chk))
         p(ix, n) = p(:, n);
         XX(:, ix, n) = XX(:, :, n);
      else
         printf("*W* mfs_flutter: mode pairing problem found ");
         printf("at k = %7.4f\n", kred(n));
      endif

  endfor

# Compute v, g and f

  g        = zeros(nx, nk);
  ix       = find(real(p(:, 1)));
  g(ix, :) = -imag(p(ix, :)) ./ real(p(ix, :));
  v1       = real(p) .* (1 + g.^2);
  v        = sqrt(abs(v1));
  f        = v * diag(kred) / (pi * cref);

# Store results

  modno = 1 : nx;

  cmp.flutter = struct("method", "k", "nofmod", nx, "modno", modno,
                       "nofpnt", nk, "nsel", nx, "kred", kred, "v", v,
                       "g", g, "f", f, "msel", msel, "XX", XX);

endfunction
