% CONDK Condition of the matrix K
%
global a y
f1='%8.0f %12.4e %12.4e %6.0f\n';
a=pi/5; r0=1; R=10;
eps0=.5e-8;
fprintf('eps0=%12.4e\n',eps0)
disp('       n     condK       max err0     max n')
B=10; 
ab=r_jacobi01(10);
for ng=2:3:8
%  xw=gauss(ng,ab);
%  x=xw(:,1);
  k=1:ng;
%  xg=r0+(R-r0)*x;
  xg=r0+(R-r0)*(k-1)/(ng-1);
  K=zeros(ng);
  err0m=0; nm=0;
  for i=1:ng
    for j=1:ng
      if j<=i
        x=xg(i); y=xg(j);
        phiB0=(4/pi)*(exp(-2*a*B)/(1-exp(-2*a*B)))*((.535/ ...
        (2*a))*y^(-3/4)*((B+1/(2*a))^2+1/(4*a^2))+(1/ ...
        (2*a+pi/2))*sqrt(2*pi/y)*exp(-y)*exp(-pi*B/2) ...
        *(B+1/(2*a+pi/2))); 
        phiB1=(4/pi)*(exp(-(2*a+pi/2)*B)/(1-exp(-a*B))) ...
        *((.535/(2*a+pi/2))*y^(-3/4)*(B+1/(2*a+pi/2)) ...
        +(1/(2*a+pi))*sqrt(2*pi/y)*exp(-y)*exp(-pi*B/2));
        [Phi,err0,n]=RMKLinvT(@fMKLinv,x,B,phiB0,phiB1,eps0);       
        K(i,j)=Phi;
        if j~=i, K(j,i)=Phi; end
        if err0>err0m, err0m=err0; end
        if n>nm, nm=n; end
      end
    end
  end
  condK=cond(K);
  fprintf(f1,ng,condK,err0m,nm)
end