function [yvec_new] = nordsieck_intpol(xvec,yvec,xvec_new)
%
%  This function performs 3rd-order Nordsieck interpolation for
%  equidistantly spaced data.  An error results if the data are not
%  equidistantly spaced.
%
%  We first set the T-matrix.  We use Eq.(4.96) from the book.
%
T = [   6   0   0   0
       11 -18   9  -2
        6 -15  12  -3
        1  -3   3  -1 ] / 6;
%
%  We check whether the data are equidistantly spaced.
%
[n,m] = size(xvec);
nm = n*m;
d = xvec(2:nm) - xvec(1:nm-1);
if min(d) ~= max(d),
  error('Data not equidistantly spaced')
 else
  dx = d(1);
end
x0 = xvec(1);
%
%  We normalize the data.
%
xvec = 0:nm-1;
xvec_new = (xvec_new - x0*ones(size(xvec_new)))/dx;
%
%  We loop over all new x values.
%
[n1,m1] = size(xvec_new);
nm1 = n1*m1;
yvec_new = zeros(n1,m1);
for i=1:nm1,
  x = xvec_new(i);
  %
  %  We start by determining, which four values to use in the
  %  interpolation.  We work with the indices.  We distinguish
  %  between the first segment (c=1), the last segment (c=3),
  %  and all intermediate segments (c=2).
  %
  ix = floor(x);
  if ix == 0,
    indx = 4:-1:1;
    c = 1;
   elseif ix >= nm-2,
    indx = nm:-1:nm-3;
    c = 3;
  else
    indx = ix+3:-1:ix;
    c = 2;
  end
  %
  %  Now we need to determine the (normalized) dxnew value.
  %
  d = x - ix;
  if c == 1,
    if d < 0.5,
      dxn = (3 - d)/3;
    else
      dxn = (3 - d)/2;
    end
   elseif c == 2,
    if d < 0.5,
      dxn = (2 - d)/2;
     else
      dxn = 2 - d;
    end
   else
    dxn = 1 - d;
  end
  %
  %  We compute the H matrix.
  %
  dxn2 = dxn*dxn;
  dxn3 = dxn2*dxn;
  H = diag([ 1 , dxn , dxn2 , dxn3 ]);
  %
  %  Now, we change the sampling rate.
  %
  v = yvec(indx);
  vnew = T\H*T*v;
  %
  %  Finally, we interpolate.
  %
  if c == 1,
    if d < 0.5,
      y = vnew(4);
     else
      y = vnew(3);
    end
   elseif c == 2,
    if d < 0.5,
      y = vnew(3);
     else
      y = vnew(2);
    end
   else
    y = vnew(2);
  end
  if x == nm-1,
    y = vnew(1);
  end
  yvec_new(i) = y;
end
%
return