function recur3(E,x1,y1,z1,l,d,movieflag)
% recur3(E,x1,y1,z1,l,d) - recursively descend the octree with one
%  corner at (x1,y1,z1) of size l to max depth d, evaluating the
%  HP index of function E for each cube (E is a string).  If the
%  function is 0 on a cube, then stop descending.
% $Id: recur3.m,v 1.4 2001/09/07 21:02:39 smann Exp smann $

persistent bcount;
persistent movieframe;
persistent movie;
gs = 9;
if nargin < 4
  fprintf(2,'Format: recur3(''E'',x1,y1,z1,[l=1,d=-2,movieflag=0])\n');
  return;
end
if nargin < 5
  l = 1;
end
if nargin < 6
  d = -2;
end
if nargin < 7
  movieflag = 0;
end
if size(bcount)==0
  bcount =0;
  movieframe = 1;
end
if d < 0
fprintf(1,'recur3: $Id: recur3.m,v 1.4 2001/09/07 21:02:39 smann Exp smann $\n');
  movieframe = 1;
  clear movie;
  bcount = 0;
  d = -d;
  fid = fopen('trace','w+');
  fclose(fid);
  clf;
end

fprintf(1, 'hpindex11(%s,%g,%g,%g,%g,%d) bcount=%d\n',E,x1,y1,z1,l,gs,bcount);
h = gbox(x1,y1,z1,x1+l,y1+l,z1+l);
if d>2 | bcount<40 | mod(bcount,4)==0 | movieflag
  axis equal;drawnow;
  if movieflag
    imwrite(frame2im(getframe),sprintf('Movie/frame%04d.png',movieframe),'png');
    movieframe = movieframe+1;
  end
end
bcount = bcount+1;
if d > 0
  ind = abs(hpindex11(E,x1,y1,z1,l,gs,d));
else
  ind = abs(hpindex11(E,x1,y1,z1,l,2*gs,d));
end

fprintf(1,'recur3(E,%g,%g,%g, %g, %d => \n      [%g %g %g %g %g %g %g])\n',...
	x1,y1,z1,l,d,ind(1),ind(2),ind(3),ind(4),ind(5),ind(6),ind(7));
delete(h);
box(x1,y1,z1,x1+l,y1+l,z1+l);
if d == 0 
  if max(ind) > .2
fprintf(1,'draw!\n');
    if ind(1) > .2
      cube(x1,y1,z1,x1+l,y1+l,z1+l,'r');
    elseif max(ind(2:7)) > .2
      cube(x1,y1,z1,x1+l,y1+l,z1+l,'c');
    else
      cube(x1,y1,z1,x1+l,y1+l,z1+l,'y');
    end
    drawnow;
    if movieflag
      imwrite(frame2im(getframe),sprintf('Movie/frame%04d.png',movieframe),'png');
      movieframe = movieframe+1;
    end
    bcount = bcount+1;
    fid = fopen('trace','a');
    fprintf(fid,'%g %g %g %g: [%g %g %g %g %g %g %g])\n',...
	    x1,y1,z1,l,ind(1),ind(2),ind(3),ind(4), ...
	    ind(5),ind(6),ind(7));
    fclose(fid);
  end
else
  if max(ind) > .1
    recur3(E,x1,y1,z1,l/2,d-1,movieflag);
    recur3(E,x1+l/2,y1,z1,l/2,d-1,movieflag);
    recur3(E,x1+l/2,y1+l/2,z1,l/2,d-1,movieflag);
    recur3(E,x1,y1+l/2,z1,l/2,d-1,movieflag);

    recur3(E,x1,y1,z1+l/2,l/2,d-1,movieflag);
    recur3(E,x1+l/2,y1,z1+l/2,l/2,d-1,movieflag);
    recur3(E,x1+l/2,y1+l/2,z1+l/2,l/2,d-1,movieflag);
    recur3(E,x1,y1+l/2,z1+l/2,l/2,d-1,movieflag);
  end
end
