%NEUTRAL_EX5
%
f0='%8.0f %4.0f %14.3e\n';
f1='%8.0f %4.0f %20.14f\n';
%disp('       n    m        error')
disp('       n    m            I') 
global imod m dig
imod=1;
for dig=48
  digits(dig);
  z=-vpa(1/1000); alpha=-1/2;
  Iexact=gamma(alpha+1)*exp(-subs(z))*gammainc(-subs(z),alpha+1,'upper');
  N=200; m=0;
  ab0=sr_laguerre(dig,N);
  for n=40:40:N
    xw=gauss(n,subs(ab0));
    x=xw(:,1); w=xw(:,2);
    I=sum(w.*(x-subs(z)).^alpha);
%    I=sum(w.*((x-subs(z)).^alpha).*atan(x));
    err=abs(I-Iexact);
    fprintf(f0,n,m,err)
%    fprintf(f1,n,m,I)
  end
  fprintf('\n')
  imod=1;
  for m=10:10:50
%  for m=50
    ab0=sr_laguerre(dig,N);
    ab=mod_d(N,m,z,ab0,@rho01_mu);
    ab=subs(ab);
    for n=40:40:N
      xw=gauss(n,ab);
      x=xw(:,1); w=xw(:,2);
      I=sum(w.*(x-subs(z)).^(m+alpha));
%      I=sum(w.*((x-subs(z)).^(m+alpha)).*atan(x));
      err=abs(I-Iexact);
      fprintf(f0,n,m,err)
%      fprintf(f1,n,m,I)
    end
    fprintf('\n')
  end
end
