function [plot1handles,plot2handles,plot3handles] = scarpdater_plotfs(x,z,thetadeg,a,vec1best,vec2best,vec3best,vec1min,vec2min,vec3min,vec1max,vec2max,vec3max,axistoken)

% Makes a three row subplot displaying the results of a sensitivity analysis for
%   The input and output variables are as follows:
%   x is the value of x for which each scarp elevation z is defined.
%   z is topographic elevation.
%   vec1best is the value of the vec1 value that yields the lowest RMS value.
%   vec2best is the value of the vec2 value that yields the lowest RMS value.
%   vec3best is the value of the vec3 value that yields the lowest RMS value.
%   * NOTE: vec1best, vec2best, and vec3best must all be valid numbers that
%     are uniquely contained in vec1, vec2, and vec3, respectively.
%   vec1min is the minimum vec1 value to be plotted in sensitivity analysis.
%   vec2min is the minimum vec2 value to be plotted in sensitivity analysis.
%   vec3min is the minimum vec3 value to be plotted in the senstivity analysis.
%   vec1max is the maximum vec1 value to be plotted in sensitivity analysis.
%   vec2max is the maximum vec2 value to be plotted in sensitivity analysis.
%   vec3max is the maximum vec3 value to be plotted in sensitivity analysis.
%   axistoken is a string designating whether or not the y-axis of all plots
%     should be equal.  Valid tokens are 'equal' and 'none'.
%   The output arguments plot1handles, plot2handles, and plot3handles are the
%     graphic handles to the first, second, and third rows of the sensitivity 
%     analysis subplots, respectively.

inc1 = (vec1max - vec1min) / 100;
inc2 = (vec2max - vec2min) / 100;
inc3 = (vec3max - vec3min) / 100;

analysisvec1 = vec1min:inc1:vec1max;
analysisvec2 = vec2min:inc2:vec2max;
analysisvec3 = vec3min:inc3:vec3max;

RMSkt = calcRMSkt(x+vec3best,z,thetadeg,a,analysisvec1,vec2best);
RMSb = calcRMSb(x+vec3best,z,thetadeg,a,vec1best,analysisvec2);
RMSoffset = calcRMSoffset(x,z,thetadeg,a,vec1best,vec2best,analysisvec3);

subplot(3,1,1)
plot1handles = plot(analysisvec1,RMSkt,'k-');
a1 = axis;
subplot(3,1,2)
plot2handles = plot(analysisvec2,RMSb,'k-');
a2 = axis;

subplot(3,1,3)
plot3handles = plot(analysisvec3,RMSoffset,'k-');
a3 = axis;

if(strcmp(axistoken,'equal'))
   minaxis = min([a1(3) a2(3) a3(3)]);
   maxaxis = max([a1(4) a2(4) a3(4)]);
   subplot(3,1,1)
   axis([a1(1) a1(2) minaxis maxaxis]);
   subplot(3,1,2)
   axis([a2(1) a2(2) minaxis maxaxis]);
   subplot(3,1,3)
   axis([a3(1) a3(2) minaxis maxaxis]);

end

   



%--
function RMS = calcRMSkt(x,z,thetadeg,a,kt,bdeg)

for(i=1:length(kt))
   
  Zpred = finitescarp(thetadeg,a,x,kt(i),bdeg);

  RMS(i) = sqrt(mean((Zpred - z).^2));
  
end

%--
function RMS = calcRMSb(x,z,thetadeg,a,kt,bdeg)

for(i=1:length(bdeg))
  Zpred = finitescarp(thetadeg,a,x,kt,bdeg(i));

  RMS(i) = sqrt(mean((Zpred - z).^2));
  
end

%--
function RMS = calcRMSoffset(x,z,thetadeg,a,kt,bdeg,offsetvec)

for(i=1:length(offsetvec))
  Zpred = finitescarp(thetadeg,a,x+offsetvec(i),kt,bdeg);

  RMS(i) = sqrt(mean((Zpred - z).^2));
  
end
