function [camber, contour, rc] = mfs_fitcontour(data, splopts)

# usage: [camber, contour, rc] = mfs_fitcontour(data, splopts)
#
# Input  data(:, 2)  x- and z-coordinates of points on camber
#        splopts     either number of spline intervals or
#                    structure with spline options
#
# Output camber      Piecewise polynomial defining camber
#        contour     Piecewise polynomial defining contour
#        rc          Return code: 0 = no errors
#
# Structure splopts has the following fields:
#
#   nbcam   Number of camber spline intervals
#   nbcon   Number of contour spline intervals
#   smooth  Number of loops to smooth camber
#
# The contour parameter is the angle, starting with 0 at the end and
# ending with 2*pi also at the end. Angles less than pi correspond to
# points on the upper side and angles larger than pi to points on the
# lower side of the airfoil.
#
# ---------------------------------------------------------------------

# Check arguments

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

# Initialize

  camber  = struct();
  contour = struct();
  rc      = 0;
  np      = rows(data);
  smooth  = 0;

# Check spline options

  if (isstruct(splopts))
     if (! isfield(splopts, "nbcam"))
        printf("*E* mfs_airfoil: number of camber spline intervals undefined\n");
        rc = 1; return;
     else
        nbcam = splopts.nbcam;
     end
     if (isfield(splopts, "nbcon"))
        nbcon = splopts.nbcon;
        if (nbcon <= 0)
           printf("*E* mfs_airfoil: number of contour spline intervals must be positive\n");
           rc = 1; return;
        end
        if (nbcon > np)
           printf("*E* mfs_airfoil: number of contour spline intervals (%d)\n", 
                  nbcon);
           printf("                 exceeds number of data points (%d)\n", np);
           rc = 1; return;
        end
        if (nbcam > nbcon)
           printf("*E* mfs_airfoil: number of camber spline intervals (%d)\n",
                  nbcam);
           printf("                 exceeds number of contour spline intervals (%d)\n", 
                  nbcon);
           rc = 1; return;
        end
     else
        nbcon = min([2 * nbcam, np]);
     end
     if (isfield(splopts, "smooth"))
        smooth = splopts.smooth;
     end
  else
     nbcam = splopts;
     nbcon = min([2 * nbcam, np]);
  end

  if (nbcam <= 0)
     printf("*E* mfs_airfoil: number of camber spline intervals must be positive\n");
     rc = 1; return;
  end

# Scale the data

  data   = data / data(end, 1);

# Find spline fitting the contour

  theta = acos(2 * data(:, 1) - 1);
  dx    = diff(data(:, 1));
  ix    = find(dx > 0) + 1;
  theta(ix) = 2 * pi - theta(ix);

  contour = splinefit(theta, data(:, 2), nbcon);

# Find camber

  xca = sort(data(:, 1));
  dxc = diff(xca);
  dxm = max(dxc);
  tol = dxm * 1e-2;
  ixc = find(dxc > tol);
  xc  = [xca(ixc); xca(end)];
  xc(1) = 0;

# Find spline fitting the camber

  psio = acos(2 * xc - 1);
  psiu = 2 * pi - psio;
  zc   = 0.5 * (ppval(contour, psio) + ppval(contour, psiu));

  if (smooth)
     for k = 1 : smooth
         zc(2 : end - 1) = (zc(1 : end - 2) + zc(2 : end - 1) + ...
                            zc(3 : end)) / 3;
     end
  end

  camber = splinefit(xc, zc, nbcam);

end
