function [ix,iy,iz]=xyzindex(varargin)
% [ix]=xy2index(x,xGr)
% [ix,iy]=xyzindex(x,y,xGr,yGr)
% [ix,iy,iz]=xyzindex(x,y,z,xGr,yGr,zGr)
% xGr,yGr,zGr are mesh grid coordinates, Z may be a complete Z grid
% of size [Nx,Ny,Nz+1]
% x,y,z are point coordinate
% computes the cell indices of points in 1D, 2D or 3D grid
%
% TO 090317 090509 100510


% Copyright 2009 Theo Olsthoorn, TU-Delft and Waternet, without any warranty
% under free software foundation GNU license version 3 or later

switch length(varargin)
    case 2,
        x=varargin{1}; xGr=varargin{2}; Nx=length(xGr)-1;
        ix=NaN(size(x));
        I=find(x>=min(xGr) & x<=max(xGr));  % these are inside the grid
        sx=sign(xGr(end)-xGr(1));           % ascending or descending grid
        for i=1:length(I);
            ii=I(i);
            if sx>0, ix(ii)=find(xGr<=x(ii),1,'last'); else ix=find(x(ii)<=xGr,'last',1); end
        end
        ix(I)=min(Nx,ix(I));
    case 4
        x=varargin{1}; xGr=varargin{3}; Nx=length(xGr)-1;
        y=varargin{2}; yGr=varargin{4}; Ny=length(yGr)-1;
        if ~all(size(x)==size(y)), error('x and y must be of same size'); end
        ix=NaN(size(x));
        iy=NaN(size(y));
        I=find(x>=min(xGr) & x<=max(xGr) & y>=min(yGr) & y<=max(yGr));
        sx=sign(xGr(end)-xGr(1));
        sy=sign(yGr(end)-yGr(1));
        for i=1:length(I);
            ii=I(i);
            if sx>0, ix(ii)=find(xGr<=x(ii),1,'last'); else ix(ii)=find(x(ii)<=xGr,1,'last'); end
            if sy>0, iy(ii)=find(yGr<=y(ii),1,'last'); else iy(ii)=find(y(ii)<=yGr,1,'last'); end
        end
        ix(I)=min(Nx,ix(I)); iy(I)=min(Ny,iy(I));
    case 6
        x=varargin{1}; xGr=varargin{4}; Nx=length(xGr)-1;
        y=varargin{2}; yGr=varargin{5}; Ny=length(yGr)-1;
        z=varargin{3}; zGr=varargin{6}; Nz=length(zGr)-1;
        if ~all(size(x)==size(y))
            error('x, y and z must be of same size');
        elseif length(x(:))==length(z(:)),
            is3DZ=0; % zGr as a vector in arbitrary direction
            sz=zGr(end)-zGr(1);
        else
            sizeZ=size(zGr);
            if length(sizeZ)==3 && sizeZ(1)==Ny+1 && sizeZ(2)==Nx+1
                is3DZ=1; % full 3D matrix given
                sz=zGr(1,1,end)-zGr(1,1,1);
            else
                error('size of zGr does not match those of xGr and yGr');
            end
        end
        ix=NaN(size(x));
        iy=NaN(size(y));
        iz=NaN(size(z));
        I=find(x>=min(xGr) & x<=max(xGr) & y>=min(yGr) & y<=max(yGr));
        sx=sign(xGr(end)-xGr(1));
        sy=sign(yGr(end)-yGr(1));
        for i=1:length(I)
            ii=I(i);
            if sx>0, ix(ii)=find(xGr<=x(ii),1,'last'); else ix(ii)=find(x(ii)<=xGr,1,'last'); end
            if sy>0, iy(ii)=find(yGr<=x(ii),1,'last'); else iy(ii)=find(y(ii)<=yGr,1,'last'); end
        end
        ix(I)=min(Nx,ix(I)); iy(I)=min(Ny,iy(I));
        
        % Now look after the z
        I=find(~isnan(ix));  % only points where x and y are inside grid
        for i=1:length(I)
            ii=I(i);
            if ~is3DZ
                if sz>0, iz(ii)=find(zGr(:)<=z(ii),1,'last'); else iz(ii)=find(z(ii)<=zGr(:),1,'last'); end
            else
                if sz>0, iz(ii)=find(squeeze(zGr(iy(ii),ix(ii),:))<=z(ii),1,'last');
                else     iz(ii)=find(z(ii)<=squeeze(zGr(iy(ii),ix(ii),:)),1,'last');
                end
            end
        end    
        iz(I)=min(Nz,iz(I));
    otherwise
        error('wrong number of input arguments, must be 2, 4 or 6');
end
