function out = strainfield(X)

% This function calculates the local strain field of a distorted fcc crystal lattice.
% X is an N by 3 list of the x,y,z coordinates of all N particles of the
% crystal.
% Use the vector "field" below to set the region where the strain field should be determined (xmin, xmax, ymin, ... )
% The locality of the strain field calculation can be changed by changing
% the NNacceptance value below. 
% The program outputs a 9 by N list containing all the elements of the
% strain tensor (in the order: epsxx, epsxy, epsxz, epsyx, epsyy, ...)
% Additional smoothing of the strain distribution
% may be necessary; this can be achieved using the function averagestrainfield.
%
% The strain field calculation of this routine is based on an algorithm proposed by
% M.L.Falk and J.S.Langer, Phys.Rev.E, 57,7192 (1998)
%
%
% Peter Schall and Dan Bonner, 2004;
%
%
xmin = 0; xmax = 10;
ymin = 0; ymax = 10;
zmin = 0; zmax = 10;

d0 = 1/sqrt(2);
N = size(X,1);
correct_zcompression = 0;
NNacceptance = 1.2 * d0;

%all 12 neighbors in fcc lattice
fcc = zeros(12:3);
fcc(1,:)=[d0,0,0];  fcc(2,:)=[0,d0,0];  fcc(3,:)=[-d0,0,0]; fcc(4,:)=[0,-d0,0];
fcc(5,:)=[d0/2,d0/2,d0/sqrt(2)];      fcc(6,:)=[d0/2,d0/2,-d0/sqrt(2)];
fcc(7,:)=[-d0/2,d0/2,d0/sqrt(2)];     fcc(8,:)=[-d0/2,d0/2,-d0/sqrt(2)];
fcc(9,:)=[d0/2,-d0/2,d0/sqrt(2)];     fcc(10,:)=[d0/2,-d0/2,-d0/sqrt(2)];
fcc(11,:)=[-d0/2,-d0/2,d0/sqrt(2)];   fcc(12,:)=[-d0/2,-d0/2,-d0/sqrt(2)];

strainepsilon = zeros(N,9);
currentfcc = fcc;

for n0 = 1:N
    FalkX = zeros(3,3);
    FalkY = zeros(3,3);
    deform = zeros(3,3);
    if X(n0,1)>xmin & X(n0,1)<xmax & X(n0,2)>ymin & X(n0,2)<ymax & X(n0,3)>zmin & X(n0,3)<zmax
        p1 = X(n0,:);
        if correct_zcompression == 1
            currentfcc(:,3) = fcc(:,3) * (1 + (0.00105*(p1(3)-6.5)));
        end
        r = sqrt((X(:,1)-p1(1)).^2 + (X(:,2)-p1(2)).^2 + (X(:,3)-p1(3)).^2);
        x1 = [X,r];
        x1(x1(:,4) > NNacceptance,:) = [];
        x1(x1(:,4) == 0,:) = [];
        NN = size(x1,1);
        Rn0 = zeros(NN,3);
        Rn0(:,1) = x1(:,1) - p1(1);  Rn0(:,2) = x1(:,2) - p1(2);  Rn0(:,3) = x1(:,3) - p1(3);
        Rn1 = zeros(NN,3);
        for n = 1:NN
            d = sqrt((currentfcc(:,1)-Rn0(n,1)).^2 + (currentfcc(:,2)-Rn0(n,2)).^2 + (currentfcc(:,3)-Rn0(n,3)).^2);
            [dummy,bestapproach] = min(d);
            Rn1(n,:) = currentfcc(bestapproach,:);
        end
        % get Falk's X ------------------------------
        for i = 1:3
            for j = 1:3
                for n = 1:NN
                    FalkX(i,j) = FalkX(i,j) + (Rn0(n,i) * Rn1(n,j));
                end
            end
         end
         %get Falk's Y ----------------------------------
        for i = 1:3
            for j = 1:3
                for n = 1:NN
                    FalkY(i,j) = FalkY(i,j) + (Rn1(n,i) * Rn1(n,j));
                end
            end
        end
        Yinv = inv(FalkY);
        % Calculate affine deformation -------------------------------
        for i = 1:3
            for j = 1:3
                for k = 1:3
                    deform(i,j) = deform(i,j) + FalkX(i,k) * Yinv(j,k);
                end
            end
        end
        deform = deform - eye(3);
    % calculate the strain (=symmetric component)
        deform = 1/2 * (deform + deform');
        strainepsilon(n0,:) = deform(1:9);
    else
        strainepsilon(n0,:) = 1000;
    end
end   
out = strainepsilon;