function F = mfs_vlmhwake(wake, kred, nofpan, npcol, kzero)

# usage: F = mfs_vlmhwake(wake, kred, nofpan, npcol, kzero)
#
# Input  wake     Structure array with wake information
#        kred     Reduced frequency
#        nofpan   Number of panels
#        npcol    Number of panel columns
#        kzero    Zero threhold for reduced frequencies
#
# Output F(:, :)  Matrix with influence coefficients of wake
#
# For each wake segment, structure wake contains the following data:
#
#   .nrows             Number of rows in matrix F
#   .ncols             Number of columns in matrix F
#   .rows(2)           Indices of first and last row in matrix C
#   .cols(2)           Indices of first and last column in matrix C
#   .s(nx)             Integration grid
#   .FS(nrows, ntep)   Factors of integral IS
#   .FF1(nrows, ncols) Factors of integral IF1
#   .FF2(nrows, ncols) Factors of integral IF2
#   .Y(:, nx)          Function values of IS, IF1 and IF2
#                      (if tmp = 0)
#   .ixY(:)            (1) Index of last entry of IS
#                      (2) Index of last entry of IF1
#                      (if tmp = 0)
#   .nsemi             Number of functions of semi-infinite vortex
#                      lines (if tmp > 0)
#   .nfini             Number of functions of finite vortex lines
#                      (if tmp > 0)
#
# Matrix F contains the contribution of the wake to the velocity at the
# control points.
#
# --------------------------------------------------------------------

# Copyright(c) 2024 by Johannes Wandinger

# Initialize F

  F = zeros(nofpan, npcol);

  if (kred < kzero) return; end

  w = 2 * kred; iw = i * w;

# Loop over wake segments and add contributions

  for wseg = wake

      nrows = wseg.nrows;
      m1    = wseg.rows(1); m2 = wseg.rows(2);
      ncols = wseg.ncols;
      n1    = wseg.cols(1); n2 = wseg.cols(2);
      iy1   = wseg.ixY(1); iy2 = wseg.ixY(2);

  % Perform the integrations

      I = mfs_qosfun1(wseg.s, wseg.Y, w);

  % Process contributions of semiinfinite vortex lines

      FS = wseg.FS .* reshape(I(1 : iy1), nrows, ncols + 1);

      F(m1 : m2, n1 : n2) -= iw * FS(:, 1 : ncols);
      F(m1 : m2, n1 : n2) += iw * FS(:, 2 : ncols + 1);

  % Process contributions of finite vortex lines

      FF  = wseg.FF1 .* reshape(I(iy1 + 1 : iy2), nrows, ncols);
      FF += wseg.FF2 .* reshape(I(iy2 + 1 : end), nrows, ncols);

      F(m1 : m2, n1 : n2) += iw * FF;

  endfor

endfunction
