function [cmp, rc] = mfs_flutter_pk(cmp, v, rho, pairs, opts)

# usage: [cmp, rc] = mfs_flutter_pk(cmp, v, rho, pairs, opts)
#
# Input  cmp       Stucture with aeroelastic component
#        v(:)      List of velocities
#        rho       Mass density of air
#        pairs(:)  List of eigenvalue pairs to be considered 
#        opts      Structure with options
# Output cmp       Structure with aeroelastic component
#
# The function performs a flutter analysis using the pk-method.
#
# --------------------------------------------------------------------

# Copyright(c) 2024 by Johannes Wandinger

# Check arguments

  if (nargin != 5 || nargout != 2)
     print_usage();
  endif

# Initialize

  rc     = 0;
  nv     = length(v);
  v      = sort(v);
  q      = 0.5 * rho * v.^2;
  cref   = cmp.aero.cref;
  qk     = 0.5 * cref * q ./ v;
  npairs = length(pairs);
  nofmod = 2 * npairs;

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

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

# Build some constant matrices

  nx  = length(msel);
  IXX = eye(nx, nx);

# Get modal stiffness matrix

  KXX = cmp.solid.modes.Kmodal(msel, msel);

# Initialize results

  XX    = zeros(nx, nofmod, nv);
  pp    = zeros(nofmod, nv);
  kconv = zeros(1, nv);

# 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 velocities

  kmin  = eps;
  ktest = 0.5 * cmp.solid.modes.omega(pairs) * cref / v(1);

  for n = 1 : nv

     printf("           processing velocity %3d of %3d\n", n, nv);
     fflush(stdout);

     w2k   = 0.5 * cref / v(n);

     for np = 1 : npairs

         niter = 1;
         mode  = 2 * (pairs(np) - 1) + 1;
         m     = 2 * (np - 1) + 1;
         kt    = ktest(np);

         do

            if (opts.msg > 0)
               printf("\nMode pair %3.0d, iteration no. %3.0d:\n", np, niter);
            endif

            A = flutter_pkmat(nx, kt, q(n), qk(n), IXX, KXX, CL, CU, unit,
                              S, D1X, D2X, S1X, S2X, cmp.aero.panels.nofpan,
                              cmp.aero.pancols.npcol, wake, opts.ktol);
            [x, p] = eig(A);
            [w, ix] = sort(abs(diag(imag(p))));
            kred    = w2k * w(mode);

            conv = abs(kred - kt) < opts.ktol;

            if (! conv)
               if (niter > 1)
                  mk = (kred - k1) / (kt - kt1);
                  knew = (k1 - mk * kt1) / (1 - mk);
                  k1 = kred; kt1 = kt;
                  kt = knew;
               else
                  k1 = kred; kt1 = kt;
                  kt = kred;
               endif
            endif

         until (conv || ++niter > opts.mxiter)

         if (conv)
            printf("             mode %2d converged, k = %7.4f, niter = %2d\r",
                   mode, kred, niter);
            fflush(stdout);
            if (opts.msg > 0) printf("\n"); end
            pd = diag(p);
            ixc = [m, m + 1]; ixg = [mode, mode + 1];
            pp(ixc, n) = pd(ix(ixg));
            XX(:, ixc, n) = x(1 : nx, ix(ixg));
            kconv(np) = kred;
         else
            printf("\n*E* mfs_flutter: k-iteration did not converge\n");
            printf("    mode = %d, niter = %d\n", m, niter);
            rc = 1; return;
         endif

     endfor

     if (n < nv)
        ktest = kconv * v(n) / v(n + 1);
     endif
     printf("\n");

  endfor     % Loop over velocities

# Converged reduced frequencies and mode numbers

  kred  = 0.5 * cref * imag(pp) ./ v;
  modno = zeros(1, nofmod);
  modno([1 : 2 : nofmod]) = 2 * (pairs - 1) + 1; 
  modno([2 : 2 : nofmod]) = 2 * pairs;

# Sort results

  for n = 1 : nv
      for m = 1 : 2 : 2 * npairs
          if (abs(kred(m, n)) < kmin)
             if (real(pp(m, n) > real(pp(m + 1, n))))
                ix1 = [m + 1, m]; ix2 = [m, m + 1];
                pp(ix1, n) = pp(ix2, n);
                XX(:, ix1, n) = XX(:, ix2, n);
             endif
          else
             if (imag(pp(m, n) < 0))
                ix1 = [m + 1, m]; ix2 = [m, m + 1];
                pp(ix1, n) = pp(ix2, n);
                XX(:, ix1, n) = XX(:, ix2, n);
             endif
          endif
      endfor
  endfor

# Store results

  pimag = imag(pp);
  f     = pimag / (2 * pi);
  kred  = 0.5 * cref * pimag ./ v;

  cmp.flutter = struct("method", "pk", "nofmod", nofmod, "modno", modno, 
                       "nofpnt", nv, "nsel", nx, "kred", kred, "v", v,
                       "a", real(pp), "f", f, "msel", msel, "XX", XX);

endfunction

function A = flutter_pkmat(nx, kred, q, qk, IXX, KXX, CL, CU, unit, S,
                           D1X, D2X, S1X, S2X, nofpan, npcol, wake,
                           kzero)

# usage: A = flutter_pkmat(nx, kred, q, qk, IXX, KXX, CL, CU, unit,
#                          D1X, D2X, S1X, S2X, nofpan, npcol, wake,
#                          kzero)
#
# Input  nx         Dimension of modal basis
#        kred       Reduced frequency
#        q          Dynamic pressure
#        qk         0.5 * q * c / v
#        IXX        Unit matrix
#        KXX        Modal stiffness matrix
#        CL, CU     Lower and upper factor of C
#        unit       Unit matrix of size npcol by npcol
#        S          Vortex summation matrix
#        D1X, D2X   Modal downwash matrices multiplied by C-1
#        S1X, S2X   Modal load integrator matrices
#        nofpan     Number of panels
#        npcol      Number of panesl columns
#        wake       Structure with wake parameters
#        kzero      Zero threshold for reduced frequencies
#
# Output A          System matrix of eigenvalue problem
#
# The function builds the matrix of the eigenvalue problem to be
# solved in the pk-method.
#
# --------------------------------------------------------------------

# Initialize

  A = zeros(2 * nx, 2 * nx);

  ixul = 1 : nx;           % Indices of upper rows and left columns
  ixlr = nx + 1 : 2 * nx;  % Indices of lower rows and right columns

# Upper hyperrow

  A(ixul, ixlr) = -IXX;

# Lower hyperrow

  if (kred > kzero)

     ik = i * kred;

     G  = D1X + ik * D2X;
     F  = mfs_vlmhwake(wake, kred, nofpan, npcol, kzero);
     F  = CU \ (CL \ F);
     R  = unit - S * F;
     dG = R \ (S * G);
     G += F * dG;

     QXX = (S1X + ik * S2X) * G;

     A(ixlr, ixul) = KXX - q * real(QXX);
     A(ixlr, ixlr) = (qk / kred) * imag(QXX);

  else

     A(ixlr, ixul) = KXX - q * S1X * D1X;
     A(ixlr, ixlr) = qk * (S2X * D1X + S1X * D2X);

  end

endfunction
