% cCTD_yoyo computes the yoyo info...
%
%      ctd_time:  	uniformly sampled time, 0.2s period
%      ctd_temp:   	water temperature sampled on time base
%      ctd_pres:  	pressure (bar) sampled on time base
%      ctd_depth:		depth (m) sampled on time base
%      ctd_salt:   	salinity sampled on time base
%      ctd_dens: 		density sampled on time base
%
% Yoyos require depth rate > yoyo_discriminator m/s (typically 0.1 m/s)
%
% cCTD_yoyo assumes:
%    fix_depth.m has been run
%    cCTD_salt.m has been run
%  
% J. G. Bellingham
% April 4, 1999

yoyo_discriminator = 0.1;

% Use depth rate computed by fix_depth.m
ok = find(~isnan(d_rate));
ctd_d_rate = interp1q(data(ok,m_present_time), d_rate(ok), ctd_time);

% Classify yoyos...
dwn = find(diff(ctd_d_rate > yoyo_discriminator) == 1);

% Make sure that yoyos are long enough to really be full yoyos
okok = find(diff(dwn) > 100);
okok(2:(length(okok)+1)) = okok+1;
okok(1) = 1;
dwn = dwn(okok);
 
% Make an array with yoyo numbers...
ctd_yoyo = NaN * ctd_time;
if length(dwn) > 2
   ctd_yoyo(1:dwn(1)) = 1;
   for i = 2:length(dwn)
      ctd_yoyo(dwn(i-1):dwn(i)) = i;
   end;
   ctd_yoyo(dwn(length(dwn)):length(ctd_time)) = i+1;
end;

% Classify going up or down
ctd_going_up = find(ctd_d_rate < -yoyo_discriminator);
ctd_going_down = find(ctd_d_rate > yoyo_discriminator);

% -------------------------------------------------------------------
% Do upcast/downcast processing...

ctd_bfrq = NaN * ctd_time;
ctd_vort = NaN * ctd_time;
ctd_p_ave = NaN * ctd_time;
ct = smooth(ctd_temp, 10);
cs = smooth(ctd_salt, 10);
cp = smooth(ctd_pres, 10);

for i = 1:max(yoyo),
   seg1 = find((ctd_yoyo == i) & (ctd_d_rate < -yoyo_discriminator));
   seg2 = find((ctd_yoyo == i) & (ctd_d_rate > yoyo_discriminator));
   
   if length(seg1) > 1
      [a,b,c] = ...
         sw_bfrq(cs(seg1),ct(seg1), cp(seg1)*10,45);  % Default latitude is Boston....
      
      ctd_bfrq(seg1(2:(length(seg1)-1))) = (a(2:length(a)) + a(1:(length(a)-1)))/2;
      ctd_vort(seg1(2:(length(seg1)-1))) = (b(2:length(b)) + b(1:(length(b)-1)))/2;
      ctd_p_ave(seg1(2:(length(seg1)-1))) = (c(2:length(c)) + c(1:(length(c)-1)))/2;
   end;
   
   if length(seg2) > 1
      [a,b,c] = ...
         sw_bfrq(cs(seg2),ct(seg2), cp(seg2)*10,45);  % Default latitude is Boston....
      
      ctd_bfrq(seg2(2:(length(seg2)-1))) = (a(2:length(a)) + a(1:(length(a)-1)))/2;
      ctd_vort(seg2(2:(length(seg2)-1))) = (b(2:length(b)) + b(1:(length(b)-1)))/2;
      ctd_p_ave(seg2(2:(length(seg2)-1))) = (c(2:length(c)) + c(1:(length(c)-1)))/2;
   end;
   
   % [ctd_bfrq(seg2), ctd_vort(seg2), ctd_p_ave(seg2)] = ...
   % sw_bfrq(ctd_salt(seg2),ctd_temp(seg2), ctd_pres(seg2)*10,45);  % Default latitude is Boston....
end;

clear cyc; clear tens; clear offset; clear d_off; clear bad; clear d_r;  
clear temp; clear ok; clear dwn; clear tmplt; clear ok1; 
clear seg1; clear seg2;

