function [cp, xv, cL, cM] = mfs_vortex2d(varargin)

# usage: [cp, xv, cL, cM] = mfs_vortex2d(varargin)
#
#    [cp, xv, cL, cM] = mfs_vortex2d(x, camber, alpha)
#    [cp, xv, cL, cM] = mfs_vortex2d("steady", x, camber, alpha)
#    [cp, xv, cL, cM] = mfs_vortex2d("harmonic", x, kred, motion)
#
# Input  x(:)       x-coordinates of interval boundaries
#        camber     Structure of piecewise polynomial defining camber
#        alpha(:)   List of angles of attack in degrees (optional)
#        kred(:)    List of reduced frequencies
#        motion     Structure defining motion
#
# Output cp(:, :)  Pressure coefficient (rows correspond to angles of
#                  attack or reduced frequencies)
#        xv(:)     Coordinates of vortex points
#        cL(:)     Lift coefficient
#        cM(:)     Moment coefficient with respect to quarter point
#
# The function computes the pressure coefficient, the lift coefficient
# and the moment coefficient 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)
#
# Pitch motion is about the quarter point. Positive pitch angle
# is clockwise.
#
# ----------------------------------------------------------------------

# Copyright (c) 2020 by Johannes Wandinger

# Check arguments

  if (length(varargin) < 2 || nargout != 4)
     print_usage();
  end

  if (ischar(varargin{1}))
     type = varargin{1};
     args = varargin(2 : end);
  else
     type = "steady";
     args = varargin;
  end

  nargs = length(args);

# Initialize output

  cp = 0; xv = 0; cL = 0; cM = 0;

# Check coordinates

  x = args{1};
  if (! isvector(x))
     error("mfs_vortex2d: x must be a vector\n"); 
  end
  [m, n] = size(x);
  if (n == 1)
     x = x';
  elseif (m != 1)
     error("mfs_vortex2d: x must be a vector\n");
  end

  if (! issorted(x))
     error("mfs_vortex2d: x not in ascending order\n");
  end

# Branch according to type

  switch type

  case "steady"

     if (nargs > 3)
        error("mfs_vortex2d: too many input arguments\n");
     end
     camber = args{2};
     if (! isstruct(camber))
        error("mfs_vortex2d: camber must be a structure\n"); 
     end
     if (nargs < 3)
        alpha = 0;
     else
        alpha = args{3};
        if (! isvector(alpha))
           error("mfs_vortex2d: alpha must be a vector\n"); 
        end
        alpha = alpha * (pi / 180);
     end

     mfs_paths("add", "mfs_vortex2d.m", "aero");
     [cp, xv, cL, cM] = mfs_vortex2ds(x, camber, alpha);
     mfs_paths("remove");

  case "harmonic"

     if (nargs != 3)
        error("mfs_vortex2d: wrong number of input arguments\n");
     end
     kred = args{2};
     if (! isvector(kred))
        error("mfs_vortex2d: kred must be a vector\n"); 
     end
     if (! isempty(find(kred < 0)))
        error("mfs_vortex2d: Reduced frequencies must be positive\n");
     end
     motion = args{3};
     if (! isstruct(motion))
        error("mfs_vortex2d: motion must be a strucgture\n");
     end
     if (isfield(motion, "flap"))
        if (! iscell(motion.flap))
           error("mfs_vortex2d: flap must be a cell array\n");
        end
        if (length(motion.flap) != 2)
           error("mfs_vortex2d: flap must consist of 2 cells\n");
        end
        if (motion.flap{1} < 0 || motion.flap{1} > 1)
           error("mfs_vortex2d: bad value of flap position (%f)\n",
                 motion.flap{1});
        end
     end

     mfs_paths("add", "mfs_vortex2d.m", "aero");
     [cp, xv, cL, cM] = mfs_vortex2dh(x, kred, motion);
     mfs_paths("remove");

  otherwise
     error("mfs_vortex2d: type %s not supported\n", type);
  end

end
