function [pancols, ls, rc] = mfs_new_pancols(msg, ls, panels, nodes, rtolp)

# usage: [pancols, ls, rc] = mfs_new_pancols(msg, ls, panels, nodes, rtolp)
#
# Input  msg        File handle of message file
#        ls         Structure with lifting surfaces
#        panels     Structure with panels
#        nodes      Structure with nodal points
#        rtolp      Relative tolerance for point coordinates
#
# Output pancols    Structure with information on panel columns
#        ls         Structure with lifting surfaces
#                   (panel column group identifier added)
#        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.
#
# A panel column group contains all panel columns sharing the same
# trailing edge.
#
# Fields of structure pancols:
#
#   npcol                Number of panel columns
#   ng                   Number of panel column groups
#   panids{npcol}        Cell array with lists of panels in each
#                        column
#   gcx(2, ng)           Index of first and last panel column within
#                        panel column group
#   tepcor(3, 2, npcol)  Coordinates of the two trailing edge points
#                        of the panel containing the trailing edge
#
# --------------------------------------------------------------------

# Check arguments

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

# Initialize

  pancols = struct();
  npcol   = 0;
  rc      = 0;

# Find lifting surfaces that are aligned in x-direction
# -----------------------------------------------------

# Initialize

  nofls = ls.nofls;
  plc   = zeros(nofls, 6);   % Leading each point coordinates
  ptc   = zeros(nofls, 6);   % Trailing edge point coordinates
  lsidx = zeros(2, nofls);   % Lifting surface indices
  lsact = ones(1, nofls);    % Active lifting surfaces

  lsidx(1, :) = 1 : nofls;

  surfs = ls.surfs;
  cmax  = 0;

  for n = 1 : nofls
      plc(n, :) = [surfs(n).P(1, :), surfs(n).P(2, :)];
      ptc(n, :) = plc(n, :);
      ptc(n, [1, 4]) += surfs(n).c;
      cmax      = max([cmax, surfs(n).c]);
  end

  tolp = cmax * rtolp;       % Tolerance for coincident points
  tol  = tolp^2;

# Find lifting surface behind lifting surface considered

  for n = 1 : nofls
      coort = ptc(n, :);
      for k = 1 : nofls
          if (k != n && lsact(k));
             d = coort - plc(k, :);
             if (d' * d < tol)
                lsidx(2, n) = k;
                lsact(k) = 0;
                break;
             end
          end
      end
  end

# Find lifting surface columns

  lsact = ones(1, nofls);
  ncol  = 0;    % counts lifting surface columns

  for n = 1 : nofls
      if (lsact(n))
         ncol++;
         lscol{ncol} = n; 
         nx = lsidx(2, n);
         while(nx)
            lscol{ncol} = [lscol{ncol}, nx];
            lsact(nx) = 0;
            nx = lsidx(2, nx);
         end
      end
  end

# Check compatibility of lifting surfaces in each column

  for n = 1 : ncol
      lsc = lscol{n};
      nls = length(lsc);
      ny  = surfs(lsc(1)).ny;
      for l = 2 : nls
          if (surfs(lsc(l)).ny != ny)
             ix1 = lsc(1); ix2 = lsc(l);
             fprintf(msg, "*E* Lifting surfaces %5d and %5d ",
                     ls.ids(ix1), ls.ids(ix2));
             fprintf(msg, "have incompatible y-divisons\n");                     
             rc = 1;
          end
      end
      npcol += ny;
  end

  if (rc) return; end

# Process panels within each column
# ---------------------------------

# Get panels within each column

  panids = cell(npcol, 1);
  gcx    = zeros(2, ncol);
  col1   = 1;

  for n = 1 : ncol

      lsc  = lscol{n};
      ix1  = lsc(1);
      ny   = surfs(ix1).ny;
      col2 = col1 + ny - 1;

      surfs(ix1).gid = n;

      pids = [];
      for l = 1 : length(lsc)
          ixls = lsc(l);
          p1   = surfs(ixls).p1;
          pend = surfs(ixls).pend;
          nx   = surfs(ixls).nx;
          pids = [pids; reshape(p1 : pend, nx, ny)];
          surfs(ixls).gid = n;
      end

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

      gcx(1, n) = col1; gcx(2, n) = col2;

      col1 = col2 + 1;

  end

  pancols.npcol  = npcol;  % Number of panel columns
  pancols.ng     = ncol;   % Number of panel groups
  pancols.panids = panids; % Panels within each column
  pancols.gcx    = gcx;    % Index of first and last panel column in group

# 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;

# Store updated liftung surface data (panel column group identifier added)

  ls.surfs = surfs;

end
