function [a,b,c,d,x0,E] = gg303_profilefit(x,y)
% function [a,b,c,d,x0,E] = gg303_profilefit(x,y)
% Find parameters for a "best fit profile"
% The equation fit to the profile is:
% y(x) = a*1 + b*(x-x0) + c*atan((x-x0)/d)
% where 
% a = y-shift
% b = scaling of the surface slope
% c = scaling for the offset
% d = length scale for the scarp width
% x0 = x-offset of profile
% E = square of misfit matrix
% x = x-coordinates of survey points
% y = y-coordinates of survey points
% a,b, and c are scales for the basis functions
% 1, x-x0, and atan((x-x0)/d), respectively.

x = x(:);	% Ensure that x is a column vector
y = y(:);	% Ensure that y is a column vector

% Set reasonable ranges for x0 and d (the non-linear terms)
n = length(x);
x0 = -20:0.1:10;
%x0 = -20:1:10;
nx = length(x0);
d = 0.1:0.1:30;
%d = 1:1:30;
nd = length(d);
for i = 1:nx
    for j = 1:nd
        A = [ones(n,1) x-x0(i) atan((x-x0(i))./d(j))];
        P = A\y;
        e = A*P - y;
        E(i,j) = e'*e;      % square of the misfit error
    end
end

%find the values of i and j where E is a minimum. 
imin = find( min(E') == min(min(E')) );  % Find the row in which Emin occurs;
jmin = find( min(E) == min(min(E)) );    % Find the column in which Emin occurs;

%find x0 and d associated with the minimum value of E
x0 = x0(imin);
d = d(jmin);

%find a, b, and c
A = [ones(n,1) x-x0 atan((x-x0)./d)];
P = A\y;
a = P(1);
b = P(2);
c = P(3);