function [lsout, panels, nodes, rc] = mfs_new_ls(msg, points, input_ls)

# usage: [lsout, panels, nodes, rc] = mfs_new_ls(msg, points, input_ls)
#
# Input  msg         File handle of message file
#        points      Structure with point data
#        input_ls    Structure array with definition of lifting surfaces
#
# Output lsout       Structure with lifting surface data
#        panels      Structure with vortex panel data
#        nodes       Structure with panel node data
#        rc          Return code: 0 = no errors
#                                 1 = errors
#
# The function builds structures ls, panels and nodes from the input
# data.
#
# The following fields are defined:
#
# ls.nofls                   Number of lifting surfaces
#   .ids(nofls)              Identifiers in ascending order
#   .surfs(nofls)            Lifting surfaces
#                .nx         Number of panels in x-direction
#                .ny         Number of panels in y-direction
#                .P(2, 3)    Coordinates of points P1 and P2
#                .c(2)       Chord at P1 and P2
#                .camber{2}  Camber at P1 and P2
#                .nvec(3)    Normal vector
#                .tx(:)      Parametric positions in x-direction
#                .ty(:)      Parametrci positions in y-direction
#                .p1         Index of first panel in structure panels
#                .pend       Index of last panel in structure panels
#
# panels.nofpan                 Number of panels
#       .corner(4, :)           Indicies of panel corner nodes
#       .A(:, 3)                Coordinates of vortex points A
#       .B(:, 3)                Coordinates of vortex points B
#       .C(:, 3)                Coordinates of control points C
#       .nvec(:, 3)             Normal vectors
#       .area(:)                Areas
#       .slope(:)               Slopes (from camber)
#
# nodes.nofnod                  Number of nodes
#      .ncoor                   Number of coordinates
#      .maxcor(3)               Maxima of the coordinates
#      .mincor(3)               Minima of the coordinates
#      .ids(:)                  Identifiers in ascending order
#      .coor(:, 3)              Coordinates
#
# --------------------------------------------------------------------

  leg_fields = {"id", "points", "chord", "alpha", "camber", ...
                "nx", "ny", "typex", "typey"};
  req_fields = {"id", "points", "chord", "nx", "ny"};

# Check arguments

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

# Initialize

  rc     = 0;
  lsout  = struct();
  panels = 0;
  nodes  = 0;

  ppzero = mkpp([0, 1], 0);

# Check input data

  if ((rc = mfs_checkfields1(msg, input_ls, req_fields, 
                             "Lifting Surfaces")))
     return;
  endif
  if ((rc = mfs_checkfields4(msg, input_ls, leg_fields, 
                             "Lifting Surfaces")))
     return;
  endif

  nofls    = length(input_ls);
  isalpha  = isfield(input_ls, "alpha");
  iscamber = isfield(input_ls, "camber");
  istypex  = isfield(input_ls, "typex");
  istypey  = isfield(input_ls, "typey"); 

# Check if identifiers of all lifting surfaces are defined

  for nls = 1 : nofls
      if (isempty(input_ls(nls).id))
         rc = 1;
         fprintf(msg, "*E* Identifier of lifting surface no. ");
         fprintf(msg, "%5d undefined\n", nls);
      endif
  endfor

  if (rc) return; end

# Check remaining data

  for ls = input_ls

      if (isempty(ls.points))
         rc = 1;
         fprintf(msg, "*E* Lifting surface %5d: ", ls.id)
         fprintf(msg, "Points undefined\n");
      elseif (length(ls.points) != 2)
         rc = 1;
         fprintf(msg, "*E* Lifting surface %5d: ", ls.id);
         fprintf(msg, "Incorrect definition of points\n");
      endif

      if (isempty(ls.chord))
         rc = 1;
         fprintf(msg, "*E* Lifting surface %5d: ", ls.id);
         fprintf(msg, "Chord length undefined\n");
      elseif (length(ls.chord) > 2)
         rc = 1;
         fprintf(msg, "*E* Lifting surface %5d: ", ls.id);
         fprintf(msg, "Incorrect definition of chords\n");
      endif

      if (isempty(ls.nx))
         rc = 1;
         fprintf(msg, "*E* Lifting surface %5d: ", ls.id);
         fprintf(msg, "Number of panels in x-direction undefined\n");
      endif

      if (isempty(ls.ny))
         rc = 1;
         fprintf(msg, "*E* Lifting surface %5d: ", ls.id);
         fprintf(msg, "Number of panels in y-direction undefined\n");
      endif

      if (isalpha && ! isempty(ls.alpha))
         if (length(ls.alpha) > 2)
            rc = 1;
            fprintf(msg, "*E* Lifting surface %5d: ", ls.id);
            fprintf(msg, "Incorrect definition of alpha\n");
         endif
      endif

      if (iscamber && ! isempty(ls.camber))
         if (iscell(ls.camber))
            if (length(ls.camber) > 2)
               rc = 1;
               fprintf(msg, "*E* Lifting surface %5d: ", ls.id);
               fprintf(msg, "Incorrect definition of camber\n");
            endif
         else
            if (length(ls.camber) != 1)
               rc = 1;
               fprintf(msg, "*E* Lifting surface %5d: ", ls.id);
               fprintf(msg, "Incorrect definition of camber\n");
            endif
         endif
      endif

  endfor

  if (rc) return; end

# Sort the lifting surfaces

  ids   = zeros(nofls, 1);

  for nls = 1 : nofls
      ids(nls) = input_ls(nls).id;
  endfor

  [ids, ixl] = sort(ids);

# Process the lifting surfaces

  for nls = 1 : nofls

      ls = input_ls(ixl(nls));

      surfs(nls).nx = ls.nx;
      surfs(nls).ny = ls.ny;
      if (istypex && ! isempty(ls.typex))
         typex{nls} = ls.typex;
      else
         typex{nls} = "cos";
      endif
      if (istypey && ! isempty(ls.typey))
         typey{nls} = ls.typey;
      else
         typey{nls} = "linear";
      endif

      surfs(nls).P = zeros(2, 3);
      ip = lookup(points.ids, ls.points(1), "m");
      if (ip)
         surfs(nls).P(1, :) = points.coor(ip, :);
      else
         rc = 1;
         fprintf(msg, "*E* Lifting surface %5d: ", ls.id);
         fprintf(msg, "Point %5d does not exist\n", ls.points(1));
      endif
      ip = lookup(points.ids, ls.points(2), "m");
      if (ip)
         surfs(nls).P(2, :) = points.coor(ip, :);
      else
         rc = 1;
         fprintf(msg, "*E* Lifting surface %5d: ", ls.id);
         fprintf(msg, "Point %5d does not exist\n", ls.points(2));
      endif

      c = ls.chord;
      if (length(c) == 1)
         c(2) = c(1);
      endif
      surfs(nls).c = c;
      if (c(1) < 0 || c(2) < 0)
         rc = 1;
         fprintf(msg, "*E* Lifting surface %5d: ", ls.id);
         fprintf(msg, "Chord length is negative\n");
      endif

      if (isalpha && ! isempty(ls.alpha))
         a = ls.alpha * pi / 180;
         if (length(a) == 1)
            a(2) = a(1);
         endif
         surfs(nls).a = a;
      else
         surfs(nls).a = zeros(1, 2);
      endif

      if (iscamber && ! isempty(ls.camber))
         if (iscell(ls.camber))
            camber = ls.camber;
         else
            camber{1} = ls.camber;
            camber{2} = ls.camber;
         endif
         surfs(nls).camber = camber;
      else
         surfs(nls).camber = {ppzero, ppzero};
      endif

      nvec = [0; surfs(nls).P(1, 3) - surfs(nls).P(2, 3);
                 surfs(nls).P(2, 2) - surfs(nls).P(1, 2)];
      lenn = norm(nvec);
      if (lenn > 0)
         surfs(nls).nvec = nvec / lenn;
      else
         rc = 1;
         fprintf(msg, "*E* Lifting surface %5d: ", ls.id);
         fprintf(msg, "Normal vector is zero\n");
      endif

  endfor

  if (rc) return; end

# Build the panels

  nofpan = 0;
  incnp  = 0;

  for nls = 1 : nofls

      irc = 0;

      switch typex{nls}
      case "linear"
         tx = linspace(0, 1, surfs(nls).nx + 1);
      case "cos"
         th = linspace(0, pi, surfs(nls).nx + 1);
         tx = 0.5 * (1 - cos(th));
      otherwise
         irc = 1;
         fprintf(msg, "*E* Lifting Surface %5d: ", ids(nls));
         fprintf(msg, "Unknown division type \"%s\"\n", typex{nls});
      endswitch

      switch typey{nls}
      case "linear"
         ty = linspace(0, 1, surfs(nls).ny + 1);
      case "cos"
         th = linspace(0, pi, surfs(nls).ny + 1);
         ty = 0.5 * (1 - cos(th));
      case "cos>"
         th = linspace(0.5 * pi, 0, surfs(nls).ny + 1);
         ty = cos(th);
      case "<cos"
         th = linspace(0, 0.5 * pi, surfs(nls).ny + 1);
         ty = 1 - cos(th);
      otherwise
         irc = 1;
         fprintf(msg, "*E* Lifting Surface %5d: ", ids(nls));
         fprintf(msg, "Unknown division type ""%s""\n", typey{nls});
      endswitch

      if (irc)
         rc = 1;
         continue;
      endif

      surfs(nls).tx = tx;
      surfs(nls).ty = ty;

      newnp = (surfs(nls).nx + 1) * (surfs(nls).ny + 1);

      PL = zeros(3, length(ty));     % Points along leading edge
      for k = 1 : 3
          PL(k, :) = (1 - ty) * surfs(nls).P(1, k) + ty * surfs(nls).P(2, k);
      endfor
      c  = (1 - ty) * surfs(nls).c(1) + ty * surfs(nls).c(2);
      dx = tx' * c;

      nx1 = surfs(nls).nx + 1;
      xx = repmat(PL(1, :), nx1, 1) + dx;
      yy = repmat(PL(2, :), nx1, 1);
      zz = repmat(PL(3, :), nx1, 1);
      x1 = reshape(xx, 1, newnp);
      y1 = reshape(yy, 1, newnp);
      z1 = reshape(zz, 1, newnp);
      xyz = [x1; y1; z1];

      tym  = 0.5 * (ty(1 : surfs(nls).ny) + ty(2 : surfs(nls).ny + 1));
      txC  = 0.25 * tx(1 : surfs(nls).nx);
      txC += 0.75 * tx(2 : surfs(nls).nx + 1);
      ppd1 = ppder(surfs(nls).camber{1});
      ppd2 = ppder(surfs(nls).camber{2});
      sl1  = ppval(ppd1, txC) - surfs(nls).a(1);
      sl2  = ppval(ppd2, txC) - surfs(nls).a(2);

      surfs(nls).p1 = nofpan + 1;

      for m = 1 : surfs(nls).ny

          for l = 1 : surfs(nls).nx

              nofpan++;

              ixc1 = l + (m - 1) * nx1;
              ixc2 = ixc1 + 1;
              ixc3 = l + m * nx1;
              ixc4 = ixc3 + 1;

              corner(1, nofpan) = ixc1 + incnp;
              corner(2, nofpan) = ixc2 + incnp;
              corner(3, nofpan) = ixc3 + incnp;
              corner(4, nofpan) = ixc4 + incnp;

              A(:, nofpan)  = 0.75 * xyz(:, ixc1);
              A(:, nofpan) += 0.25 * xyz(:, ixc2);
              B(:, nofpan)  = 0.75 * xyz(:, ixc3);
              B(:, nofpan) += 0.25 * xyz(:, ixc4);
              D(:, nofpan)  = 0.25 * xyz(:, ixc1);
              D(:, nofpan) += 0.75 * xyz(:, ixc2);
              E(:, nofpan)  = 0.25 * xyz(:, ixc3);
              E(:, nofpan) += 0.75 * xyz(:, ixc4);
              C(:, nofpan) = 0.5 * (D(:, nofpan) + E(:, nofpan));

              nvec(:, nofpan) = surfs(nls).nvec;

              area(nofpan) = mfs_tzarea(xyz(:, ixc1), xyz(:, ixc2), ...
                                        xyz(:, ixc3), xyz(:, ixc4), ...
                                        nvec(:, nofpan));

              slope(nofpan)  = (1 - tym(m)) * sl1(l);
              slope(nofpan) += tym(m) * sl2(l);

          endfor

      endfor

      surfs(nls).pend = nofpan;

      if (incnp)
         coor = [coor; xyz'];
      else
         coor = xyz';
      endif
      incnp = incnp + newnp;

  endfor

  dx = 2 * (C(1, :) -0.5 * (A(1, :) + B(1, :)));
  dy = area ./ dx;

  if (rc) return; end

# Build the structures

  maxcor = max(coor);
  mincor = min(coor); 
  idnod  = 1 : incnp;

  nodes = struct("nofnod", incnp, "ncoor", 3, "maxcor", maxcor,
                 "mincor", mincor, "ids", idnod, "coor", coor); 

  lsout = struct("nofls", nofls, "ids", ids, "surfs", surfs);

  panels = struct("nofpan", nofpan, "corner", corner, "A", A,
                  "B", B, "C", C, "nvec", nvec, "area", area,
                  "width", dy, "slope", slope);

endfunction
