function [y1, yblock] = ave_fit(x, y_av, x1)

%
% function [y1, yblock] = ave_fit(x, y_av, x1)
%
% Fit a piece-wise parabolic curve to a function for
% which point-to-point averages are known.  For example,
% Given the average geostrophic velocities between CTD
% stations, generate a smooth curve with the same
% average value between stations.
%
% Inputs:
%   x      : vector of length n; boundaries (CTD stations)
%   y_av   : vector of length (n-1); average between boundaries
%   x1     : new grid on which smooth function will be evaluated
%
% Outputs:
%   y1     : piece-wise parabolic fit at x1 grid
%   yblock : piece-wise constant fit at x1 grid
%
% If invoked with no input parameters, ave_fit will run a simple
% demo/test routine.
%
% Notes: at interior nodes, the function and its first derivative
%        are continuous.  The first derivatives are zero at the
%        end points.
%
%        This should be extended to handle matrix y_av,
%        as required for contouring geostrophic currents.
%        It may need more error checking.
%        Upon generalization, this might be renamed paraspline,
%        or something similar.

% Shuiming Chen, around 1997
% EF 2002/01/12
%    Edited, largely recoded; added modification of demo by Kathy Donohue.

if nargin == 0,
   ave_fit_demo;
   help ave_fit
   return
end

if nargin ~= 3,
   help ave_fit
   error('Three input arguments are required.')
end

if ~isempty(find(isnan(x1)))
   help ave_fit
   error('The third argument, output grid (x1), contains NaNs.')
end

nx1 = length(x1);
nx = length(x);
nsegs = length(y_av);  % number of segments
if nsegs ~= (nx - 1),
   help ave_fit
   error('The length of y_av must be one less than the length of x.')
end

ncoefs = 3 * nsegs;    % number of coefficients: 3 per segment

rr = zeros(ncoefs,1);  % RHS of linear equation.
rr(1:3:end) = y_av(:);

cc = zeros(ncoefs);    % Model matrix.

x = x(:);              % Ensure it is a column vector.
dx = diff(x);
dx_i = dx(1:(end-1));  % Inner values of dx.
% Coefficients for the average constraint:
idiag = 1 + (0:(ncoefs-1)) * (ncoefs + 1);
isub = idiag(1:3:ncoefs);
cc(isub) = dx.^2 / 3;            % diag
cc(isub + ncoefs) = dx / 2;   % one column right
cc(isub + ncoefs * 2) = 1;    % two columns right

% Coefficients for the continuity constrant:
isub = idiag(2:3:(ncoefs-3));
cc(isub - ncoefs) = dx_i.^2;
cc(isub)          = dx_i;
cc(isub + ncoefs) = 1;
cc(isub + 4 * ncoefs) = -1;

% Coefficients for the continuous deriv:
isub = idiag(3:3:(ncoefs-3));
cc(isub - 2 * ncoefs) = 2 * dx_i;
cc(isub - ncoefs)     = 1;
cc(isub + 2 * ncoefs) = -1;

% Boundary conditions: zero deriv at start and end of sequence.
cc((ncoefs-1), 2) = 1;
cc(ncoefs, (ncoefs - [2 1])) = [ 2 * dx(end), 1];

ss = cc\rr;

aa = ss(1:3:end);    % quadratic
bb = ss(2:3:end);    % linear
cc1= ss(3:3:end);    % constant

y1 = NaN * zeros(size(x1));
yblock = y1;

% We use inclusive intervals (>= and <= on both ends) so we can
% be sure to cover all the points with a single operation; a
% point in x1 that coincides with a point in x will be calculated
% twice, but this adds nearly nothing to the computation time.
for ii = 1:nsegs
   mseg = (x1 >= x(ii) & x1 <= x(ii+1)) | (x1 <= x(ii) & x1 >= x(ii+1));
   if (sum(mseg) > 0)
      y1(mseg) = aa(ii) * (x1(mseg)-x(ii)).^2 + bb(ii) * (x1(mseg) - x(ii)) + cc1(ii);
      yblock(mseg) = y_av(ii);
   end
end


%===================================================================

function ave_fit_demo

% An x-grid for plotting and for defining a pseudo-data function;
% think of it as being in kilometers.
x = 1:0.05:200;
nx = length(x);

% Pseudo-data; looks a bit like cross-track velocity component
% from an ADCP section across an eddy field:
u = cos(2*pi/20.*x) + sin(2*pi/12.*x) + cos(2*pi/80.*x) + ...
    0.1 * randn(size(x));

% Here are the x-locations of CTD stations:
xx = [x(1):5:mean(x) mean(x)+10:10:x(end-1) x(end)];
nxx = length(xx);

% uavg will be the average of u between the stations (xx points);
% it is what one would calculate from geostrophy if one knew
% absolute dynamic height at the stations.
uavg = NaN .* ones(1, nxx-1);
uavg_block = NaN .* zeros(nx, 1);
for ii = 1:(nxx - 1),
   iwant = find(x>= xx(ii) & x < xx(ii+1) );
   uavg(ii) = mean(u(iwant));
end

[uu, uavg_block] = ave_fit(xx, uavg, x);

plot(x, u, x, uavg_block, x, uu, xx, 0.*xx, '.');
legend('U', 'Avg U', 'Fit to avg U', 'Stations');
title('Built-in test of ave\_fit.m')




