function [f, g] = refavg(u, ilayer, n_iterations)
%
% function [vprofile, trefl] = refavg(u, ilayer, n_iterations)
%
% Improved single-ping profile averaging scheme, using a
% reference layer, with reduced sensitivity to data gaps within
% that reference layer.  Therefore the reference layer can be
% larger than one would choose for the simplest ref layer
% averaging method, such as was used in DAS 2.48.
%
% arguments:
%  u (Nbins x Ntimes) : measured vel relative to platform
%  ilayer (1 x Nref)  : indices of ref layer (eg. 5:10 or 1:30)
%  n_interations (small integer > 1, say 2-3)
%
% output:
%  vprofile Nbins x 1  :  time-averaged profile relative to ref layer
%  trefl    1 x Ntimes :  reference layer time series
%
% The estimated velocity profile relative to the ship, for each
% ping, and assuming there is no change in ocean velocity during
% the time period, is:
%   u_rel = repmat(vprofile, 1, Ntimes) + repmat(trefl, Nbins, 1);
% so the residuals are, of course
%   u_residual = u - u_rel;
%
% EF 2000, on board Kaiyo, KY0006.
% bug fix and doc by JH
% EF 2001/12/31, augmented doc.
%

[nd, nt] = size(u);

uu = u;
f = zeros(nd, 1);
g = zeros(1,nt);


for it = 1:n_iterations
   [ff, gg] = ref_est(uu, ilayer);
   %fprintf('%f %f %f %f\n', mean(ff), std(ff), mean(gg), std(gg))
   uu = uu - ff(:, ones(1, nt)) - gg(ones(1, nd), :);
   f = f + ff;
   g = g + gg;
end


if (~isreal(f))
   nani = find(isnan(f));
   f(nani) = NaN + i*NaN;

   nani = find(isnan(g));
   g(nani) = NaN + i*NaN;
end


%---------------------------------------------------------------

function [f, g] = ref_est(u, ilayer);


nl = length(ilayer);
[nd, nt] = size(u);

[f, sf, nf] = mstdgap(u, 2);
f = f.' - mstdgap(f(ilayer));    %was mean(ilayer)

[g, sg, ng] = mstdgap(u(ilayer, :) - f(ilayer, ones(1, nt)));

%% This is slow; it can be sped up by using a simplified
%% local version of mstdgap, and probably eliminating the
%% function calls.  The gaps and numbers of points stay the
%% same for each iteration--no need to keep recalculating,
%% (Asking for the medians in mstdgap makes it much slower.)



