clear all
close all
[X,Y] = ndgrid((-50:50),(-20:20));
Z = X+1i*Y;
Z = Z(:);
n = length(Z);
[r1,r2,r3,r4] =  exponential_integral(Z);
r = [r1(:),r2(:),r4(:),r3(:)];
load 'expint.ref';
err = abs(repmat(dref(:),1,4)-r)./abs(repmat(dref(:),1,4));
disp('          Z          series    Ramanuj cont-fr nr-cont-fr')
for i = 1:n
    disp(sprintf('%9.2e%s%9.2ei %.2e %.2e %.2e %.2e %.2e',...
		 real(Z(i)),cplx(imag(Z(i))),imag(Z(i)),...
		 err(i,1),err(i,2),...
		 err(i,3),err(i,4)))
end
disp('          Z          series    Ramanuj cont-fr nr-cont-fr')
[~,idx] = min(err,[],2);
hold on
c = {'*r','ob','xk','+g'};
l = {'series','Ramanujan','nr-cont-frac','cont-frac'};
counter = 0;
for i = 1:5
  if (any(idx==i))
    counter = counter+1;
    plot3(X(idx==i),Y(idx==i),idx(idx==i),c{i})
    leg{counter} = l{i};
  end
end
view(0,90)
legend(leg)
title('best results')
axis equal
