function [pancols, rc] = mfs_new_wake(msg, wake, ls, panels, nodes)

# usage: [pancols, rc] = mfs_new_wake(msg, wake, ls, panels, nodes)
#
# Input  msg        File handle of message file
#        wake       Structure array with definition of wake groups
#        ls         Structure with lifting surfaces
#        panels     Structure with panels
#        nodes      Structure with nodal points
#
# Output pancols    Structure with information on panel columns
#        rc         Return code: 0 = no errors
#                                1 = errors
#
# The function builds the information on panel columns that is needed
# to process the wake in unsteady aerodynamics. A panel column contains
# all adjacent panels with the same y-coordinates, from the leading
# edge to the trailing edge.
#
# Fields of structure pancols:
#
#   npcol                Number of panel columns
#   tepcor(3, 2, npcol)  Coordinates of the two trailing edge points
#                        of the panel containing the trailing edge
#   panids{npcol}        Cell array with lists of panels in each
#                        column
#
# --------------------------------------------------------------------

  fields = {"id", "ls"};   % Legal and required fields

  deps = eps * 1e3;    % Distance tolerance

# Check arguments

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

# Initialize

  pancols = struct();
  npcol   = 0;

# Check fields

  if ((rc = mfs_checkfields1(msg, wake, fields, "Wakes")))
     return;
  end
  if ((rc = mfs_checkfields4(msg, wake, fields, "Wakes")))
     return;
  end

# Get lifting surface indices

  ngroup = length(wake);

  for ng = 1 : ngroup
      if (isempty(wake(ng).id))
         fprintf(msg, "*E* wake(%d): identifier is missing\n", ng);
         rc = 1;
         continue;
      end
      lsix = lookup(ls.ids, wake(ng).ls, "m");
      if (isempty(lsix))
         fprintf(msg, "*E* Wake group %5d: ", wake(ng).id);
         fprintf(msg, "Definition of lifting surfaces missing\n");
         rc = 1;
         continue
      end
      ixb  = find(lsix == 0);
      if (! isempty(ixb))
         for k = 1 : length(ixb)
             fprintf(msg, "*E* Wake group %5d: ", wake(ng).id);
             fprintf(msg, "Lifting surface %5d does not exist\n",
                     wake(ng).ls(ixb(k)));
         end
         rc = 1;
      end
      wake(ng).ls = lsix;
  end

  if (rc) return; end

# Check if the lifting surfaces are correctly aligned

  for ng = 1 : ngroup
      nls = length(wake(ng).ls);
      ny  = ls.surfs(wake(ng).ls(1)).ny;
      for l = 2 : nls
          ix1 = wake(ng).ls(l - 1); 
          ix2 = wake(ng).ls(l);
          if (ls.surfs(ix2).ny != ny)
             fprintf(msg, "*E* Wake group %5d: ", wake(ng).id);
             fprintf(msg, "Lifting surfacess have incompatible ");
             fprintf(msg, "y-divisions (%d != %d)\n",
                     ny, ls.surfs(ix2).ny);
             rc = 1;
          end
          tol = max(ls.surfs(ix1).c) * deps;
          P1 = ls.surfs(ix1).P;
          P2 = ls.surfs(ix2).P;
          d = norm(P2 - P1, "rows")' - ls.surfs(ix1).c;
          if (max(abs(d)) > tol)
             fprintf(msg, "*E* Wake group %5d: ", wake(ng).id);
             fprintf(msg, "Lifting surfaces not correctly aligned\n");
             rc = 1;
          end
      end
      npcol += ny;
  end

  if (rc) return; end

# Get panels within each column

  panids  = cell(npcol, 1);

  col1 = 1;

  for ng = 1 : ngroup

      ix1  = wake(ng).ls(1);
      ny   = ls.surfs(ix1).ny;
      col2 = col1 + ny - 1;

      pids = [];
      for l = 1 : length(wake(ng).ls)
          ixls = wake(ng).ls(l);
          p1   = ls.surfs(ixls).p1;
          pend = ls.surfs(ixls).pend;
          nx   = ls.surfs(ixls).nx;
          pids = [pids; reshape(p1 : pend, nx, ny)];
      end

      panids(col1 : col2)  = num2cell(pids, 1);

      col1 = col2 + 1;

  end

  pancols.npcol   = npcol;
  pancols.panids  = panids;

# Get trailing edge points

  tepcor = zeros(3, 2, npcol);

  for n = 1 : npcol
      pids = panids{n};
      ixp  = pids(end);
      ixc1 = panels.corner(2, ixp);
      ixc2 = panels.corner(4, ixp);
      tepcor(:, 1, n) = nodes.coor(ixc1, :)';
      tepcor(:, 2, n) = nodes.coor(ixc2, :)';
  end

  pancols.tepcor = tepcor;

end
