function [cmps, cmpa, nfound] = mfs_diverg(cmp, nvec, options)

# usage: [cmps, cmpa, nfound] = mfs_diverg(cmp, nvec, options)
#
# Input  cmp     Structure with data of aeroelastic component
#        nvec    Number of vectors
#        options Structure with options (optional)
#
# Output cmps    Structure with data of solid component
#        cmpa    Structure with data of aerodynamic component
#        nfound  Number of positive real eigenvalues found
#
# The function computes the dynamic pressure at divergence and the
# corresponding displacements.
#
# The eigenvalue problem is solved for nvec eigenvectors. Subsequently,
# positive real eigenvalues are searched for.
#
# ---------------------------------------------------------------------

# Copyright (c) 2024 by Johannes Wandinger

  t0 = clock();

# Check arguments

  if (nargin < 2 || nargin > 3 || nargout < 2 || nargout > 3)
     print_usage();
  endif

  if (! isstruct(cmp))
     error("mfs_diverg: first argument must be a structure\n");
  endif
  if (! isscalar(nvec) || nvec <= 0)
     error("mfs_diverg: second argument must be a positive scalar\n");
  endif 

# Check component type and status

  if (! strcmp(cmp.type, "aeroelastic"))
     error("mfs_diverg: component type must be \"aeroelastic\"\n");
  endif
  if (! isfield(cmp, "splines"))
     error("mfs_diverg: definition of splines missing\n");
  endif
  if (! isfield(cmp.splines, "Shg"))
     error("mfs_diverg: spline matrices missing\n");
  endif

# Initialize

  cmps   = cmp.solid;
  cmpa   = cmp.aero;
  nfound = 0;

  time_stamp = ignore_function_time_stamp("all");
  mfs_paths("add", "mfs_diverg.m", "aero", "aeroelastic", "util");

# Get information on degrees of freedom

  ndofl = cmp.solid.dofs.ndofl;
  ndofd = cmp.solid.dofs.ndofd;
  dofl  = cmp.solid.dofs.dofl;
  if (ndofd)
     dofd  = cmp.solid.dofs.dofd;
  endif

# Get stiffness matrix

  Kll = mfs_matpart(cmp.solid.stiff.K, cmp.solid.dofs, [1, 0]);

# Get aerodynamic matrix

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

  C = mfs_vlmmat(cmp.aero.ls, cmp.aero.panels, symy, ys); 

  G   = 2 * diag(cmp.aero.panels.width);
  qvh = G * (C \ cmp.splines.D1h);
  qhh = cmp.splines.Svh' * qvh;

# Solve the eigenvalue problem

  Shl = mfs_matpartc(cmp.splines.Shg, cmp.solid.dofs);
  Klh = Kll \ (Shl' * qhh);

  AF = @(x) Klh * (Shl * x);
  [vec, eval]  = eigs(AF, ndofl, nvec, "lm");

# Extract physically relevant data

  d     = diag(eval);
  neval = length(d);
  ix    = find(imag(d) == 0);

  nreal = length(ix);

  if (nreal)
     qda    = 1 ./ real(d(ix));
     veca   = real(vec(:, ix));
     ix     = find(qda > 0);
     nfound = length(ix);
     if (nfound)
        qd = qda(ix); 
        [qd, ix] = sort(qd);
        vecd = veca(:, ix);
     endif
  endif

# Store results

  divergs = struct("ndiv", nfound);
  diverga = divergs;

  if (nfound)

     divergs.qdyn = qd;
     divergs.disp = zeros(cmp.solid.dofs.ndofg, nfound, "double");
     divergs.disp(dofl, :) = vecd;
     if (ndofd)
        divergs.disp(dofd, :) = ...
           cmp.solid.dofs.C(dofd, :) * divergs.disp;
     endif

     diverga.qdyn = qd;
     diverga.disp = cmp.splines.Snh * (cmp.splines.Shg * divergs.disp);

  endif

  cmps.diverg = divergs;
  cmpa.diverg = diverga;

# Terminate

  printf("           Eigenvalues of static divergence analysis:\n");
  printf("             Found = %2.0d, Real = %2.0d, Positive = %2.0d\n",
         neval, nreal, nfound);

  mfs_paths("remove");

  ignore_function_time_stamp(time_stamp);

  elapsed_time = etime(clock(), t0);
  printf("%10.4f seconds needed to perform divergence analysis ",
         elapsed_time)
  printf("of component %s\n", inputname(1));
  fflush(stdout);

endfunction
