function ntz = bztest(rk,b)
%
% This function determines if vector rk lies within or at the border of the 
% first Brillouin zone. rk-components of reciprocal lattice vectors
% (b1,b2,b3) = (x,y,z)*b. cark are the cartesian components of rk.
% ntz=-1 if within the Brillouin zone
% ntz=0 if on the border of the Brillouin zone
% ntz=1 if outside of the Brillouin zone
%
sign = ones(1,3);
for I = 1:3
  if rk(I) < 0
    sign(I) = -1;
  end
end
cark(1:3) = b(1:3,1:3)*rk(1:3).';
nbz = 0;
for I1 = 1:7
  if I1 <= 3
    c(1:3) = b(1:3,I1)*sign(I1);
  end
  if I1 > 3
    if I1 < 7
      if I1 == 4
        K1 = 1;
        K2 = 2;
      end
      if I1 == 5
        K2 = 3;
      end
      if I1 == 6
        K1 = 2;
      end
%
% Vectors b1+b2, b1+b3, b2+b3, with adapted signs
%
      c(1:3) = b(1:3,K1)*sign(K1) + b(1:3,K2)*sign(K2);
    else
%
% Vector b1+b2+b3, with adapted signs.
%
      c(1:3) = b(1:3,1:3)*sign(1:3).';
    end
  end
%
% Test if the projection of rk on vector c is less than half the length
% of c.
%
  bnorm = 0.5*sum(c(1:3).^2);
%
% The adapted signs ensure tests in the right octant.
%
  ckinp = cark(1:3)*c(1:3).';
  ckinp = ckinp - bnorm;
  if ckinp < -0.001
    nbz = nbz + 1;
  end
  if ckinp > 0.001
%
% Outside the Brillouin zone.
%
    ntz = 1;
    return
  end
end
if nbz == 7
%
% When all seven tests are fulfilled, then the rk lies inside the BZ.
%
  ntz = -1;
  return
else
%
% Otherwise, rk lies on the border of the Brillouin zone.
%
  ntz = 0;
  return
end
