% testscript.m
%
% loads output from test1.f to check rate of convergence

load -force out1
load -force out2
load -force out4
load -force out8
load -force out16
load -force out32
load -force out64
load -force out128
load -force out256
load -force out512

 err = [mean(abs(out1(:,3) - out1(:,2))); ...
        mean(abs(out2(:,3) - out2(:,2))); ...
        mean(abs(out4(:,3) - out4(:,2))); ...
        mean(abs(out8(:,3) - out8(:,2))); ...
        mean(abs(out16(:,3) - out16(:,2))); ...
        mean(abs(out32(:,3) - out32(:,2))); ...
        mean(abs(out64(:,3) - out64(:,2))); ...
        mean(abs(out128(:,3) - out128(:,2))); ...
        mean(abs(out256(:,3) - out256(:,2))); ...
        mean(abs(out512(:,3) - out512(:,2)))];

%lng = length(out1);
%err = [out1(lng,4); ...
%       out2(lng,4); ...
%       out4(lng,4); ...
%       out8(lng,4); ...
%       out16(lng,4); ...
%       out32(lng,4); ...
%       out64(lng,4); ...
%       out128(lng,4); ...
%       out256(lng,4); ...
%       out512(lng,4)];


stp = 0.04 * 0.5.^(0:(length(err)-1))';

loglog(stp, err, 'o');
hold on;
%loglog(stp, err);
range = (1:length(err));
p = polyfit(log(stp(range)),log(err(range)),1);
xfit = [stp(range(1)),stp(range(length(range)))];
loglog(xfit, exp(polyval(p, log(xfit))));
hold off;
fprintf('fitted order = %f\n', p(1));