function [data, rc] = mfs_new_spline_tb(msg, splinedef, nodes, nset,
                                        surf, ptol)

# usage: [data, rc] = mfs_new_spline_tb(msg, splinedef, nodes, nset,
#                                       surf, ptol)
#
# Input  msg           File handle of message file
#        splinedef     Structure with spline definition
#        nodes         Structure with data of the solid nodes
#        nset          Structure with nodal point sets
#        surf(:)       Lifting surface data
#        ptol          Angular tolerance
#
# Output data          Structure with spline data
#        rc            Return code : 0 means no errors
#                                    1 means errors
#
# The function creates the spline data needed to compute a
# torsion-bending spline:
#
#   data.T(3, 3)      Transformation matrix to the spline system
#       .ixs(:)       Indices of solid nodes
#       .scoor(2, :)  Coordinates of solid nodes in spline coordinate
#                     system
#       .nbreaks      Number of spline breaks
#       .breaks(:)    Spline breaks
#       .pid(:)       Panel indices in structure panel
#
# --------------------------------------------------------------------

# Check arguments

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

# Initialize

  rc   = 0;
  data = struct();
  ns   = length(surf);

# Check if normal vectors of lifting surfaces are parallel

  nv1 = surf(1).nvec;
  for n = 2 : ns
      if ((1 - nv1' * surf(n).nvec) > ptol)
         rc = 1;
         fprintf(msg, "*E* Spline %2d: Lifting surface normals ",
                 splinedef.id)
         fprintf(msg, "are not parallel\n");
      end
  end

  if (rc) return; end
  
# Transformation matrix to the spline system

  ex = [1; 0; 0];
  T  = [ex, cross(nv1, ex), nv1]';

  data.T = T;

# Solid nodes indices

  if (isfield(splinedef, "nodes") && ! isempty(splinedef.nodes))
     if (ischar(splinedef.nodes))
        setnam = splinedef.nodes;
        if (! isfield(nset, setnam))
           fprintf(msg, "*E* Node set %s referenced from ", setnam);
           fprintf(msg, "spline %4d does not exist\n",
                   splinedef.id);
           rc = 1;
        else
           ixs    = nset.(setnam);
           nofnod = length(ixs);
        end
     else
        nodids = sort(splinedef.nodes);
        nofnod = length(nodids);
        ixs    = lookup(nodes.ids, nodids, "m");
        ixb    = find(ixs == 0);
        if (! isempty(ixb))
           for ix = ixb
               rc = 1;
               fprintf(msg, "*E* Solid node %6d referenced from ",
                       nodids(ix));
               fprintf(msg, "spline %4d does not exist\n",
                       splinedef.id);
           end
        end
     end
  else
     nofnod = nodes.nofnod;
     ixs = 1 : nodes.nofnod;
  end

  if (rc) return; end

  data.ixs = ixs;

# Transform the solid node coordinates to the spline system

  coor  = nodes.coor(ixs, :)';
  scoor = T(1 : 2, :) * coor;

  data.scoor = scoor;

# Spline breaks

  ysmax = max(scoor(2, :)); ysmin = min(scoor(2, :));

  for n = 1 : ns
     acoor(:, 1) = T(1 : 2, :) * surf(n).P(1, :)';
     acoor(:, 2) = T(1 : 2, :) * surf(n).P(2, :)';
     ymax = max([ysmax, acoor(2, :)]);
     ymin = min([ysmin, acoor(2, :)]);
  end

  if (isfield(splinedef.data, "breaks"))
     breaks  = sort(splinedef.data.breaks);
     nbreaks = length(breaks);
     if (ymin < breaks(1) || ymax > breaks(end) )
        rc = 1;
        fprintf(msg, "*E* Spline breaks of spline %4d do not cover ",
                splinedef.id);
        fprintf(msg, "coordinate range\n    [%12.5e, %12.5e]\n",
                ymin, ymax);
        return
     end
  elseif (isfield(splinedef.data, "nbreaks"))
     nbreaks = splinedef.data.nbreaks;
     breaks  = linspace(ymin, ymax, nbreaks);
  else
     rc = 1;
     fprintf(msg, "*E* Definition of spline breaks missing ");
     fprintf(msg, "for spline %4d\n", splinedef.id);
     return;
  end

  data.nbreaks = nbreaks;
  data.breaks  = breaks;

# Panel indices

  pid = [];
  for n = 1 : ns
      pid = [pid, surf(n).p1 : surf(n).pend];
  end
  data.pid = pid;

end
