function out = NNscalarproduct(X,xmin,xmax,ymin,ymax,zmin,zmax)

% This routine calculates dot products of nearest neighbor vectors 
% to identify stacking faults in an fcc crystal. 
% Use the "showdefects.m" routine to visualize the stacking fault in 
% three dimensions.
% X is an N by 3 list of x,y,z coordinates of all N particles of the crystal.
% xmin,xmax,ymin,ymax,zmin,zmax define the region, in which the program
% calculates the nearest neighbor configuration (you should exclude the edges).
% The output is a 7 by N list containing the number of nearest neighbors,
% and the 6 smallest dotproducts of nearest neighbor vectors.
% (that's how stacking faults are identified: unlike the perfect fcc
% lattice where particles have 6 opposing nearest neighbors, at the
% stacking fault, particles have merely 3 nearest neighbors)
% Please set the nearest neighbor distance d0 and the nearest neighbor
% acceptance range below.
%
% Peter Schall, 2004
%
%
% Set the nearest neighbor distance d0 and the acceptance region 
% for nearest neighbors here
d0 = 1.63;
min = 0.8 * d0;
max = 1.26 * d0;

% initialization of variables
[N,M] = size(X);
out = zeros(N,7);

for i = 1:N
    x = X(i,1);  y = X(i,2);  z = X(i,3);
    if (x>=xmin) & (x<=xmax) & (y>=ymin) & (y<=ymax) & (z>=zmin) & (z<=zmax)
        R = sqrt((X(:,1)-x).*(X(:,1)-x) + (X(:,2)-y) .* (X(:,2)-y) + (X(:,3)-z) .* (X(:,3)-z));
        XX = [X R];
        XX(XX(:,4) < min,:) = [];
        XX(XX(:,4) > max,:) = [];
        XX(XX(:,4) == 0,:) = [];
        [NN,MM] = size(XX);
        if NN > 3
           nnvectors = zeros(NN,3);
           products = zeros(NN*NN,1);
           opposingneighbors = 0;
           for j = 1:NN
               nnvectors(j,:) = XX(j,1:3) - [x,y,z];
           end
           for j = 1:NN
                for k = (j+1):NN
                   products(j*NN+k,1) = dot(nnvectors(j,:),nnvectors(k,:)) / (norm(nnvectors(j,:)) * norm(nnvectors(k,:)));
              end
           end  
           products = sortrows(products,1);
           out(i,1) = NN;
           out(i,2:7) = (products(1:6,1))';
         else
            out(i,1:7) = [0,0,0,0,0,0,0];
         end;
    else
        out(i,1:7) = [0,0,0,0,0,0,0];
    end
end