function out = averagestrainfield(X,epsilon)

% This function smoothes the local strain field by averaging the epsilon tensor 
% of the particles with that of its nearest neighbors;
% X is a N by 3 list of x,y,z coordinates of all N particles, epsilon is a
% N by 9 list of the elements of the strain tensor of the particles,
% calculated using the "strainfield.m" routine.
% You can change the range of the averaging by changing the NNacceptance value below.
% Output is a N by 9 list of the elements of the averaged strain tensor.
%
% Peter Schall, 2004
%
%
N = size(X,1);
d0 = 1/sqrt(2);
NNacceptance = 1.2 * d0;
averagestrain = zeros(N,9);

for n0 = 1:N
    if epsilon(n0,1) ~= 0
        p1 = X(n0,:);
        r = sqrt((X(:,1)-p1(1)).^2 + (X(:,2)-p1(2)).^2 + (X(:,3)-p1(3)).^2);
        x1 = [X,epsilon,r];
        x1(x1(:,13) > NNacceptance,:) = [];
        x1(x1(:,4) == 0,:) = [];
        NN = size(x1,1);
        localepsilon = x1(:,4:12);
        averagestrain(n0,:) = sum(localepsilon) / NN;
    end
end
out = averagestrain;