function n = hpindex11(E,x1,y1,z1,l,gs,depth);
% hpindex11(E,x1,y1,z1,l,gs,depth) - compute the HP index for
% function E over the cube with corner (x1,y1,z1) of length l.
% Each edge of the cube is subdivided into gs edges, with gs^2
% samples per face.  depth is the depth in the recursion, and
% currently only used for debugging statements.
% $Id: hpindex11.m,v 1.3 2001/09/10 11:19:39 smann Exp smann $

% Cell array p{} contains full vectors over faces.  We normalize
% these vectors, but retain normalization factors in array s()
% for later use.  We sum over the faces and normalize to get
% the Hopft-Poincare index.  This will detect point singularities.
%
% To detect line singularites, the vector field is restricted
% (projected) to each of the faces of the cube.  The projected
% vector field is stored in ee{}.  For each face, we sum over
% the boundary of the face.  This gives the Hopft-Poincare
% index of the vector field flowing through the face (eg, curve
% singularities) with a catch: potentially the projection step
% introduced a singularity in the vector field.  To detect this
% case, we first see if there is a zero vector on the face (in
% which case, we have our singularity); otherwise, we compare 
% the max/min lengths of the vectors over the face, and if the
% ratio is large enough, we assume there is a face singularity.
%
% We detect surface singularities in a manner similar to face
% singularities: we project the vector field onto the edges of
% the cube.  Now life is simpler: if the vectors change direction,
% we have a singularity, although again we have to avoid the
% new singularities introduced by the projection.

debug = 0;
x2 = x1+l; y2 = y1+l; z2 = z1+l;
xc = (x1+x2)/2;     yc = (y1+y2)/2;     zc = (z1+z2)/2;
dx = (x2-x1)/gs; dy = (y2-y1)/gs; dz = (z2-z1)/gs;
sum = 0;
for i=0:gs
  for j=0:gs
    [xx,yy,zz] = feval(E,x1+i*dx,y1+j*dy,z1);
    [pp,s(i+1,j+1)] = znormalize(xx+yy+zz);
    p{i+1,j+1} = pp;

    if (i==0|i==gs)|(j==0|j==gs)
      ee{i+1,j+1} = znormalize(xx+yy);
    end

    if (i==0)
      fy(1,j+1) = sign(inner(yy,e2));
      fys(1,j+1) = s(i+1,j+1);
    elseif (i==gs)
      fy(2,j+1) = sign(inner(yy,e2));
      fys(2,j+1) = s(i+1,j+1);
    end
    if (j==0)
      fx(1,i+1) = sign(inner(xx,e1));
      fxs(1,i+1) = s(i+1,j+1);
    elseif (j==gs)
      fx(2,i+1) = sign(inner(xx,e1));
      fxs(2,i+1) = s(i+1,j+1);
    end
  end
end
for i=1:gs
  for j=1:gs
    sum = sum+p{i,j}^p{i+1,j}^p{i+1,j+1};
    sum = sum+p{i,j}^p{i+1,j+1}^p{i,j+1};
    if debug
      p{i,j}^p{i+1,j}^p{i+1,j+1}+p{i,j}^p{i+1,j+1}^p{i,j+1}
    end
  end
end
if debug
  sum
end


fsum1 = 0;
B = e1^e2;
for i=1:gs
  fsum1 = fsum1+ee{i,1}^ee{i+1,1};
  fsum1 = fsum1+ee{1,i+1}^ee{1,i};
  fsum1 = fsum1+ee{i+1,gs+1}^ee{i,gs+1};
  fsum1 = fsum1+ee{gs+1,i}^ee{gs+1,i+1};
end
fsum1 = fsum1/2/pi/B;
if abs(fsum1) > .2
  mins = min(min(s));
  if mins > 1e-9
    maxs = max(max(s));
      fprintf(1,'1: min/max = %g/%g = %g\n',maxs,mins,maxs/mins);
    if maxs/mins < 10
      fsum1 = 0;
    end
  end
end

for i=0:gs
  for j=0:gs
    [xx,yy,zz] = feval(E,x1+i*dx,y1+j*dy,z2);
    [pp,s(i+1,j+1)] = znormalize(xx+yy+zz);
    p{i+1,j+1} = pp;

    if (i==0|i==gs)|(j==0|j==gs)
      ee{i+1,j+1} = znormalize(xx+yy);
    end

    if (i==0)
      fy(3,j+1) = sign(inner(yy,e2));
      fys(3,j+1) = s(i+1,j+1);
    elseif (i==gs)
      fy(4,j+1) = sign(inner(yy,e2));
      fys(4,j+1) = s(i+1,j+1);
    end
    if (j==0)
      fx(3,i+1) = sign(inner(xx,e1));
      fxs(3,i+1) = s(i+1,j+1);
    elseif (j==gs)
      fx(4,i+1) = sign(inner(xx,e1));
      fxs(4,i+1) = s(i+1,j+1);
    end

  end
end
for i=1:gs
  for j=1:gs
    sum = sum-p{i,j}^p{i+1,j}^p{i+1,j+1};
    sum = sum-p{i,j}^p{i+1,j+1}^p{i,j+1};
    if debug
      -(p{i,j}^p{i+1,j}^p{i+1,j+1}+p{i,j}^p{i+1,j+1}^p{i,j+1})
    end
  end
end
if debug
  sum
end

fsum2 = 0;
for i=1:gs
  fsum2 = fsum2+ee{i,1}^ee{i+1,1};
  fsum2 = fsum2+ee{1,i+1}^ee{1,i};
  fsum2 = fsum2+ee{i+1,gs+1}^ee{i,gs+1};
  fsum2 = fsum2+ee{gs+1,i}^ee{gs+1,i+1};
end
fsum2 = fsum2/2/pi/B;
if abs(fsum2) > .2
  mins = min(min(s));
  if mins > 1e-9
    maxs = max(max(s));
      fprintf(1,'2: min/max = %g\n',maxs/mins);
    if maxs/mins < 10
      fsum2 = 0;
    end
  end
end


for i=0:gs
  for j=0:gs
    [xx,yy,zz] = feval(E,x1,y1+i*dy,z1+j*dz);
    [pp,s(i+1,j+1)] = znormalize(xx+yy+zz);
    p{i+1,j+1} = pp;
    if (i==0|i==gs)|(j==0|j==gs)
      ee{i+1,j+1} = znormalize(yy+zz);
    end

    if (j==0)
      fz(1,i+1) = sign(inner(zz,e3));
      fzs(1,i+1) = s(i+1,j+1);
    elseif (j==gs)
      fz(2,i+1) = sign(inner(zz,e3));
      fzs(2,i+1) = s(i+1,j+1);
    end
  end
end
for i=1:gs
  for j=1:gs
    sum = sum+p{i,j}^p{i+1,j}^p{i+1,j+1};
    sum = sum+p{i,j}^p{i+1,j+1}^p{i,j+1};
    if debug
      p{i,j}^p{i+1,j}^p{i+1,j+1}+p{i,j}^p{i+1,j+1}^p{i,j+1}
    end
  end
end
if debug
  sum
end
fsum3 = 0;
B = e2^e3;
for i=1:gs
  fsum3 = fsum3+ee{i,1}^ee{i+1,1};
  fsum3 = fsum3+ee{1,i+1}^ee{1,i};
  fsum3 = fsum3+ee{i+1,gs+1}^ee{i,gs+1};
  fsum3 = fsum3+ee{gs+1,i}^ee{gs+1,i+1};
end
fsum3 = fsum3/2/pi/B;
if abs(fsum3) > .2
  mins = min(min(s));
  if mins > 1e-9
    maxs = max(max(s));
      fprintf(1,'3: min/max = %g\n',maxs/mins);
    if maxs/mins < 10
      fsum3 = 0;
    end
  end
end


for i=0:gs
  for j=0:gs
    [xx,yy,zz] = feval(E,x2,y1+i*dy,z1+j*dz);
    [pp,s(i+1,j+1)] = znormalize(xx+yy+zz);
    p{i+1,j+1} = pp;
    if (i==0|i==gs)|(j==0|j==gs)
      ee{i+1,j+1} = znormalize(yy+zz);
    end

    if (j==0)
      fz(3,i+1) = sign(inner(zz,e3));
      fzs(3,i+1) = s(i+1,j+1);
    elseif (j==gs)
      fz(4,i+1) = sign(inner(zz,e3));
      fzs(4,i+1) = s(i+1,j+1);
    end

  end
end
for i=1:gs
  for j=1:gs
    sum = sum-p{i,j}^p{i+1,j}^p{i+1,j+1};
    sum = sum-p{i,j}^p{i+1,j+1}^p{i,j+1};
    if debug
      -(p{i,j}^p{i+1,j}^p{i+1,j+1}+p{i,j}^p{i+1,j+1}^p{i,j+1})
    end
  end
end
if debug
  sum
end
fsum4 = 0;
for i=1:gs
  fsum4 = fsum4+ee{i,1}^ee{i+1,1};
  fsum4 = fsum4+ee{1,i+1}^ee{1,i};
  fsum4 = fsum4+ee{i+1,gs+1}^ee{i,gs+1};
  fsum4 = fsum4+ee{gs+1,i}^ee{gs+1,i+1};
end
fsum4 = fsum4/2/pi/B;
if abs(fsum4) > .2
  mins = min(min(s));
  if mins > 1e-9
    maxs = max(max(s));
      fprintf(1,'4: min/max = %g\n',maxs/mins);
    if maxs/mins < 10
      fsum4 = 0;
    end
  end
end



for i=0:gs
  for j=0:gs
    [xx,yy,zz] = feval(E,x1+i*dx,y1,z1+j*dz);
    [pp,s(i+1,j+1)] = znormalize(xx+yy+zz);
    p{i+1,j+1} = pp;
    if (i==0|i==gs)|(j==0|j==gs)
      ee{i+1,j+1} = znormalize(xx+zz);
    end
  end
end
for i=1:gs
  for j=1:gs
    sum = sum-p{i,j}^p{i+1,j}^p{i+1,j+1};
    sum = sum-p{i,j}^p{i+1,j+1}^p{i,j+1};
    if debug
      -(p{i,j}^p{i+1,j}^p{i+1,j+1}+p{i,j}^p{i+1,j+1}^p{i,j+1})
    end
  end
end
if debug
  sum
end
fsum5 = 0;
B = e3^e1;
for i=1:gs
  fsum5 = fsum5+ee{i,1}^ee{i+1,1};
  fsum5 = fsum5+ee{1,i+1}^ee{1,i};
  fsum5 = fsum5+ee{i+1,gs+1}^ee{i,gs+1};
  fsum5 = fsum5+ee{gs+1,i}^ee{gs+1,i+1};
end
fsum5 = fsum5/2/pi/B;
if abs(fsum5) > .2
  mins = min(min(s));
  if mins > 1e-9
    maxs = max(max(s));
    fprintf(1,'HPindex11: %g/%g = %g (%d)\n',maxs,mins,maxs/mins,depth);
    if maxs/mins < 10
      fsum5 = 0;
    end
  end
end


for i=0:gs
  for j=0:gs
    [xx,yy,zz] = feval(E,x1+i*dx,y2,z1+j*dz);
    [pp,s(i+1,j+1)] = znormalize(xx+yy+zz);
    p{i+1,j+1} = pp;
    if (i==0|i==gs)|(j==0|j==gs)
      ee{i+1,j+1} = znormalize(xx+zz);
    end
  end
end
for i=1:gs
  for j=1:gs
    sum = sum+p{i,j}^p{i+1,j}^p{i+1,j+1};
    sum = sum+p{i,j}^p{i+1,j+1}^p{i,j+1};
    if debug
      p{i,j}^p{i+1,j}^p{i+1,j+1}+p{i,j}^p{i+1,j+1}^p{i,j+1}
    end
  end
end
if debug
  sum
end

fsum6 = 0;
for i=1:gs
  fsum6 = fsum6+ee{i,1}^ee{i+1,1};
  fsum6 = fsum6+ee{1,i+1}^ee{1,i};
  fsum6 = fsum6+ee{i+1,gs+1}^ee{i,gs+1};
  fsum6 = fsum6+ee{gs+1,i}^ee{gs+1,i+1};
end
fsum6 = fsum6/2/pi/B;
if abs(fsum6) > .2
  mins = min(min(s));
  if mins > 1e-9
    maxs = max(max(s));
      fprintf(1,'6: min/max = %g\n',maxs/mins);
    if maxs/mins < 10
      fsum6 = 0;
    end
  end
end

% Compute surface singularities by checking along lines of cube.
seps=1e-09;
sratio=40;

%fx
%fxs
fxhp = [];
for i=1:4
  if min(fxs(i,:)) < seps
    fprintf(1,'ax %d %g\n',i,min(fxs(i,:)));
    fxhp = [fxhp,1];
  else
    if (max(1-fx(i,:))==2 & min(1-fx(i,:))==0)&...
	  max(fxs(i,:))/min(fxs(i,:)) > sratio
      fprintf(1,'bx %d %g\n',i,max(fxs(i,:))/min(fxs(i,:)));
      fxhp = [fxhp,1];
    else
      fprintf(1,'cx %d %g\n',i,max(fxs(i,:))/min(fxs(i,:)));
      fxhp = [fxhp,0];
    end
  end
end
%fxhp

fyhp = [];
for i=1:4
  if min(fys(i,:)) < seps
    fprintf(1,'ay %d %g\n',i,min(fys(i,:)));
    fyhp = [fyhp,1];
  else
    if (max(1-fy(i,:))==2 & min(1-fy(i,:))==0)&...
	  max(fys(i,:))/min(fys(i,:)) > sratio
      fprintf(1,'by %d %g\n',i,max(fxs(i,:))/min(fxs(i,:)));
      fyhp = [fyhp,1];
    else
      fprintf(1,'cy %d %g\n',i,max(fys(i,:))/min(fys(i,:)));
      fyhp = [fyhp,0];
    end
  end
end

fzhp = [];
for i=1:4
  if min(fzs(i,:)) < seps
    fprintf(1,'az %d %g\n',i,min(fzs(i,:)));
    fzhp = [fzhp,1];
  else
    if (max(1-fz(i,:))==2 & min(1-fz(i,:))==0)&...
	  max(fzs(i,:))/min(fzs(i,:)) > sratio
      fprintf(1,'bz %d %g\n',i,max(fxs(i,:))/min(fxs(i,:)));
      fzhp = [fzhp,1];
    else
      fprintf(1,'cz %d %g\n',i,max(fzs(i,:))/min(fzs(i,:)));
      fzhp = [fzhp,0];
    end
  end
end

clear p,ee;

% bundle singularity values into array for return
n = [I3*sum/6/(4/3)/pi, fsum1, fsum2, fsum3, fsum4, fsum5, fsum6, ...
     fxhp, fyhp, fzhp]
