function [cp, xv, cL, cM] = mfs_vortex2dh(x, kred, motion)

# usage: [cp, xv, cL, cM] = mfs_vortex2dh(x, kred, motion)
#
# Input  x(:)      Coordinates of interval boundaries
#        kred(:)   List of reduced frequencies
#        motion    Structure defining motion
#
# Output cp(:, :)  Pressure coefficient
#        xv(:)     Coordinates of vortex points
#        cL(:)     Lift coefficients
#        cM(:)     Moment coefficients
#
# The function computes the harmonic aerodynamics of an airfoil using
# the discrete vortex method.
#
# Fields of structure motion:
#
#    Name    Type        Meaning
#  ------------------------------------------------------------------
#    heave   complex     Amplitude of heave motion
#    pitch   complex     Amplitude of pitch motion in radians
#    flap    cell array  (1) Starting position of flap in percent
#                            of chord (real)
#                        (2) Amplitude of flap angle in radians
#                            (complex, positive down)
#
# ---------------------------------------------------------------------

# Check arguments

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

# Initialize output

  mk = length(kred);      % Number of reduced frequencies
  nc = length(x) - 1;     % Number of intervals
  cp = zeros(mk, nc);
  xv = zeros(nc);
  cL = zeros(mk);
  cM = zeros(mk);

# Compute coordinates of vortex and control points

  dx = diff(x);
  xv = x(1 : nc) + 0.25 * dx;
  xc = x(1 : nc) + 0.75 * dx;

# Initialize some variables

  c = x(end) - x(1);    % Chord length
  G = zeros(nc, mk);    % Array with vortices: cols = reduced freq.

# Build the right-hand side

  rhs = zeros(nc, mk);

  if (isfield(motion, "heave"))
     rhs = 2 * i * (motion.heave / c) * ones(nc, 1) * kred;
  endif
  if (isfield(motion, "pitch"))
     rhs -= motion.pitch * (1 + 2 * i * (xc' / c - 0.25) * kred);
  endif
  if (isfield(motion, "flap"))
     xf  = c * motion.flap{1};
     ixf = find(xc >= xf);
     rhs -= motion.flap{2} * (1 + 2 * i * (xc(ixf)' - xf) * kred / c);
  endif

  rhs = 2 * pi * rhs;

# Build the system matrix of the steady discrete vortex method

  C = 1 ./ (xv - xc');

# Build the wake matrices

  ikc  = 2 * i * kred / c;
  xred = (c - xc') * ikc;
  F    = exp(xred) .* expint(xred) * diag(ikc);

# Prepare solution

  WF = C \ [rhs, F];
  W  = WF(:, 1 : mk);
  F  = WF(:, mk + 1 : end);
  sw = sum(W);
  sf = sum(F);
  G  = W;

# Loop over reduced frequences

  for m = 1 : mk
      if (kred(m) > 0)
         r = 1 - sf(m);         
         G(:, m) += F(:, m) * (sw(m) / r);
      endif
  endfor

# Compute aerodynamic coefficients

  Gs  = [zeros(1, mk); cumsum(G)];
  dxT = dx';
  cpT = 2 * (G ./ dxT + 2 * i * Gs(1 : nc, :) * diag(kred) / c);
  cp  = transpose(cpT);

  xq = x(1) + 0.25 * c;

  cL = (cp * dxT) / c;
  cM = (((xq - xv) .* cp) * dxT) / c^2;

endfunction
