% EXAMPLE3_34 Computations for Example 3.34. 
%
clear mc mp iq idelta irout AB Z eps0 Nmax ab0 M
global mc mp iq idelta irout AB Z eps0 Nmax ab0 M
f1='%4.0f %4.0f %22.14e %12.4e %7.0f %7.0f\n';
disp('                           M=2*N')
%disp('                      M=2*floor((N+1)/2)')
%disp('                             M=4')
%disp('                             M=2')
fprintf('\n')
mc=1; mp=0; iq=1; idelta=2; irout=1; 
AB=[-1 1];
eps0=100*eps; Nm=10; Nmax=100; nu0=Nm+5; om=1.001;
ab0=r_jacobi(Nmax);
exact=12.9292568500027;
disp('   n    m        integral            error        Ncap    kount')
for N=1:Nm
  M=2*N;
%  M=2*floor((N+1)/2);
%  M=4;
%  M=2;
  M=M-2;
  sgn=1;
  for m=1:M
    sgn=-sgn;
    Z(m,1)=sgn/floor((m+3)/2); Z(m,2)=1;
  end
  [ab1,Ncap,kount]=r_mod(N,ab0);
  x1=om;
  hp=log(abs((x1+1)/(x1-1)));
  hn=-hp;
  if M>0
    for n=1:M
      p=1;
      for m=1:M
        if m~=n
          p=(1-Z(m,1)/Z(n,1))*p;
        end
      end
      c(n)=1/p;
    end
    hp=sum(c(1:M)'.*log(abs((x1+1)*(Z(1:M,1)+1)./((x1-1)*...
      (Z(1:M,1)-1))))./(1+x1*Z(1:M,1)));
    hn=sum(c(1:M)'.*log(abs((x1-1)*(Z(1:M,1)+1)./((x1+1)*...
      (Z(1:M,1)-1))))./(1-x1*Z(1:M,1)));
  end
  abmod=chri8(N,ab1,x1,eps0,nu0,Nmax,hp,hn,2);
  abmod(1,2)=-(om^2)*abmod(1,2);
  M=M+2;
  Z(M-1,1)=1/om; Z(M-1,2)=1;
  Z(M,1)=-1/om; Z(M,2)=1;
  xw=gauss_rational(N,abmod);
  int=sum(xw(:,2).*fex34(xw(:,1),om));
  err=abs((int-exact)/exact);
  if ~all(N-[1 4 7 10])
    fprintf(f1,N,M,int,err,Ncap,kount)
  end
end
