digits(dig);
Nmax=100; if gam<=1/2, eps0=100*eps0; end
sc=vpa(sqrt(1+(gam/2)^2));
c=subs(sc);
sab0=sr_jacobi01(dig,Nmax,-1/2,0);
r0=vpa((2/gam)*sqrt(1+sc));
rho0=vpa(2*r0*atan(r0));
if gam<=1/2, ab=schri4(dig,Nmax,sab0,sc,rho0); end
N=0;
I0=0; I=1;
err=abs(I-I0);
while err>eps0
  N=N+5; 
  if N>Nmax
    error('N exceeds Nmax')
  end
  I0=I;
  if gam>1/2
    xw=gauss(N,subs(sab0)); x=xw(:,1); w=xw(:,2);
    if n==2*floor(n/2)
      fun=-.25*(w.*(chebpol(x,n)./((x+c).*sqrt(1+x)))./ ...
          (x-c));
    else
      fun=.25*(w.*(sinh(bet*x)./(1+cosh(bet*x))).* ...
          (chebpol(x,n)./((x+c).*sqrt(1+x)))./(x-c));
    end
  else
    xw=gauss(N,subs(ab)); x=xw(:,1); w=xw(:,2);
    if n==2*floor(n/2)
      fun=-.25*(w.*(chebpol(x,n)./((x+c).*sqrt(1+x))));
    else
      fun=.25*(w.*(sinh(bet*x)./(1+cosh(bet*x))).* ...
          (chebpol(x,n)./((x+c).*sqrt(1+x))));
    end
  end
  I=sum(fun); M=max(abs(fun));  
  can=log10(M/abs(I));
  err=abs(I-I0);
  if abs(I)>1, err=err/abs(I); end
end
