function [camber, contour]  = mfs_airfoil(type, varargin)

# usage: pp = mfs_airfoil(type, varargin)
#
# Input  type       Type of airfoil
#        varargin   Airfoil data
#
# Output camber     Piecewise polynomial of camber
#        contour    Piecewise polynomial of contour
#
# The function returns the camber line of an aifoil as piecewise
# polynomial.
#
# Airfoil types:
#
# NACA:
#
#    varargin{1} = Series
#
#    Series = 4:  varargin{2} = digits 1 and 2
#                 varargin{3} = x1 (optional)
#                 varargin{4} = x2 (optional)
#
#    Series = 5:  varargin{2} = digits 2 and 3
#                 varargin{3} = x1 (optional)
#                 varargin{4} = x2 (optional)
#
# fit:
#
#    varargin{1} = subtype      "camber" or "contour"
#    varargin{2} = data(:, 2)   x- and z-coordinates of camber or
#                               contour
#    varargin{3} = splopts      either number of spline intervals
#                               or structure with spline options
#    varargin{4} = x1           Start position (optional)
#    varargin{5} = x2           End position (optional)
#
# Start and end positions are given as fractions of the chord length.
# The default is 0 for the start positions and 1 for the end position.
#
# Structure splopts has the following fields:
#
#   nbcam   Number of camber spline intervals
#   nbcon   Number of contour spline intervals
#   smooth  Number of loops to smooth data
#
# Smoothing, if requested, is done by averaging over 3 adjacent
# points.
#
# If the last argument is "plot", the camber line is plotted. The
# plot also contains the fitted points and the contour if available.
#
# --------------------------------------------------------------------

# Copyright (c) 2017 by Johannes Wandinger

# Initialize

  camber   = struct();
  contour  = struct();
  rc       = 0;
  cut      = 0;
  plotflag = 0;
  dataflag = 0;

# Check arguments

  if (nargin < 2 || nargout > 2)
     print_usage();
  end

  if (! ischar(type))
     error("mfs_airfoil: first argument must be a string\n");
  end

# Check if plots are requested

  nargs = length(varargin);

  if (strcmp(varargin(nargs), "plot"))
     plotflag = 1;
     nargs   -= 1;
  end
  args = varargin(1 : nargs);

# Branch on airfoil type

  switch type

  case "NACA"

     if (nargs > 2)
        if (nargs < 4)
           error("mfs_airfoil: incomplete definition of x1 and x2\n");
        else
           x1  = args{3}; x2 = args{4};
           if (x1 < 0 || x2 < x1 || x2 > 1)
              error("mfs_airfoil: bad definition of x1 or x2\n");
           end
           dx = x2 - x1;
           cut = 1;
        end
     end

     nplot = 40;

     switch args{1}

     case 4
       if (nargs < 2)
          error("mfs_airfoil: NACA series %d: insufficient data supplied\n",
                args{1});
       end
       mfs_paths("add", "mfs_airfoil.m", "aero");
       camber = mfs_naca4(args{2});
       mfs_paths("remove");
       if (cut)
          D = diag([dx^2, dx, 1]);
       end

     case 5
       if (nargs < 2)
          error("mfs_airfoil: NACA series %d: insufficient data supplied\n",
                args{1});
       end
       mfs_paths("add", "mfs_airfoil.m", "aero");
       [camber, rc] = mfs_naca5(args{2});
       mfs_paths("remove");
       if (cut)
          D = diag([dx^3, dx^2, dx, 1]);
       end

     otherwise
        error("mfs_airfoil: NACA series %d not supported\n", args{1});
     end

  case "fit"

     if (nargs < 3)
        error("mfs_airfoil: insufficient data supplied for option \"fit\"\n");
     end

     if (! ischar(args{1}))
        error("mfs_airfoil: type \"fit\": second argument must be a string\n");
     end

     data = args{2};
     if (! ismatrix(data))
        error("mfs_airfoil: third argument must be a matrix\n");
     end
     if (columns(data) != 2)
        error("mfs_airfoil: data matrix must have 2 columns\n");
     end

     if (nargs > 3)
       if (nargs < 5)
           error("mfs_airfoil: incomplete definition of x1 and x2\n");
        else
           x1  = args{4}; x2 = args{5};
           if (x1 < 0 || x2 < x1 || x2 > 1)
              error("mfs_airfoil: bad definition of x1 or x2\n");
           end
           dx  = x2 - x1;
           D   = diag([dx^3, dx^2, dx, 1]);
           cut = 1;
        end
     end

     splopts = args{3};

     switch args{1}

     case "camber"
       mfs_paths("add", "mfs_airfoil.m", "aero");
       [camber, rc] = mfs_fitcamber(data, splopts);
       mfs_paths("remove");
       nplot = 4 * rows(data);

     case "contour"
       mfs_paths("add", "mfs_airfoil.m", "aero");
       [camber, contour, rc] = mfs_fitcontour(data, splopts);
       mfs_paths("remove");
       nplot = 2 * rows(data);

     otherwise
        error("mfs_airfoil: airfoil subtype \"%s\" not supported\n", args{1});
     end

     dataflag = 1;

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

  if (rc)
     error("mfs_airfoil ended with errors\n");
  end

# Cut the camber if requested

  if (cut)
     cb     = camber.breaks;
     ix     = lookup(cb, [x1, x2]);
     ixs    = ix(1) : min(ix(2) + 1, camber.pieces + 1);
     b      = (cb(ixs) - x1) / dx;
     c      = camber.coefs(ixs(1 : end - 1), :) * D;
     c      = c / dx;
     camber = mkpp(b, c);
  end

# Plot data if requested

  if (plotflag)

     phi = linspace(0, pi, nplot);
     xc  = 0.5 * (1 - cos(phi));
     zc  = ppval(camber, xc);

     hold on

     plot(xc, zc, "color", "red");

     if (! cut)
         if (dataflag)
            data = data / data(end, 1);
            plot(data(:, 1), data(:, 2), "color", "green", 
                 "marker", "o", "linestyle", "none");
            if (numfields(contour))
               nplot = 2 * nplot;
               phi   = linspace(0, 2 * pi, nplot);
               x     = 0.5 * (1 + cos(phi));
               z     = ppval(contour, phi);
               plot(x, z, "color", "blue");
            end
         end
     end

     hold off

     grid;
     xlabel("x/c"); ylabel("z/c");

  end

end
