function [K,Kq,E] = gg303_profile(x,y)
% function [K,Kq,E] = gg303_profile(x,y)
% Plots profile for GG303 field trip data and a modified arctan profile
% Input parameters
% x = horizontal distance along porfile
% y = relative elevation
% Output values
% K = curvature from actual data by fitting of circles
% Kq = curvature of best fit modified arctan profile
% E = square of misfit error
% Example
% [x,y] = gg303_2003_profile_data;
% [K,Kq,E] = gg303_profile(x,y);

l = length(x);
m = 2:l-1;

K = curv(x,y);

% Find arctangent fit and the square of the misfit error
%q = 2*atan(0.3*(x+7))-1.5;     % This was my original "eye-ball best fit"

[a,b,c,d,x0,E] = gg303_profilefit(x,y);
q = a*1 + b*(x-x0) + c*atan((x-x0)/d);
Kq = curv(x,q);
E = (q-y)'*(q-y);

% plot the topographic and curvature profiles
subplot(3,1,1)
plot(x,y,x,q)
xlabel('Distance from base (m)')
ylabel('Relative Elevation (m)')
axis([min(x) max(x) min(y) max(y)]);
title('Topography')

subplot(3,1,2)
plot(x,y,x,q)
xlabel('Distance from base (m)')
ylabel('Relative Elevation (m)')
axis('equal');
title('Topography')

subplot(3,1,3)
plot(x(m),K,x(m),Kq)
xlabel('Distance from base (m)')
ylabel('Curvature (m^{-1})')
%axis([min(x) max(x) min(K) max(K)]);
title('Curvature')
legend('From actual data','From best fit atan curve')