function lonlat = intrpnav(nav, refU, time)
%
% interpolates missing fixes using the relative reference layer velocities
% (assuming that water velocity is constant over the interval)

nprfs = length(time);

lonlat = nav;
igood = find(~isnan(lonlat(1,:)));
ngood = length(igood);
if ngood == nprfs
  return;
end

%--- drop bad fixes at either end ---
lonlat = lonlat(:,igood(1):igood(ngood));

%--- bad values found, identify gaps to interpolate over ---
if (length(find(diff(igood) ~= 1)) ~= 0)
  df = diff(igood);
  df_1 = find(df ~= 1);
  gap = df(df_1)-1;
  ngap = length(gap);
  fprintf('%d navigation gaps to be filled by INTRPNAV\n', ngap);
  cumgap = cumsum(gap);
  ind0 = [df_1(1)  [(df_1(2:ngap) + cumgap(1:(ngap-1)))]];
  indn = ind0+gap+1;
  dt = diff(time)*86400;
  lonlat = [];
  i0 = igood(1);
  save intrp
  for ii = 1:length(gap)
    lonlat = [lonlat nav(:,i0:ind0(ii))];
    alat = 0.5*(nav(2,indn(ii)) + nav(2,ind0(ii)));
    dxmeter = lon_to_m((nav(1,indn(ii)) - nav(1,ind0(ii))),alat);
    dymeter = lat_to_m((nav(2,indn(ii)) - nav(2,ind0(ii))),alat);
    jj=[1:(gap(ii)+1)]';
    cumsumx = cumsum([refU(1,ind0(ii)+jj).*dt(ind0(ii)+jj-1)]);
    cumsumy = cumsum([refU(2,ind0(ii)+jj).*dt(ind0(ii)+jj-1)]);
    cumsumt = cumsum([dt(ind0(ii)+jj-1)]);
    dx = (dxmeter/cumsumt(gap(ii)+1) + cumsumx(gap(ii)+1)/cumsumt(gap(ii)+1) ...
         - cumsumx(1:gap(ii))./cumsumt(1:gap(ii))) .* cumsumt(1:gap(ii));
    dy = (dymeter/cumsumt(gap(ii)+1) + cumsumy(gap(ii)+1)/cumsumt(gap(ii)+1) ...
         - cumsumy(1:gap(ii))./cumsumt(1:gap(ii))) .* cumsumt(1:gap(ii));
    dlonlat = [m_to_lon(dx,alat);  m_to_lat(dy,alat)];
    [mlonlat nlonlat] = size(lonlat);
    lonlat = [lonlat lonlat(:,nlonlat)*ones(1,length(dlonlat(1,:)))+dlonlat];
    i0 = indn(ii);
  end
  lonlat = [lonlat nav(:,i0:igood(ngood))];
end
lonlat=[NaN*ones(2,(igood(1)-1)) lonlat NaN*ones(2,(length(nav(1,:))-igood(ngood)))];
