% EXAMPLE3_33  Computations for Example 3.33.
%
clear mc mp iq idelta irout AB Z eps0 Nmax ab0 M
global mc mp iq idelta irout AB Z eps0 Nmax ab0 M
f2='%4.0f %5.0f %24.14e %14.4e %7.0f %7.0f\n';
%disp('                           M=2*N')
%disp('                            M=N')
disp('                             M=2')
%disp('                              M=1')
%disp('                             M=0')
fprintf('\n')
mc=1; mp=0; iq=1; idelta=2; irout=1; 
AB=[0 1];
eps0=100*eps; Nm=10; Nmax=100; om=.001;
ab0=r_jacobi01(Nmax,0,-.5);
exact=96.7036882293812;
disp('   n     m          integral              error        Ncap    kount')
for N=1:Nm
%  M=2*N;
%  M=N;
   M=2;
%  M=1;
%  M=0;
  M=M-1;
  if M>0
    for m=1:M
      Z(m,1)=1/m; Z(m,2)=1;
    end
  end
  [ab1,Ncap,kount]=r_mod(N,ab0);
  x1=-om;
  rho0=-2*atan(1/sqrt(abs(x1)))/sqrt(abs(x1));
  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
%
% the following statement, including the computation of the preceding rho0,
% should preferably be done in higher precision because of severe
% cancellation in the summation.
%
    rho0=sum(c(1:M)'.*(rho0+2*sqrt(Z(1:M,1)).*atan(sqrt(Z(1:M,1))))./ ...
      (1+x1*Z(1:M,1)));
  end
  if M<0
    xw=gauss(N,ab1);
    M=M+1;
  else
    z=complex(x1,0); nu0=0;
    abmod=chri4(N,ab1,z,eps0,nu0,Nmax,rho0,0);
    abmod(1,2)=om*abmod(1,2);
    M=M+1;
    Z(M,1)=1/om; Z(M,2)=1;
    xw=gauss_rational(N,abmod);
  end
  int=sum(xw(:,2).*fex33(xw(:,1),om));
  err=abs((int-exact)/exact);
  if ~all(N-[1 4 7])
    fprintf(f2,N,M,int,err,Ncap,kount)
  end
end
