function [x,y,t] = roms_cruisetrack(lonpath,latpath,g,dkm,speed)
% [x,y,t] = roms_cruisetrack(lonpath,latpath,grd,dkm,speed)
% 
% Calculate a longitude/latitude trajectory analogous to a cruise track
%
% Inputs:
%    lonpath,latpah = a sequence of way points defining the track
%    grd = roms grid structure from roms_get_grid
%    dkm = nominal station interval in kilometers (default = 30)
%    speed = vessel speed for computing time (default = 15 knots)
%
% Outputs:
%    x,y = lon/lat coordinates of the sampling stations
%    t (optional) = sampling time in days with t(1) = 0
%
% John Wilkin

if nargin < 4
  % hi-res xbt sample interval is about 30 km
  dkm = 30;
end
if nargin < 5
  speed = 15;
end

if nargin < 3
  % load the roms grid structure
  g = grd('eac4');
end

% compute the distance along the trajectory of the way points
dspath  = sw_dist(latpath,lonpath,'km');
spath   = cumsum([0; dspath(:)]);
distance = spath(end);

% interpolate path to nominal sampling interval dkm
num = round(distance/dkm);
lontrk = interp1(spath,lonpath,0:(distance/num):distance);
lattrk = interp1(spath,latpath,0:(distance/num):distance);

% time along path in days is computed assuming a ship speed of 15 kn
% starting from zero
km_per_day = speed*111.12/60 * 24;
timtrk = (0:(distance/num):distance)/km_per_day;

% remove points that fall outside the model domain (defined by the psi
% points grid)
xu = [g.lon_psi(1,:)'; g.lon_psi(:,end); ...
  g.lon_psi(end,end:-1:1)'; g.lon_psi(end:-1:1,1)];
xv = [g.lat_psi(1,:)'; g.lat_psi(:,end); ...
  g.lat_psi(end,end:-1:1)'; g.lat_psi(end:-1:1,1)];
in = inpolygon(lontrk,lattrk,xu,xv);
lontrk = lontrk(in);
lattrk = lattrk(in);
timtrk = timtrk(in);

% remove points that fall on land
m = griddata(g.lon_rho,g.lat_rho,g.mask_rho_nan,lontrk,lattrk);
lontrk(find(isnan(m)==1)) = [];
lattrk(find(isnan(m)==1)) = [];
timtrk(find(isnan(m)==1)) = [];

x = lontrk;
y = lattrk;
if nargout > 2
  t = timtrk;
end

