function [f, g, resid] = refavg(u, ilayer, n_iterations)
%
% function [vprofile, trefl, resid] = 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)
%                     : defaults to all and forces default of n_iterations)
%  n_interations (small integer > 1, say 2-3)
%                     : defaults to 2
%
% 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.
%
% EF 2003/12/24 added ufake function, rescued from obsolete
%    directory, for future testing.
%

[nd, nt] = size(u);


if nargin < 3
   n_iterations = 2;
end
if nargin < 2
   ilayer = 1:nd;
end

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

%%
fM = repmat(f, 1, nt);
gM = repmat(g, nd, 1);
resid = u - (fM + gM);


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

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.)

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

function [u, dd, tt] = ufake(dims, ampt, ampd, amptr, ampdr)
% function [u, f, g] = ufake(dims, ampt, ampd, amptr, ampdr)
%
%  Generate a set of u profiles that perfectly match
%  the form assumed by refavg: that is, a single function
%  of z (f; down a column) times a single function of time
%  (g; increasing with column number)
%
%     dims  2-d array with size of u matrix
%     ampt  magnitude of the non-random part of the time dependence;
%               it is a step function, ampt in the middle third,
%               zero elsewhere.
%     ampd  magnitude of the non-random part of the depth function, a
%               half-cosine
%     amptr random noise amplitude in the time function
%     ampdr                               depth function
%
%  To simulate the case where the depth range is reduced when
%  the ship is underway (middle third), the bottom half of the
%  profile is NaN during the middle third.
%

nr = dims(1); nc = dims(2);
u = zeros(dims);

tt = zeros(1, nc);
ii = round(nc/3):round(2*nc/3);
tt(ii) = ampt;
tt = tt + randn(1, nc) * amptr;



dd = zeros(nr, 1);
d = (1:nr).'/nr;
dd(:) = ampd * cos(d*pi) + ampdr * randn(nr, 1);

u = dd(:,ones(1,nc)) +   tt(ones(1,nr),:);

u((round(nr/2):nr), ii) = NaN;




