function parms = linfit(x, y, res, A0, A1, B0, B1, its);
%Creates and graphs a matrix (As Figure 1) giving \chi^2
%as a function of A and B parameters in the formula
%f(x) = Ax + B.
%Syntax: linfit(x, y, res, A0, A1, B0, B1, its)
%   x, y    column vectors of data points
%   res     grid resolution
%   A0..B1  parametre range to check
%   its     number of recursions to take
% example:
% load logdata;
% linfit(logdata(:,1), logdata(:,2), 128, -10, 0, 0, 50, 3)

if(res <= 1)
    error('Resolution must be greater then one.');
end
if((its > 10) | (res > 512))
    warning('This might be a bit extreme. Remember: Ctrl-C to abort.');
end

a_inc = (A1 - A0) / (res - 1);
b_inc = (B1 - B0) / (res - 1);
disp('Interpolating between a0 - a1 and b0 - b1 with r steps.');
disp([A0, A1, B0, B1, res]);
chimtx(res, res) = 0;
b = B0;
for j = 1:res
    a = A0;
    for i = 1:res
        chimtx(i, j) = sum((a*x + b - y) .^ 2);
        a = a + a_inc;
    end
    b = b + b_inc;
end
clear u;
u(res,1) = 0;
for i = 1:res
    u(i, 1) = A0 + a_inc * (i - 1);
end
clear v;
v(res,1) = 0;
for j = 1:res
    v(j, 1) = B0 + b_inc * (j - 1);
end
[min_as, a_indices] = min(chimtx);
[chi2, b_index] = min(min_as);
a_index = a_indices(1, b_index);
a_fit = u(a_index,1);
b_fit = v(b_index,1);
disp('Best fit point to Ax + B in the range was');
disp([a_fit, b_fit]);
disp('With \chi^2 of');
disp(chi2);
if its > 0
    linfit(x, y, res, a_fit - a_inc, a_fit + a_inc, b_fit - b_inc, b_fit + b_inc, its - 1);
else
    figure(1);
    surface(u, v, chimtx);
    colormap(hot);
    shading interp;
    axis tight;
    view(3);
    %alpha interp;
    %alpha(gradient(z));
    title '\chi^2 vs. Linear Parametres AX + B';
    xlabel 'A';
    ylabel 'B';
    zlabel '\chi^2';
end

% return parametres
parms = [a_fit, b_fit, chi2];
