function a = pwl_approx_1d_matrix ( nd, xd, yd, nc, xc )

%*****************************************************************************80
%
%% pwl_approx_1d_matrix() returns the matrix for the PWL approximant controls.
%
%  Discussion:
%
%    The value of the piecewise linear approximant, using control points XC
%    and control values YC, evaluated at the point XD, can be represented by
%
%      YD = A * YC
%
%    where A is a matrix whose values depend on XC and XD.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    08 January 2022
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer ND, the number of data points.
%    ND must be at least 1.
%
%    real XD(ND,1), the data points.
%
%    real YD(ND,1), the data values.
%
%    integer NC, the number of control points.
%    NC must be at least 1.
%
%    real XC(NC,1), the control points.
%
%  Output:
%
%    real A(ND,NC), the matrix.
%

%
%  Destroy all row vectors!
%
  xd = xd(:);
  yd = yd(:);
  xc = xc(:);

  a = zeros ( nd, nc );

  [ k, k ] = histc ( xd, xc );
%
%  Use extrapolation, if necessary.
%
  k ( k == 0 ) = 1;
  k ( k == nc ) = nc - 1;
  t = ( xd - xc(k,1) ) ./ ( xc(k+1,1) - xc(k,1) );
%
%  The obscure vectorized form of the following statements 
%  currently escapes me.
%
  for i = 1 : nd
    a(i,k(i)) = 1 - t(i);
    a(i,k(i)+1) = t(i);
  end

  return
end
