function [xi, yi, zi, vi] = grid3d(X, Y, Z, V, xInc, yInc, zInc)
% GRID3D   - Grids FSI data onto an evenly spaced 3-D mesh
%
% Use as: [XI, YI, ZI, VI] = grid3d(X, Y, Z, V, xInc, yInc, zInc)
%
% Inputs: X,Y,Z,V = Vectors of 3-d data. X,Y, and Z specify
%         the spatial coordinates. V contains the data values.
%         xInc, yInc, zInc = Scalar specifying the distance 
%         between nodes on the nesh grid.
%
% Output: XI,YI,ZI,VI = 3d arrays of gridded data. XI,YI and ZI
%         specify the spatial coordinates while VI contains the data values
%
% Note: VI is determined by binning the data within a each z level while
%       looping through each z-level.

% Brian Schlining
% 25 Jul 2000

fprintf(1, 'Gridding data ...');

% Calculate the 3-d grid points
xMin = min(X);
xMax = max(X);
yMin = min(Y);
yMax = max(Y);
zMin = min(Z);
zMax = max(Z);
[xi, yi, zi] = meshgrid([xMin:xInc:xMax], [yMin:yInc:yMax], [zMin:zInc:zMax]);
vi = ones(size(xi))* NaN;

nMax = length([zMin:zInc:zMax]);
n    = 0;
bad = [];

for zn = min(Z):zInc:max(Z)                   % Loop through all depths by 2m increments
   fprintf(1, '.');
   n = n + 1;
   good = find((Z >= zn) & (Z < (zn + zInc)));  % Find all points in this 2m bin 
   try 
      vi(:,:,n) = griddata(X(good), Y(good), V(good), xi(:,:,n), yi(:,:,n));
   catch
      vi(:,:,n) = NaN;
      bad = [bad; n];                         % Record the bad layers to be fixed later
   end
end

m = size(vi, 3);              % 3rd dimension length of vi
inc = 1;                      % Pointer 
for i = 1:length(bad)         % Loop through each bad layer and try to fix it
   n = bad(i);                % The current bad layer to be fixed
   if n == 1 | n == m
      vi(:, :, n) = NaN;
      warning('NaN''s present in vi')
   else
      while ismember(n+inc, bad)
         inc = inc + 1;
         if n + inc > m
            vi(:, :, n:m) = NaN;
            warning(['vi(:,:,' num2str(n) ') through vi(:,:,' num2str(m) ') are NaN''s'])
            return
         end
      end
      
      if inc > 1
         n = [n:n+inc];
      end
      
      vi(:,:,n) = interp3(xi(:,:,[n(1)-1 n(1)+inc]), ...
         yi(:,:,[n(1)-1 n(1)+inc]), ...
         zi(:,:,[n(1)-1 n(1)+inc]), ...
         vi(:,:,[n(1)-1 n(1)+inc]), ...
         xi(:,:,[n]), ...
         yi(:,:,[n]), ...
         zi(:,:,[n]));
      
      inc = 1;
   end
   
end

zi = zi + zInc/2;


