function [cp, xv, cL, cM] = mfs_vortex2ds(x, camber, alpha)

# usage: [cp, xv, cL, cM] = mfs_vortex2ds(x, camber, alpha)
#
# Input  x(:)      Coordinates of interval boundaries
#        camber    Piecewise polynomial structure defining camber
#        alpha(:)  List of angles of attack (in rad)
#
# Output cp(:, :)  Pressure coefficient
#        xv(:)     Coordinates of vortex points
#        cL(:)     Lift coefficients
#        cM(:)     Moment coefficients with respect to quarter point
#
# The function computes the steady aerodynamics of an airfoil using
# the discrete vortex method.
#
# ---------------------------------------------------------------------

# Check arguments

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

# Initialize output

  mr = length(alpha);
  nc = length(x) - 1;
  cp = zeros(mr, nc);
  xv = zeros(nc);
  cL = zeros(mr);
  cM = zeros(mr);

# Compute coordinates of vortex and control points

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

# Build the righthandside

  slope = ppder(camber);
  rhs   = 2 * pi * (ppval(slope, xc)' - alpha);

# Build matrix of influence coefficients

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

# Solve for vortices

  G = C \ rhs;

# Compute aerodynamic coefficients

  c  = x(end) - x(1);
  xq = x(1) + 0.25 * c;

  cL = 2 * sum(G) / c;
  cM = 2 * (xq - xv) * G / c^2;
  cp = 2 * G' ./ dx;

end
