function ab=sr_Fermi(i)
%SR_FERMI Recurrence coefficients for polynomials orthogonal 
%with respect to Fermi type weight functions.
%   SR_FERMI(I) generates to 32 decimal digits the first
%   80 recurrence coefficients for polynomials orthogonal
%   with respect to the weight function W(T)=1/(EXP(T)+1)
%   on [0,Inf] if I=1, for polynomials orthogonal with
%   respect to the weight function W(T)=[1/(EXP(T)+1)]^2
%   on [0,Inf] if I=2, and for polynomials orthogonal with
%   respect to the weight function W(T)=T^{-1/2}/(EXP(T)+1)
%   on [0,Inf] if I=3, and places them into the 80x2 array AB.
   
syms mom ab 
if i==1
  digits(98); dig=98;
  mom(1)=log(vpa(2,dig));
  for k=2:160
    mom(k)=(1-2^(-vpa(k-1,dig)))*gamma(vpa(k,dig))*zeta(vpa(k,dig));
  end
elseif i==2
  digits(135); dig=135;
  mom(1)=log(vpa(2,dig))-1/2;
  mom(2)=vpa(pi,dig)^2/12-log(vpa(2,dig));
  for k=3:160
    mom(k)=gamma(vpa(k,dig))*((1-2^(vpa(1-k,dig)))*zeta(vpa(k,dig)) ...
    -(1-2^(vpa(2-k,dig)))*zeta(vpa(k-1,dig)));
  end
elseif i==3
  digits(97); dig=97;
  for k=1:160
    mom(k)=(1-2^(vpa(-k+3/2,dig)))*gamma(vpa(k-1/2))*zeta(vpa(k-1/2));
  end
else
  error('i not in range')
end 
ab=schebyshev(dig,80,mom);
ab=vpa(ab,32);
