function [velocityTotal, velocityLatitude, velocityLongitude] = horizontalVelocity(pLatitude, pLongitude, pInterval, pMethod, pNop)

% 29 October 2003, HHFurey.   This routine does not account for
% discontinuous segments of track being passed, and therefore returns
% erroneous results.  Not true - data passed in is bad - some segments that
% are discontinuous are passed through with NaNs as filler, but others are not!
% Hmmm ... working on r023 as a test case, this has three discontinous segments
% of track.


if nargin < 4
  pMethod='average';
end

if nargin < 5
  pNop=1;
end
%save debug_02
%                       nautic miles in kilometers
nm_in_km=60*1.852;

velocityTotal=-9.99*ones(size(pLatitude))*NaN;
velocityLatitude=-9.99*ones(size(pLatitude))*NaN;
velocityLongitude=-9.99*ones(size(pLatitude))*NaN;

if strcmp(pMethod,'average');    % centered mean between point ix+nop and ix-nop, centered at ix
  for ix=1:length(pLatitude)
    if ~((ix-pNop < 1) | (ix+pNop > length(pLatitude)))
      if ~any(isnan(pLatitude([ix-pNop:ix+pNop])))
        dkm = artoa.data.calculateGeodist([pLatitude(ix-pNop) pLongitude(ix-pNop)] , [pLatitude(ix+pNop) pLongitude(ix+pNop)]);   %[m/s]
        velocityTotal(ix)=(dkm*1000)/(pInterval*2*pNop);
        velocityLatitude(ix)=((pLatitude(ix+pNop)-pLatitude(ix-pNop))*nm_in_km*1000)/(pInterval*2*pNop);
        velocityLongitude(ix)=((pLongitude(ix+pNop)-pLongitude(ix-pNop))*nm_in_km*1000*cos(pLatitude(ix)/180*pi))/(pInterval*2*pNop);
      end
    end
  end
elseif strcmp(pMethod,'backward');  % shifted means between points ix and ix-nop
  for ix=1:length(pLatitude)
    if ~(ix-pNop < 1)
      if ~any(isnan(pLatitude([ix-pNop,ix])))
        dkm=artoa.data.calculateGeodist([pLatitude(ix-pNop) pLongitude(ix-pNop)], [pLatitude(ix) pLongitude(ix)]);
        velocityTotal(ix)=(dkm*1000)/(pInterval*pNop);
        velocityLatitude(ix)=((pLatitude(ix)-pLatitude(ix-pNop))*nm_in_km*1000)/(pInterval*pNop);
        velocityLongitude(ix)=((pLongitude(ix)-pLongitude(ix-pNop))*nm_in_km*1000*cos(pLatitude(ix)/180*pi))/(pInterval*pNop);
      end
    end
  end
elseif strcmp(pMethod,'forward');   % shifted means between points ix and ix+nop
   for ix=1:length(pLatitude)
    if ~(ix+pNop > length(pLatitude))
      if ~any(isnan(pLatitude([ix,ix+pNop])))
        dkm=artoa.data.calculateGeodist([pLatitude(ix+pNop) pLongitude(ix+pNop)], [pLatitude(ix) pLongitude(ix)]);
        velocityTotal(ix)=(dkm*1000)/(pInterval*pNop);
        velocityLatitude(ix)=((pLatitude(ix+pNop)-pLatitude(ix))*nm_in_km*1000)/(pInterval*pNop);
        velocityLongitude(ix)=((pLongitude(ix+pNop)-pLongitude(ix))*nm_in_km*1000*cos(pLatitude(ix)/180*pi))/(pInterval*pNop);
      end
    end
 end
 elseif strcmp(pMethod,'spline');  % spline at 2hour intervals
	%  spline latitude and longitudes and calculate velocities
   %  dummies must be excluded for the spline
   resolution = 2;                   % this is the resolution at which spline will be performed
   step = pInterval/3600/24;          % interval lenght between successive data points [days]
   time = 0:step:length(pLatitude)*step;   % creates virtural time basis for raw data
   time = time(:);
   index = find(pLatitude<999);            % find dummies, might be different value within artoa
	pLatitude   = pLatitude(index); pLongitude=pLongitude(index); time = time(index);  % filter out dummies

	% new timebase: every 2 hours including first and last time
    if isempty(time)
        dummy = 3;
    end
	time_min = time(1);
   time_max = time(size(time,1)); % save for later
   %time_spl  = linspace(time_min, time_max,(time_max-time_min)*24/resolution+1 );  %  time basis for spline [days]
   %Corrected by T. Reynaud: 22/06/2005
   time_spl = time_min:resolution/24:time_max;
   time_spl  = time_spl(:);
   
	lat_spl_f = csape(time,pLatitude,[2 1]); 
	lat_spl   = fnval(lat_spl_f,time_spl);
	v_spl_f   = fnder(lat_spl_f);
	v_spl     = fnval(v_spl_f,time_spl) *60*1.852*1000/(3600*24) ; % [m/s]  
   
	lon_spl_f = csape(time,pLongitude,[2 1]); 
	lon_spl   = fnval(lon_spl_f,time_spl);
	u_spl_f   = fnder(lon_spl_f);
	u_spl     = fnval(u_spl_f,time_spl) .*cos(lat_spl/180*pi)*60*1.852*1000/(3600*24) ; % [m/s]  
   
   % subsample to original rate:
   %save debug_tr
   for i = 1:length(time)
%      indi(i) = find(time_spl == time(i));
%Modified by T.Reynaud: 22/03/2005
      indi(i) = find(abs(time_spl-time(i)) <= 0.000000000001);
   end
   
   %figure
   %title('control plot for spline function')
   %plot(lon_spl,lat_spl); hold on
   %plot(lon,lat,'rx');
   %zoom on;
   
   velocityLongitude(index) = u_spl(indi);
   velocityLatitude(index) = v_spl(indi);
   velocityTotal(index) = sqrt(u_spl(indi).^2 + v_spl(indi).^2);
   
end


%%%%%%%%%%%%% idea for faster calculation of 'cetered mean'

%v = NaN*ones(size(data(:,7)))
%v(2:length(data)-1,1) = data(1:length(data)-2,6)
%v(2:length(data)-1,2) = data(3:length(data),6)
%v = diff(v')'
%v = v*60*1.852*1000*100/3600/24