function [x,y,z,xm,ym,zm,dx,dy,dz,Nx,Ny,Nz,M]=modelsize3D(x,y,z)
%[x,y,xm,ym,dx,dy,Nx.Ny,M]=modelsize(x,y)
% makes values in x and y unique and sorted with correct orientation
% then compute the other values
% xm,ym cell center cordinates
% dx, dy cell widths (positivie)
% Nx, Ny size grid 
% M=ones(Ny,Nx)
% TO 2004


% Copyright 2009 Theo Olsthoorn, TU-Delft and Waternet, without any warranty
% under free software foundation GNU license version 3 or later
x=unique(x(:)');
y=unique(y(:)');

if ndims(z)<3
    z=sort(unique(z(:)),'descend');
end

dx=abs(diff(x));
dy=abs(diff(y));
Nx=length(dx);
Ny=length(dy);
Nz=length(z)-1;

xm=0.5*(x(1:end-1)+x(2:end));
ym=0.5*(y(1:end-1)+y(2:end));

dz=zeros(Nz,1);
zm=zeros(Nz,1);

k=0;
for i=1:Nz
    k=k+1;
    dz(i)=abs(z(k)-z(k+1));
    zm(i)=0.5*(z(k)+z(k+1));
    if ICBD(i), k=k+1; end
end

M=ones(Ny,Nx,Nz);
