function I=cellIndex(varargin)
% indices=cellIndex(ix,iy,Nx,Ny)
% indices=cellIndex(ix,iy,iz,Nx,Ny,Nz)
% compute the global indices of the cell indices in the array with given dimension Nx Ny Nz
% ix,iy,iz are vectors of equal length with index pairs
% TO 090317 090509 100521


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

if length(varargin)==4
    ix=varargin{1}; iy=varargin{2}; Nx=varargin{3}; Ny=varargin{4};
    if length(ix)==1, ix=ones(size(iy))*ix; end
    if length(iy)==1, iy=ones(size(ix))*iy; end
    I=NaN(size(ix(:)));
    if max(iy)>Ny || min(iy)<1, error('iy out of bound'); end
    if max(ix)>Nx || min(ix)<1, error('ix out of bound'); end
    for i=1:length(ix)
        I(i)=Ny*(ix(i)-1)+iy(i);
    end
elseif length(varargin)==6
    ix=varargin{1}; iy=varargin{2}; iz=varargin{3};
    if length(ix)==1 && length(iy)==1, ix=ones(size(iz))*ix; iy=ones(size(iz))*iy; end
    if length(ix)==1 && length(iz)==1, ix=ones(size(iy))*ix; iz=ones(size(iy))*iz; end
    if length(iy)==1 && length(iz)==1, iy=ones(size(ix))*iy; iz=ones(size(ix))*iz; end
    I=NaN(size(ix(:)));
    Nx=varargin{4}; Ny=varargin{5}; Nz=varargin{6};
    if max(ix)>Nx || min(ix)<1, error('ix out of bound'); end
    if max(iy)>Ny || min(iy)<1, error('iy out of bound'); end
    if max(iz)>Nz || min(iz)<1, error('iz out of bound'); end
    for i=1:length(ix)
        I(i)=Ny*Nx*(iz(i)-1)+Ny*(ix(i)-1)+iy(i);
    end
else
    help('cellIndex');
    error('See help for this function');
end
