function [Tx,Ty,Tz] = surfcv4(x,y,z,f,n,const,col);
%SURFCV4 3-D surface of constant value plot (Matlab 4 version).
%   SURFCV4(X,Y,Z,F,N,CONST,col) draws a surface of constant value
%   for a function of three variables:  F(X,Y,Z) = CONST .
%   The arrays X,Y,Z define the coordinates for F and
%   must be monotonic volume arrays (as if produced by
%   meshgrid).  The surface is calculated by linear
%   interpolation of the values of F and is drawn with 
%   triangular patches using FILL3.  The color of each 
%   point on the surface is determined by its distanse 
%   from the origin.  N is the number of columns in the 
%   M-by-N-by-P volume array F.
%
%   SURFCV4(F,N,CONST) assumes X=1:N, Y=1:M, Z=1:P.
%
%   H = SURFCV4(...) returns a column vector of handles 
%   to PATCH objects, one handle per patch.
%
%   [TX,TY,TZ] = SURFCV4(X,Y,Z,F,CONST); returns 3-by-K
%   arrays specifying x,y,z coordinates of the vertices
%   of the triangular patches, where K is the number
%   of patches.  FILL3(TX,TY,TZ,C) can be used to draw
%   the surface.
%
%   Example: To visualize the surface 
%        1./(x.^2+y.^2) +  1./(y.^2+z.^2) + 1./(z.^2+x.^2) = 2
%   over the range -3 < x < 3, -3 < y < 3, -3 < z < 3, 
%
%      [x,y,z] = meshgrid(-3:.4:3, -3:.6:3, -3:.4:3);
%      f = 1./(x.^2+y.^2) +  1./(y.^2+z.^2) + 1./(z.^2+x.^2);
%      n = 16;  const = 2;
%      surfcv4(x,y,z,f,n,const);

%  Copyright (c) 1997 Ruslan L. Davidchack
%                     University of Kansas, Lawrence, KS
%                     davidchack@kuphsx.phsx.ukans.edu
%      created:  Sep 29, 1997; 
%      modified for Matlab 4:  Oct 31, 1997 

%Short description of the algorithm.
%   Function F is defined on the M-by-N-by-P 3-D grid which 
%   divides the space into (M-1)*(N-1)*(P-1) cubes.   
%   F specifies values of the function at the vertices.
%   The surface of constant value cuts a patch through a 
%   cube if F < CONST for some vertices of the cube and 
%   F > CONST for others. 
%     Note: F is slightly modified on the input to assure 
%           that it does not exactly equal CONST on a vertex,
%           and thus a vertex of a surface element belongs to
%           only one edge (see below). 
%   Every vertex of the patch lays on the edge of the 
%   cube and its coordinates can be determined by linear 
%   interpolation between the values of the function on the
%   ends of the edge. Two vertices of a patch are connected 
%   if they lay on the edges that belong to the same face of
%   a cube. If a patch has more than three vertices, it is
%   subdivided into triangular elements by additional 
%   connections between vertices that do not belong to the
%   same face.
%   The idea of the algorithm is to group the edges which
%   contain vertices of the surface patches into pairs according 
%   to their affiliation with the faces (there can be only 2 or 4   
%   patch vertices on a face).  Then the edge pairs are grouped
%   according to their affiliation with the cubes.  The "cube-edge" 
%   array, thus created, is then decomposed into triplets of patch
%   vertices.  There can be more than one patch per cube and 
%   the algorithm handles this case as well.
%Naming convensions
%   ve - vertex, ed - edge, fa - face, cu - cube
%   nv - number of vertices, ne - number of edges, etc.
%   ce - coordinates of edges, cf - coordinates of faces, etc.
%   ie - indices of edges, ia - indices of faces, etc.

%Simple case to follow the algorithm
%      [x,y,z] = meshgrid(1:4, 1:3, 1:2); n = 4; const = 15;
%      f = x.^2 + y.^2 + z.^2; 
%Different cool surfaces:
% 1.   [x,y,z] = meshgrid(-3:.2:3, -3:.2:3, -3:.2:1.6); 
%      n = 31; const = 7;
%      f = 1./(x.^2+y.^2)+1./(y.^2+z.^2)+1./(z.^2+x.^2)+x.^2+y.^2+z.^2;
%     
% 2.   [x,y,z] = meshgrid(-2:.2:2, -2:.2:2, -2:.2:2);  
%      n = 21; const = 0;
%      f = (x.^2 + y.^2 + z.^2).^2 - 4*(x.^2 - y.^2 - z.^2);
%
% 3.   [x,y,z] = meshgrid(-2:.2:2, -2:.2:2, -2:.2:2);  
%      n = 21; const = 0;
%      f = x.^2 + y.^2 - z.^2.*(2 - z)./(2 + z);                 

  if nargin == 3,
    f = x; n = y; const = z;
    [x,y,z] = meshgrid(1:n,1:round(size(f,1)/n),1:size(f,2));
  elseif nargin < 6 | nargin > 7
    disp('Wrong number of input arguments'); return;
  end,  

%Make sure F == CONST is never true
  f1 = f;  i1 = find(f == const);
  if const == 0
    f1(i1) = 1e-10.*ones(size(i1));
  else 
    f1(i1) = (1 + 1e-10).*const(ones(size(i1)));
  end,

%Number of vertices in each direction:  nv = [M N P]
%   Total number of vertices = prod(nv)
  nv = [0 n 0]; 
  [i1 nv(3)] = size(f1); nv(1) = round(i1/n);
  iv = reshape(1:prod(nv),i1,nv(3));
  [ix,iy,iz] = meshgrid(1:nv(2),1:nv(1),1:nv(3));

%Number of edges of each orientation:  
%			ne = [(M-1)*N*P  M*(N-1)*P  M*N*(P-1)]
%   Total number of edges = sum(ne)
  ne = prod(nv) - prod(reshape(nv([2 1 1 3 3 2])',3,2)');

%Determine indices of vertices that belong to each edge
%   Indexing of vertices is such as in f(:)
  ed = zeros(sum(ne),2);
%   Edges parallel to y-axis
  i1 = reshape(iv,nv(1),nv(2)*nv(3));
  ed((1:ne(1)),1) = reshape(i1(1:nv(1)-1,:),ne(1),1); 
  ed((1:ne(1)),2) = ed((1:ne(1)),1) + 1;
%   Edges parallel to x-axis
  ie = (1:ne(2)) + ne(1);
  ed(ie,1) = reshape(iv(1:(nv(1)*(nv(2)-1)),:),ne(2),1);
  ed(ie,2) = ed(ie,1) + nv(1); 
%   Edges parallel to z-axis
  ie = (1:ne(3)) + ne(1) + ne(2);
  ed(ie,1) = (1:ne(3))';
  ed(ie,2) = ed(ie,1) + nv(1)*nv(2);

%Select edges crossed by the surface 
  se = find(prod(f1(ed)' - const) < 0)';
  if isempty(se), disp('There is no surface in the specified range');
    return; end;
  ed = ed(se,:);

%Restore coordinates of the selected edges 
  ce = zeros(length(se),4);
  ce(:,4) = (se > ne(1)) + (se > ne(1)+ne(2)) + 1;
  ce(:,1:3) = [ix(ed(:,1)) iy(ed(:,1)) iz(ed(:,1))];
 
%Determine coordinates of the patch vertices  
  ie = find(diff([ce(:,4);4]) > 0); 
  ie1 = 1:ie(1);  ie2 = ie(1)+1:ie(2);  ie3 = ie(2)+1:ie(3);
  ds = (const - f1(ed(:,1)))./diff(f1(ed)')';
  le = [x(nv(2)+1,1)-x(1,1);y(2,1)-y(1,1);z(1,2)-z(1,1)];
  xv=x(nv(1)*ce(:,1),1); xv=xv(:); xv(ie2) = xv(ie2) + ds(ie2).*le(1); 
  yv = y(ce(:,2),1);   yv = yv(:); yv(ie1) = yv(ie1) + ds(ie1).*le(2); 
  zv = z(1,ce(:,3));   zv = zv(:); zv(ie3) = zv(ie3) + ds(ie3).*le(3); 

%Find faces containing the selected edges (maximum of 4 faces per edge)
  cf = [(1:length(se))' ce]; cf = cf(:,[1:5 1:5 1:5 1:5]);
%Edge-to-face coordinate transformation matrix
  etof = [0 0 0 0 0   0 -1  0  0 0   0 0 0 0  2   0  0  0 -1  2;...
	  0 0 0 0 0   0  0  0 -1 0   0 0 0 0 -1   0  0 -1  0 -1;...
          0 0 0 0 0   0  0 -1  0 0   0 0 0 0 -1   0 -1  0  0 -1];
  ie = diff([0;ie]);
  cf = cf + [(etof(1,:)'*ones(1,ie(1)))';...
             (etof(2,:)'*ones(1,ie(2)))';...
             (etof(3,:)'*ones(1,ie(3)))'];
  cf = reshape(cf',5,4*length(se));

%Discard faces that are outside the range of x,y,z 
  cf(:,find(prod(cf) == 0)) = [];  cf = cf';
  i1 = find(cf(:,5) ~= 2);  cf(i1(find(cf(i1,3) > nv(1)-1)),:) = [];  
  i2 = find(cf(:,5) < 3);   cf(i2(find(cf(i2,2) > nv(2)-1)),:) = [];  
  i3 = find(cf(:,5) > 1);   cf(i3(find(cf(i3,4) > nv(3)-1)),:) = [];

%Find indices of the selected faces
  i1 = find(cf(:,5)==1); i2 = find(cf(:,5)==2); i3 = find(cf(:,5)==3);
%Number of faces of each orientation:
%                ne = [(M-1)*(N-1)*P  M*(N-1)*(P-1)  (M-1)*N*(P-1)]
  nf = [nv(3)*(nv(2)-1)*(nv(1)-1) nv(1)*(nv(2)-1)*(nv(3)-1)...
        nv(2)*(nv(3)-1)*(nv(1)-1)];
  fa = zeros(size(cf,1),1);  
  fa(i1) = cf(i1,3)+((cf(i1,2)-1)+(cf(i1,4)-1)*(nv(2)-1))*(nv(1)-1);
  fa(i2) = cf(i2,3)+((cf(i2,2)-1)+(cf(i2,4)-1)*(nv(2)-1))*nv(1)+nf(1);
  fa(i3) = cf(i3,3)+((cf(i3,2)-1)+(cf(i3,4)-1)*nv(2))*(nv(1)-1)+...
	   nf(1) + nf(2);
  [fa ia] = sort(fa);

%!!!!!!!! Check whether fac values come in pairs.
%   Each face can have only 2 or 4 selected edges.
%   The check below can be removed after debugging process 
%   is completed
if rem(length(fa),2) ~= 0,
  disp('Odd number!'); return,
else 
  if sum(diff(reshape(fa,2,length(fa)/2))) > 0,
  disp('Odd number!'); return, end, end;
%!!!!!!!!

%Combine face coord. with indices of edges that belong to the faces
  cf = [fa(1:2:length(fa)) reshape(cf(ia,1),2,length(fa)/2)' ...
        cf(ia(1:2:length(ia)),2:5)];

%Find cubes containing selected faces
  i1 = find(cf(:,7)==1); i2 = find(cf(:,7)==2); i3 = find(cf(:,7)==3);
  cc = cf(:,[2:6 2:6]);
  ftoc = [zeros(3,7) [0 0 -1;0 -1 0;-1 0 0]];
  cc = cc + [(ftoc(1,:)'*ones(1,length(i1)))';...
             (ftoc(2,:)'*ones(1,length(i2)))';...
             (ftoc(3,:)'*ones(1,length(i3)))'];
  cc = reshape(cc',5,size(cf,1)*2);
  cc(:,find(prod(cc) == 0 | cc(4,:) > nv(1)-1 | cc(3,:) > nv(2)-1 |...
            cc(5,:) > nv(3)-1)) = [];  cc = cc';

  cu = cc(:,4) + ((cc(:,3)-1) + (cc(:,5)-1)*(nv(2)-1))*(nv(1)-1);
  [cu ic] = sort(cu); 
%Combine indices of selected cubes with indices of edges that belong
%   to them into a cube-edge index array
  cc = [cu cc(ic,1:2)];

%Combine edges into triplets composing surface patches
  ip = [];   cc = sortrows(cc);
  while 1,
%Indices of patches
    ic = find(diff([0;cc(:,1);0]) ~= 0);
%Indices of triangular patches
    ic3 = ic(find(diff(ic) == 3));
%Indices of patches with more than 3 vertices
    icm = ic(find(diff(ic) > 3));   
    
%Select triangular patches and reduce 
    ic(length(ic)) = [];   ip = [ip [cc(ic,2:3)';cc(ic+1,3)']];
    if ~isempty(ic3),
      irm = ic3*ones(1,3) + ones(length(ic3),1)*(0:2);
    else irm = []; end,
    irm = [irm(:);icm];   cc(icm+1,2) = cc(icm,3);   cc(irm,:) = [];

%Sort new cube-edge index array or quit if it is empty
    if ~isempty(cc),  cc = sortrows(cc);  else, 
      disp(['Number of patches = ' num2str(size(ip,2))]); 
      break; end

%Find cases when there is more than one patch per cube 
    irm = find(sum(abs(diff(cc)')) == 0);
    if ~isempty(irm),
      %disp(['More than one patch in ' int2str(length(irm)) ' cubes']);
      irm = [1;1]*irm + [0;1]*ones(1,length(irm)); irm = irm(:);
      cc(irm,:) = []; end,

  end;

  Tx = reshape(xv(ip),size(ip));
  Ty = reshape(yv(ip),size(ip));  
  Tz = reshape(zv(ip),size(ip));

%This is to compute area of each patch
%  x2 = Tx(2,:)-Tx(1,:); y2 = Ty(2,:)-Ty(1,:); z2 = Tz(2,:)-Tz(1,:); 
%  x3 = Tx(3,:)-Tx(1,:); y3 = Ty(3,:)-Ty(1,:); z3 = Tz(3,:)-Tz(1,:);
%A=sqrt((x2.*y3-x3.*y2).^2+(y2.*z3-y3.*z2).^2+(z2.*x3-z3.*x2).^2)./2; 

  if nargin==3 | nargin==6,
    C = sqrt(Tx.^2 + Ty.^2 + Tz.^2);
    H = fill3(Tx,Ty,Tz,C);
    else,  
      for i=1:length(Tx),
        fill3(Tx(:,i),Ty(:,i),Tz(:,i),col);
        if i==1,
          hold on;
        end;
      end;
      hold off;
  end;

  set(gca,'xlim',[x(1,1) x(nv(1)*nv(2),1)],...
          'ylim',[y(1,1) y(nv(1),1)],...
          'zlim',[z(1,1) z(1,nv(3))],'box','on');