function [Shg, Svh, Snh, D1h, D2h] = ...
         mfs_spline_tb(data, panels, anodes, ng, mxdofpnt, cref)

# usage: [Shg, Svh, Snh, D1h, D2h] = ...
#        mfs_spline_tb(data, panels, anodes, ng, mxdofpnt, cref)
#
# Input  data         Structure with spline data
#        panels       Structure with panel data
#        anodes       Structure with panel corner nodes
#        ng           Number of global dofs
#        mxdofpnt     Maximum number of dofs per nodal point
#        cref         Reference chord length
#
# Ouput  Shg(:, ng)   Contribution to matrix Shg
#        Svh(:, ;)    Contribution to matrix Svh
#        Snh(:, :)    Contribution to matrix Snh
#        D1h(:, ;)    Contribution to matrix D1h
#        D2h(:, ;)    Contribution to matrix D2h
#
# Matrix Shg relates the spline coefficients to the global
# displacements of the solid.
#
# Matrix Svh relates the displacements at the vortex points to the
# spline coefficients.
#
# Matrix Snh relates the displacements at the panel corner nodes to
# the spline coefficients.
#
# Matrix D1h relates the normal wash at the control points to the
# spline coefficients.
#
# Matrix D2h relates the displacements at the control points to the
# spline coefficients. It is only computed for cref > 0.
#
# --------------------------------------------------------------------

# Check arguments

  if (nargin != 6 || nargout != 5)
     print_usage();
  endif

# Initialize

  D2h    = [];
  pid    = data.pid;
  breaks = data.breaks;
  ixs    = data.ixs;
  T      = data.T;
  pA     = panels.A;
  pB     = panels.B;
  pC     = panels.C;
  corner = panels.corner;

# Build the spline base 

  nb = data.nbreaks;
  nh = 2 * nb;

  I      = eye(nb, nb);
  for n = 1 : nb
      B(n) = spline(breaks, I(:, n));
  endfor

# Build the system matrix

  nofnod = length(ixs);

  H = zeros(nofnod, nh);
  x = data.scoor(1, :);
  y = data.scoor(2, :);

  for n = 1 : nb
      H(:, n)      = ppval(B(n), y);
      H(:, nb + n) = -x' .* H(:, n);
  endfor

# Solve

  [Q, R] = qr(H, 0);
  Shz    = R \ Q';

# Matrix to get solid displacements normal to lifting surface

  if (mxdofpnt == 6)
     ez = [T(3, :), zeros(1, 3)];
  else
     ez = T(3, :);
  endif

  Pn = sparse(nofnod, ng);
  for n = 1 : nofnod
      ix1 = (ixs(n) - 1) * mxdofpnt + 1;
      ix2 = ix1 + mxdofpnt - 1;
      Pn(n, ix1 : ix2) = ez;
  endfor

# Matrix Shg

  Shg = Shz * Pn;

# Matrix Svh

  nv  = length(pid);
  Svh = zeros(nv, nh);
  x   = zeros(nv, 1);
  y   = zeros(nv, 1);
  n   = 1;

  for p = pid
      C = 0.5 * T * (pA(:, p) + pB(:, p));
      x(n)   = C(1);
      y(n++) = C(2);
  endfor

  for n = 1 : nb
      Svh(:, n)      = ppval(B(n), y);
      Svh(:, nb + n) = -x .* Svh(:, n);
  endfor

# Matrix D1h

  D1h = sparse(nv, nh);
  
  n = 1;
  for p = pid
      C = T * pC(:, p);
      x(n)   = C(1);
      y(n++) = C(2);
  endfor

  for n = 1 : nb
      D1h(:, nb + n) = -ppval(B(n), y);
  endfor

# Matrix D2h

  if (cref)

     D2h = sparse(nv, nh);

     for n = 1 : nb
         D2h(:, n)      = ppval(B(n), y) / cref;
         D2h(:, nb + n) = -x .* D2h(:, n);
     endfor

  endif

# Matrix Snh

  na = 3 * anodes.nofnod;
  Snh = zeros(na,  nh);

  node_used = zeros(1, anodes.nofnod);
  for p = pid
      for k = 1 : 4
          ix = corner(k, p);
          node_used(ix) = 1;
      endfor
  endfor
  nodix = find(node_used == 1);

  acoor = T * anodes.coor(nodix, :)';

  nz = length(nodix);
  uzh = zeros(nz, nh);

  for n = 1 : nb
      uzh(:, n)      = ppval(B(n), acoor(2, :));
      uzh(:, nb + n) = -acoor(1, :)' .* uzh(:, n);
  endfor

  Ez = sparse(na, nz);
  ez = T(3, :)';
  for n = 1 : nz
      ix1 = 3 * (nodix(n) - 1) + 1;
      ix2 = ix1 + 2;
      Ez(ix1 : ix2, n) = ez;
  endfor

  Snh = Ez * uzh;

endfunction
