function [ di wi ti ni ] = depth_of_parameter( varargin )

% depth_of_parameter finds the depth at which param_in.value equals
% param_target.
%
% The routine assumes that you have run yoyo_classifier, and expects the
% outputs of that function as inputs here.
%
% Input arguments:
% (param_in, param_target, depth, time, yoyo_number, whichway, delta_time, start_depth)
%
%   param_in = parameter structure with .value and .time to be analyzed
%   param_target = parameter value for which one wishes to the know the depth
%   depth, time, yoyo_number, and whichway are outputs of yoyo_classifier
%   delta_time = optional shift in time
%

% Get arguments
param_in = varargin{1};
param_target = double(varargin{2});

depth = double(varargin{3});
time = double(varargin{4});
yoyo_number = (varargin{5});
whichway = double(varargin{6});

if nargin == 7
    delta_time = varargin{7};
else
    delta_time = 0;
end;

%  Make sure time array of parameter is single valued!  Amazing I have to check for this...
ok = [1 find(diff(param_in.time) > 0)+1];

%  Interpolate parameter to depth time base
param = interp1(param_in.time(ok), double(param_in.value(ok)), time+delta_time);
n_param = length(param);

% Generic search for shallowest depth at which parameter has desired value
% for each yoyo leg.
n_yoyos = yoyo_number(length(yoyo_number));
zero_crossing = param * 0;
zero_crossing((param(1:n_param-1)-param_target).*(param(2:n_param)-param_target) <= 0) = 1;

for i = 1:n_yoyos
    zi = find(yoyo_number == i & zero_crossing == 1);
    ni(i) = length(zi);
    
    if ni(i) >= 1
        % Get the zero crossing at the shallowest depth...
        if whichway == 1
            shallow = 1;
        else
            shallow = ni(i);
        end;
        
        % 
        ok = bound(zi(shallow)-1:zi(shallow)+2, 1, n_param);
        p = polyfit(param(ok)-param_target, time(ok), 1);
        
        ti(i) = p(2);
        di(i) = interp1(time(ok), double(depth(ok)), ti(i));
        wi(i) = whichway(zi(shallow));
        ii(i) = i;
    else
        di(i) = NaN;
        wi(i) = NaN;
        ii(i) = i;
    end;
end;

std(diff(di(~isnan(di))))

hold off;
plot(ii(wi > 0.5), di(wi > 0.5), '.r');
hold on;
plot(ii(wi < -0.5), di(wi < -0.5), '.b');
plot(ii, di, 'g');
hold off;
