%NEUTRAL_EX6
%
f0='%8.0f %21.15f %12.4e\n';
f1='%8.0f %21.15f\n';
global m imod dig
disp('       n            I              error')
%disp('       n            I')
dig=28;
digits(dig);
N=80; 
z=vpa(i/10); y=imag(subs(z));
for alpha=-1/4
  m=0;
  Iexact=.5*sqrt(pi)*(2*y)^(alpha+1/2)*gamma(alpha+1) ...
    *(struve(alpha+1/2,y)-bessely(alpha+1/2,y));
  ab0=sr_laguerre(dig,N);
  for n=20:20:N
    xw=gauss(n,subs(ab0));
    x=xw(:,1); w=xw(:,2);    
    I=sum(w.*(x.^2+y^2).^alpha);
    err=abs(I-Iexact);
%    I=sum(w.*((x.^2+y^2).^alpha).*besselj(0,x));
    fprintf(f0,n,I,err)
%    fprintf(f1,n,I)
  end
  fprintf('\n')
  for m=12
    imod=1;
    ab=mod_d(N,m,z,ab0,@rho01_mu);
    if imag(subs(z))~=0
      imod=2;
      ab=mod_d(N,m,z,ab,@rho01_mu);
    end
    maximagab=max(max(subs(abs(imag(ab)))))
    for n=20:20:N
      xw=gauss(n,subs(real(ab))); x=xw(:,1); w=xw(:,2);
      I=sum(w.*(x.^2+y^2).^(alpha+m));
      err=abs(I-Iexact);
%      I=sum(w.*((x.^2+y^2).^(alpha+m)).*besselj(0,x));
      fprintf(f0,n,I,err)
%      fprintf(f1,n,I)
    end
    fprintf('\n')
  end
  fprintf('\n')
end
