
% lowpassdrifters.m
%
% lowpass the drifter velocities and
% NCEP daily wind stresses at 
% a period of dT,.
% Also, NaN-out the speeds where elon or elat
% are greater than 0.03 degrees 
% (approx. 10 cm/s max error for speed).

dT=5;	% lowpass filter width (days)
% make filter with "butter": for dT=2,
% the amplitude of a filered sine wave with
% unfiltered amplitude of 1 becomes:
%   Period (d)     Amplitude
%       4           .94
%       3           .84
%       2           .50
%       1           .06
[bf,af]=butter(2,2*median(diff(t))/dT);
fprintf('Lowpassing data at %d days.\n',dT);

q=find(diff(id) | diff(t)>t(2)-t(1));
q=[0;q;length(id)];

for i=2:length(q);
  fprintf(' %d of %d\n',i-1,length(q)-1);
  a=q(i-1)+1;b=q(i);
  if (b>a+5); % data must be 3 times filter length
      Tx(a:b)=filtfilt(bf,af,Tx(a:b));
      Ty(a:b)=filtfilt(bf,af,Ty(a:b));
      uu=fillnans(u(a:b),0);
      vv=fillnans(v(a:b),0);
      uu=filtfilt(bf,af,uu);
      vv=filtfilt(bf,af,vv);
      qj=find(isnan(u(a:b)));
      uu(qj)=NaN;
      qj=find(isnan(v(a:b)));
      vv(qj)=NaN;
      u(a:b)=uu;
      v(a:b)=vv;
  else % data too short to lowpass
      u(a:b)=NaN;
      v(a:b)=NaN;
  end;  % (if b>a+5)
 
end;  % loop over i

% NaN-out velocities where elon or elat
% is greater than 0.03 degree
q=find(elon>.03 | elat>.03);
u(q)=NaN; v(q)=NaN;

% NaN-out sst where esst > 0.1 deg.
%q=find(esst>0.1);
%sst(q)=NaN;

lowpass=1;
clear a b af bf wt dT q i j ang dist lonj qj uu vv

