param_sel = 1;

if param_sel == 1
   param = data(:, water_temp);
   spread = [8 10 11 12 13 14 16];
elseif param_sel == 2
   param = data(:, water_cond);
   spread = [3.3 3.4 3.5 3.6 3.7 3.8 3.9];
elseif param_sel == 3
   param = data(:, m_fl_voltage);
   spread = [250 300 350 400 450 500 600 700 800];
elseif param_sel == 4
   param = data(:, m_bs_voltage);
   spread = [250 300 350 400 450 500 600 700 800];
end;
 
% Find non-NaN values, and block out first 45 seconds..
oklag = find((~isnan(d) & ~isnan(param)) & time > 45);

% Find derivative of water property sensor (assumed read at 4 Hz)
p_rate = param*NaN;
oktemp = find(~isnan(param));
p_rate(oktemp(2:length(oktemp))) = diff(param(oktemp))/0.25;
temp = conv(p_rate(oktemp), [1 1 1 1 1]/5);
p_rate(oktemp) = temp(3:(length(temp)-2));

% Establish variables used for computation below...
dpdt = p_rate(oklag);
dzdt = d_rate(oklag);
param = param(oklag);

% ldepth = fix_ctd_depth(d(oklag), data(oklag, m_pitch));
ldepth = data(oklag, m_depth) + 0.30 * sin(data(oklag, m_pitch));

ltime = data(oklag, m_present_time);

clear zlag, clear tlag; clear temp; clear oktemp;

for ilag = 1:max(yoyo),
   ilag,
   for isprd = 1:length(spread),
      
      % Find closes data point to desired parameter value...on down-going leg
      valid = find((yoyo(oklag) == ilag) & (d_rate(oklag) > 0.01)); 
      if ((spread(isprd) < max(param(valid))) & spread(isprd) > min(param(valid)))
         
         % Find the closest point to the desired parameter value
         adelta = abs(param(valid)-spread(isprd));
         minp = min(adelta);
         
         center = find(adelta - minp == 0);
         center = center(1);
         
         % Make sure we don't try to access beyond the end of the valid data region
         if (center-1 > 0) & (center+1 < length(valid)+1)         
            oklag2 = valid((center -1):(center+1));
            
            % Make sure selected data set brackets desired parameter value
            if (param(oklag2(1))-spread(isprd))*(param(oklag2(3))-spread(isprd)) < 0
               
               % Plot data set if diagnostic on...
               if diag_plot == 1,
                  hold off;
                  plot(ldepth(valid), param(valid));
                  hold on
                  plot(ldepth(oklag2), param(oklag2), 'o');
                  hold on;
               end;
               
               zlag(ilag*2-1, isprd) = interp1(param(oklag2), ldepth(oklag2), spread(isprd));
               tlag(ilag*2-1, isprd) = interp1(param(oklag2), ltime(oklag2), spread(isprd));
               dzdtlag(ilag*2-1, isprd) = dzdt(valid(center));
               dpdtlag(ilag*2-1, isprd) = dpdt(valid(center));
               dpdzlag(ilag*2-1, isprd) = dpdt(valid(center))/dzdt(valid(center));
            end;
         end;
      end;
      
      % Find closes data point to desired parameter value...on down-going leg
      valid = find((yoyo(oklag) == ilag) & (d_rate(oklag) < -0.01)); 
      if (spread(isprd) < max(param(valid)) & spread(isprd) > min(param(valid)))
         
         % Find the closest point to the desired parameter value
         adelta = abs(param(valid)-spread(isprd));
         minp = min(adelta);
         
         center = find(adelta - minp == 0);
         center = center(1);
         
         % Make sure we don't try to access beyond the end of the valid data region
         if (center-1 > 0) & (center+1 < length(valid)+1)         
            oklag2 = valid((center -1):(center+1));
            
            % Make sure selected data set brackets desired parameter value
            if (param(oklag2(1))-spread(isprd))*(param(oklag2(3))-spread(isprd)) < 0
               
               % Plot data set if diagnostic on...
               if diag_plot == 1,
                  hold off;
                  plot(ldepth(valid), param(valid));
                  hold on
                  plot(ldepth(oklag2), param(oklag2), 'o');
                  hold on;
               end;
               
               zlag(ilag*2, isprd) = interp1(param(oklag2), ldepth(oklag2), spread(isprd));
               tlag(ilag*2, isprd) = interp1(param(oklag2), ltime(oklag2), spread(isprd));
               dzdtlag(ilag*2, isprd) = dzdt(valid(center));
               dpdtlag(ilag*2, isprd) = dpdt(valid(center));
               dpdzlag(ilag*2, isprd) = dpdt(valid(center))/dzdt(valid(center));
            end;
         end;
      end;
   end;
end;

tlag(find(tlag == 0)) = NaN;
zlag(find(zlag == 0)) = NaN;

figure;
n = length(zlag);
plot(tlag(1:2:n, :), zlag(1:2:n, :), '-x'); 
hold on;
plot(tlag(2:2:n, :), zlag(2:2:n, :), '-o'); 
set( gca, 'YDir', 'rev' );
xlabel('Mission time (s)');
ylabel('Isotemp depth (m)');
hold off;
zoom on

% Compute errors...

bnd_deriv = (zlag(1:2:(n-3), :) - zlag(3:2:(n-1), :)) ./ (tlag(1:2:(n-3), :) - tlag(3:2:(n-1), :));
z_up_interp = (tlag(2:2:(n-2),:) - tlag(1:2:(n-3), :)) .* bnd_deriv + zlag(1:2:(n-3), :);

z_err = zlag(2:2:(n-2),:) - z_up_interp;

% clear ldepth, clear ltime, clear param, clear spread, clear minp, clear center
% clear adelta, clear valid

% Compute hysterysis

test = diff(zlag); 
ans = sum(test .* (-1) .^ cumsum(ones(size(test))))/length(test)