% GLOBATTO_JACOBI Generalized Gauss-Lobatto quadrature formula
%
%    Given the parameters al and be of the Jacobi measure and the 
%    fixed nodes -1 and +1, both of multiplicity r, this routine 
%    generates the generalized Gauss-Lobatto quadrature rule 
%    involving r consecutive derivative values at the fixed points 
%    endl and endr and N interior points. The N interior nodes
%    and weights are stored in the last N rows of the (N+2r)x2
%    output array xw, and the fixed nodes endl and endr (repeated
%    r times) and the corresponding 2r weights associated with 
%    the r derivative values in the first 2r rows of xw. It is 
%    assumed that N>2 and r>0.
%
%    The routine is designed to avoid breakdown at very large 
%    values of N, as it would occur in the routine globatto for 
%    ab=r_jacobi.
%
function xw=globatto_jacobi(N,al,be,r)
if r<1, error('r too small'), end
if N<3, error('N too small'), end
apb=al+be+2*r;
odd=2*floor(r/2)~=r;
%
% internal nodes and weights
%
ab1=r_jacobi(N,al+r,be+r);
xw0=gauss(N,ab1);
xw(2*r+1:N+2*r,:)=xw0;
xw(2*r+1:N+2*r,2)=xw(2*r+1:N+2*r,2)./((1-xw(2*r+1:N+2*r,1) ...
  .^2).^r);
xw(1:r,1)=-ones(r,1);
xw(r+1:2*r,1)=ones(r,1);
%
% boundary weights
%
r0=floor((r+1)/2); rs=zeros(r-1,1); tau=zeros(N+r0,1);
tau(1:N)=xw0(:,1);
for e=1:2
  if e==1
    a=-1; b=1; tau(N+1:N+r0)=b; c=be;
  else
    a=1; b=-1; tau(N+1:N+r0)=b; c=al;
  end
  p=1;
  for n=1:N
    p=(1+(c+r)/n)*p;
  end
  A=zeros(r); rhs=zeros(r,1); x=zeros(r,1);
  for rho=1:r-1
    rs(rho)=sum((tau-a).^(-rho));
    if odd, rs(rho)=rs(rho)-.5*(b-a)^(-rho); end
  end
  for i=r:-1:1
    A(i,i)=factorial(i-1)*p^2*prod((tau(N+1:N+r0)-a).^2);
    if odd, A(i,i)=A(i,i)/(b-a); end
    if i<r
      for j=i+1:r
        A(i,j)=-2*sum(A(i+1:j,j).*rs(1:j-i))/(j-i);
      end
    end
  end
  ng=N+r;
  ab=r_jacobi(ng,al,be);
  xwg=gauss(ng,ab);
  for i=1:r
    s=0;
    for ig=1:ng
      t=xwg(ig,1); w=xwg(ig,2);
      Pn=.5*(apb+2)*t+.5*(al-be); Pnm1=1;
      for n=2:N
        n2=2*n;
        Pnm2=Pnm1; Pnm1=Pn;
        Pn=((n2+apb-1)*((n2+apb)*(n2+apb-2)*t ...
          +(al+r)^2-(be+r)^2)*Pnm1 ...
          -2*(n+al+r-1)*(n+be+r-1)*(n2+apb)*Pnm2) ...
          /(n2*(n+apb)*(n2+apb-2));
    end
    p=prod((tau(N+1:N+r0)-t).^2);
    if odd, p=p/(b-t); end 
    s=s+w*((t-a)^(i-1))*p*Pn^2;
    end
    rhs(i)=s;
  end
  x(r)=rhs(r)/A(r,r);
  for i=r-1:-1:1
    x(i)=(rhs(i)-sum(A(i,i+1:r).*(x(i+1:r))'))/A(i,i);
  end
  xw(1+(e-1)*r:e*r,2)=x;
end