function [u,v,spd,head] = llt2spd(lon,lat,time)

% Synthax : [u,v,spd,head] = llt2spd(lon,lat,time)
% 
% Returns the speed 'spd' and it's eastward/northward components 'u' and
% 'v' in (m/s), and the heading 'head' in degrees with 0 as east increasing
% towards the north (90 degrees). 
%
% Values are calculated from a latitude/longitude time series with 
% coordinates in decimal degrees and time in days.
% 
% Values are estimated between longitude/latitude coordinates and then
% interpolated linearly to the input 'time' vector. Begining and end
% values are extrapolated using the "pchip" method.

%% make everything column vectors
if ~iscolumn(lon)  ; lon  = lon'  ; end
if ~iscolumn(lat)  ; lat  = lat'  ; end
if ~iscolumn(time) ; time = time' ; end

%% make differential vectors
dlon = diff(lon) ;
dlat = diff(lat) ;
dt   = diff(time) ;
dl   = m_lldist(lon,lat)*1000 ; % (m)

%% calculate norm of speed
spd = dl./(dt*24*60*60) ;

%% calculate heading
head = range_pm180_2_360(atan2d(dlat,dlon)) ;

%% decompose into eastward/northward
u = spd.*cosd(head) ;
v = spd.*sind(head) ;

%% make the time vector
t = time(1:end-1)+0.5*dt(1);

%% interpolate to the input time grid
u    = interp1(t,u,time,'linear','extrap') ;
v    = interp1(t,v,time,'linear','extrap') ;
spd  = interp1(t,spd,time) ;
spd(1)   = sqrt(u(1).^2+v(1).^2);
spd(end) = sqrt(u(end).^2+v(end).^2);
head     = interp1(t,head,time) ;


end