function [mps, x, y, cog] = ll2mps(dday,lon,lat);
% [mps, x, y, cog] = ll2mps(dday,lon,lat);
% output sizes match input sizes by padding the first speed with 0 
%     and assuming the speed should be applied to the last time of 
%     the interval (consistent with shipboard adcp data)
% works on Nx1 or 1XN vectors for 3 arguments; or Nx3 (or 3XN) if one argument
% 

mps = [];
x = [];
y = [];
cog = [];

if nargin == 1
   [rows,cols]=size(dday);
   txy = dday;
   if rows == 3
      txy = txy';
   end
   dday = txy(:,1);
   lon = txy(:,2);
   lat = txy(:,3);
else
   switchback = 0;
   [rows, cols] = size(dday + lon + lat);
   if (rows == 1)
      switchback = 1;
      dday = dday';
      lon = lon';
      lat = lat';
   end
end

if (length(lat) == 1)
   mps = NaN;
   x = 0;
   y = 0;
   cog = NaN;
   return
end


% warn about bad points
badi = find(isnan(lon+lat));
if (badi > 0)
   fprintf('ll2mps: filling %d bad positions\n',length(badi))
end

dy = lat_to_m(diff(lat),mid(lat));
dx = lon_to_m(diff(lon),mid(lat));

%% assume constant speed, not constant position for missing points
dt = 86400*diff(dday);
u = fillends(dx./dt);  
v = fillends(dy./dt);  

x = cumsum(u.*dt);  x = [0; x];
y = cumsum(v.*dt);  y = [0; y];

dx = diff(x);
dy = diff(y);


ds = sqrt(dx.^2 + dy.^2);
mps = ds./(diff(dday)*86400);
cog = [90-180*angle(dx+i*dy)/pi];
cog = [cog(1);cog];
% make it so cog is 0-360
cog = rem(cog+360,360);

mps = [mps(1); mps];   %replicate the first value (using zero makes no sense)

if (switchback)
   mps = mps';
   cog = cog';
   x = x';
   y = y';
end
