% Three versions of the quarter circle, with Tres Hermanos.

format long  

pow = [1,2,3,4];
Nsteps = (3*10.^pow)
Ntrials = numel(pow) 

O = '---- circle ---- = O'
for i = 1:Ntrials
   t0 = cputime;
   results(i) = thint(@circle, 0, 1, Nsteps(i));
   execution_time = cputime - t0
end
errors = results(2:end) - results(1:end-1);
   % this error estimate presumes that there is much less error in
   % the next trial, which has 10 times the number of points evaluated. 
estimated_error = errors(end)
estimate_pi = 4*results(end)
actual___pi = pi
abs_error = estimate_pi - pi
hold off
plot(log10(Nsteps(1:end-1)), log10(abs(errors)), '+-b');
hold on

O = '---- circle2 (transformed from circle) ---- = O'
for i = 1:Ntrials
   t0 = cputime;
   results(i) = thint(@circle2, 0, 1, Nsteps(i));
   execution_time = cputime - t0
end
errors = results(2:end) - results(1:end-1);
   % this error estimate presumes that there is much less error in
   % the next trial, which has 10 times the number of points evaluated. 
estimated_error = errors(end)
estimate_pi = 4*( results(end) + 2*sqrt(2)/3 )
actual___pi = pi
abs_error = estimate_pi - pi
plot(log10(Nsteps(1:end-1)), log10(abs(errors)), '^-g');

O = '---- circle3 (transformed from circle) ---- = O'
for i = 1:Ntrials
   t0 = cputime;
   results(i) = thint(@circle3, 0, 1, Nsteps(i));
   execution_time = cputime - t0
end
errors = results(2:end) - results(1:end-1);
   % this error estimate presumes that there is much less error in
   % the next trial, which has 10 times the number of points evaluated. 
estimated_error = errors(end)
estimate_pi = 4*( results(end) + 17/15/sqrt(2) )
actual___pi = pi
abs_error = estimate_pi - pi
plot(log10(Nsteps(1:end-1)), log10(abs(errors)), 'o-m');

plot(log10(Nsteps(1:end-1)), log10(Nsteps(1:end-1).^(-3))-3.5, '--k');
plot(log10(Nsteps(1:end-1)), log10(Nsteps(1:end-1).^(-4))-3, '--k');

xlabel('log_{10} number of steps', 'fontsize', 16);
ylabel('log_{10} absolute error', 'fontsize', 16);
title('Calculating \pi', 'fontsize', 18);
legend('circle','circle2','circle3','y \propto N^{-3}','y \propto N^{-4}','location','bestoutside');
% print('../plots/test3-pi.pdf','-dpdf')

hold off