function [po,ERROR] = ortpol ( X , F , W , NTERMS  )

%*****************************************************************************80
%
%MATLAB version of SUBROUTINE ORTPOL on page 263f of conte/deboor
%**************   I N P U T   ************
%  (X(i), F(i)), i=1:length(X) are the abscissae and ordinates of
%         the given data points to be fitted
%  W  is the vector of weights to be used in the weighted scalar product
%  NTERMS-1  is the upper bound on the degree of the polynomial approximant  p  
%**************  O U T P U T  ************
%  po  is a structure containing the coefficients  D  of the polynomial  p  
%    with respect to an orthonormal polynomial basis, as well as the
%    arrays  B  and  C  for the three-term recurrence that generates that
%    orthonormal basis. 
%    The call  ORTVAL(po,x)  provides the value of  p  at  x.
%  ERROR(i)  is the difference between Y(i) and  p(X(i)).
%**************  M E T H O D  ************
%  The sequence  P0, P1, ..., PNTERMSm1 of orthogonal polynomials with
%  respect to the discrete inner product
%      (P,Q) := SUM ( P(X(I))*Q(X(I))*W(I), I=1:length(X) )
%  is generated in terms of their three-term recurrence
%     PJp1(x) = (x - B(J+1))*PJ(x) - C(J+1)*PJm1(x)
%  and the coefficient  D(J) of the weighted least squares approximation 
%  to the given data is obtained correspondingly as
%     D(J+1) = (F,PJ)/(PJ,PJ),  J=0:NTERMS-1.
%  Actually, in order to reduce cancellation,  (F,PJ)  is computed as
%  (ERROR,PJ), with  ERROR = F initially and, for each J , ERROR  is reduced
%  by  D(J+1)*P(J)  as soon as  D(J+1)  becomes available.
%
%  Reference:
%
%    Samuel Conte, Carl de Boor,
%    Elementary Numerical Analysis,
%    Third Edition,
%    SIAM, 2017,
%    ISBN: 978-1-611975-19-2.
%
W = W(:); X = X(:).'; F = F(:).'; B = X*W; C = 0; D = F*W; S = sum(W);
D = D/S; ERROR = F - D; 
if (NTERMS==1), po = struct('B',B,'C',C,'D',D); return, end
B = [B/S,zeros(1,NTERMS-1)]; PJM1 = ones(size(X)); PJ = X - B(1);
C = [C,zeros(1,NTERMS-1)]; 
D = [D,zeros(1,NTERMS-1)]; 
S = [S,zeros(1,NTERMS-1)]; 

for J=2:NTERMS
   P = PJ.*(W.'); D(J) = ERROR*(P.');
   P = P.*PJ; B(J) = X*(P.'); S(J) = sum(P);
   D(J) = D(J)/S(J);
   ERROR = ERROR - D(J)*PJ;
   if J==NTERMS, continue, end
   B(J) = B(J)/S(J); C(J) = S(J)/S(J-1);
   P = PJ; PJ = (X-B(J)).*PJ - C(J)*PJM1; PJM1 = P;
end
  po = struct('B',B,'C',C,'D',D); 

  return
end

