%NEUTRAL_EX3
%
f1='%10.2f %5.0f %6.0f %21.15f\n';
global m imod dig
%disp('     alpha     m      n            I              error')
disp('     alpha     m      n            I')
dig=32; digits(dig);
N=100; m=0;
z=vpa(i/1000); y=imag(subs(z)); alpha=-1/2;
ab0=sr_jacobi01(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);
  fprintf(f1,alpha,m,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
  imab11=subs(abs(imag(ab(1,1))));
  maximagab=max(max(subs(abs(imag(ab(2:N,:))))));
  maxIMab=max(imab11,maximagab)
  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));
    fprintf(f1,alpha,m,n,I)
  end
  fprintf('\n')
end
