function [fvals]=contour3d(x,y,z,f,nf,varargin)
// Draws 3D contour-/iso-surfaces of function f on 3D lattice x,y,z.
// Description: contour3d  draws level surfaces of a 3D function  f(x,y,z)
// on a 3D plot. The values of  f(x,y,z)  are in the 3D hypermatrix  f  
// at the lattice grid points defined by 1D or 3D matrices  x , y and  z. 
// Set any missing f-values to %nan.
// Extra arguments are passed on to plot3d to set view, legend, etc.
// 
// Uses 40.nx.ny.nz memory or thereabouts.
// Tony Roberts, July 2006.

// first get & check sizes of 3D arrays
if length(size(f))~=3, error("contour3d: f must be 3D hypermatrix"), end
if min(size(f))<2, error("contour3d: need >=2 points in each dimension"), end
[nx,ny,nz]=size(f);
if length(x(:))==nx, x=hypermat([nx,1,1],x(:)); x=x(:,ones(1,ny),ones(1,nz)); end
if length(y(:))==ny, y=hypermat([1,ny,1],y(:)); y=y(ones(1,nx),:,ones(1,nz)); end
if length(z(:))==nz, z=hypermat([1,1,nz],z(:)); z=z(ones(1,nx),ones(1,ny),:); end
if size(x)~=[nx,ny,nz],  error("contour3d: x 3D matrix does not match f"), end
if size(y)~=[nx,ny,nz],  error("contour3d: y 3D matrix does not match f"), end
if size(z)~=[nx,ny,nz],  error("contour3d: z 3D matrix does not match f"), end
// find iso-surface values to include zero
if (length(nf)==1)&(nf==round(nf))&(nf>0)
  minf=min(f);
  maxf=max(f);
  df=(maxf-minf)/nf;
  fvals=df*(ceil(minf/df):floor(maxf/df));
  maxf=fvals($)+df/2; minf=fvals(1)-df/2;
else // use contour values in nf
  fvals=nf;
  maxf=max(nf); minf=min(nf);
  df=(maxf-minf)/(length(nf)-1+%eps);
  maxf=maxf+df/2; minf=minf-df/2;
end
// make a huge list of tetrahedrons
// use these to cater for 8 different orientations of the 5 tetras
i=2*ceil((1:nx-1)/2); id=(-1).^(1:nx-1);
j=2*ceil((1:ny-1)/2); jd=(-1).^(1:ny-1);
k=2*ceil((1:nz-1)/2); kd=(-1).^(1:nz-1);
l=(nx-1)*(ny-1)*(nz-1); // number of 3D bricks
ntet=5*l;  // number of tetras
// define tetras by their lattice point labels
p=hypermat([nx,ny,nz],1:nx*ny*nz); // index of points
ps=zeros(ntet,4);
ps(:,1)=matrix([p(i,j,k);p(i+id,j,k+kd);p(i,j+jd,k+kd);p(i+id,j+jd,k);p(i+id,j,k)],ntet,1);
ps(:,2)=matrix([p(i+id,j,k);p(i,j,k+kd);p(i,j,k+kd);p(i,j+jd,k);p(i,j+jd,k)],ntet,1);
ps(:,3)=matrix([p(i,j+jd,k);p(i+id,j,k);p(i,j+jd,k);p(i+id,j,k);p(i,j,k+kd)],ntet,1);
ps(:,4)=matrix([p(i,j,k+kd);p(i+id,j+jd,k+kd);p(i+id,j+jd,k+kd);p(i+id,j+jd,k+kd);p(i+id,j+jd,k+kd)],ntet,1);
fs=[f(ps(:,1)) f(ps(:,2)) f(ps(:,3)) f(ps(:,4))];
// sort each tetra corner into decreasing order of f
[fs,ll]=sort(fs,"c");
ps(:)=ps(sub2ind([ntet,4],(1:ntet)'*ones(1,4),ll));
// only use tetras that have all vertices defined
ok=find(~isnan(sum(fs,"c")));

// build facets for each isosurface value
xf=[]; yf=[]; zf=[]; cf=[];
for col=1:length(fvals)
  fval=fvals(col);
  i=sum(fs(ok,:)>fval,"c"); // how many values>contour
  // do all the 1-3 tetras
  j=ok(find(i==1)); 
  if ~isempty(j) then k=ones(1,4); l=[2 2:4];
  rats=(fval-fs(j,l))./(fs(j,k)-fs(j,l));
  xf=[xf; rats.*matrix(x(ps(j,k)),-1,4)+(1-rats).*matrix(x(ps(j,l)),-1,4)];
  yf=[yf; rats.*matrix(y(ps(j,k)),-1,4)+(1-rats).*matrix(y(ps(j,l)),-1,4)];
  zf=[zf; rats.*matrix(z(ps(j,k)),-1,4)+(1-rats).*matrix(z(ps(j,l)),-1,4)];
  cf=[cf col(ones(j))];
  end
  // do all the 2-2 tetras
  j=ok(find(i==2)); 
  if ~isempty(j) then k=[1 1 2 2]; l=[3 4 4 3];
  rats=(fval-fs(j,l))./(fs(j,k)-fs(j,l));
  xf=[xf; rats.*matrix(x(ps(j,k)),-1,4)+(1-rats).*matrix(x(ps(j,l)),-1,4)];
  yf=[yf; rats.*matrix(y(ps(j,k)),-1,4)+(1-rats).*matrix(y(ps(j,l)),-1,4)];
  zf=[zf; rats.*matrix(z(ps(j,k)),-1,4)+(1-rats).*matrix(z(ps(j,l)),-1,4)];
  cf=[cf col(ones(j))];
  end
  // do all the 3-1 tetras
  j=ok(find(i==3)); 
  if ~isempty(j) then k=[1 1:3]; l=4*ones(1,4);
  rats=(fval-fs(j,l))./(fs(j,k)-fs(j,l));
  xf=[xf; rats.*matrix(x(ps(j,k)),-1,4)+(1-rats).*matrix(x(ps(j,l)),-1,4)];
  yf=[yf; rats.*matrix(y(ps(j,k)),-1,4)+(1-rats).*matrix(y(ps(j,l)),-1,4)];
  zf=[zf; rats.*matrix(z(ps(j,k)),-1,4)+(1-rats).*matrix(z(ps(j,l)),-1,4)];
  cf=[cf col(ones(j))];
  end
end

// draw factes according to current colormap
colorbar(minf,maxf)
plot3d(xf',yf',list(zf',cf),varargin(:))
//  ,flag=[1,3,4],ebox=[min(x) max(x) min(y) max(y) min(z) max(z)]);
fog=gcf();
fog=fog.children(1).children;
fog.hiddencolor=-1; // make front and back same color

endfunction