diff --git a/zUtil_Strain2/gtStrainGrainsGeneral.m b/zUtil_Strain2/gtStrainGrainsGeneral.m
new file mode 100644
index 0000000000000000000000000000000000000000..23bcd0340630e0e9c6e5eb63cdac93c59a2a17c8
--- /dev/null
+++ b/zUtil_Strain2/gtStrainGrainsGeneral.m
@@ -0,0 +1,89 @@
+function grain = gtStrainGrainsGeneral(grain, cryst, energy, Bmat, stdlim)
+% GTSTRAINGRAINSGENERAL grain = gtStrainGrainsGeneral(grain, cryst, energy,
+% Bmat, stdlim)
+%
+% Runs a strain calculation for all the grains in the cell input 'grain',
+% and puts the results in the field grain{i}.strain.
+% Needs the 'cryst' data as in the output of Indexter
+% (in file 4_grains/phase_##/index.mat) and the beam energy in keV.
+% It filters outliers based on their principal strains being away from the
+% median value of all grains more than 'stdlim' times its std.
+% The fields 'Bmat' and 'stdlim' are optional.
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Set input
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+tic
+
+lambda = gtConvEnergyToWavelength(energy);
+
+if ~exist('Bmat','var') || isempty(Bmat)
+ Bmat = gtCrystHKL2CartesianMatrix(cryst.latticepar);
+end
+
+if ~exist('stdlim','var') || isempty(stdlim)
+ stdlim = Inf;
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Loop through grains
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+parfor ii = 1:length(grain)
+
+ % Deformation tensor
+ DT = gtStrainSingleGrain(grain{ii}, cryst, lambda, Bmat);
+ grain{ii}.strain.DT = DT;
+
+ % Strain tensor - symmetric part
+ grain{ii}.strain.ST = 0.5*(DT + DT');
+
+ % Rotation tensor - skew part
+ grain{ii}.strain.ROT = 0.5*(DT - DT');
+
+
+ % Principal strains - eigenvalues and eigenvectors
+ [ivec, ival] = eig(grain{ii}.strain.ST);
+
+ [ivalsort, ivalinds] = sort([ival(1,1), ival(2,2), ival(3,3)], 2, 'descend');
+
+ % principal strain values ordered
+ grain{ii}.strain.pstrain = ivalsort;
+
+ % principal strain directions in columns
+ grain{ii}.strain.pstraindir = ivec(:, ivalinds);
+
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Filter outliers
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+med_ps = median(gtIndexAllGrainValues(grain,'strain','pstrain',1));
+std_ps = std(gtIndexAllGrainValues(grain,'strain','pstrain',1));
+
+
+parfor ii = 1:length(grain)
+
+ if (grain{ii}.strain.pstrain(1) > med_ps + stdlim*std_ps) || ...
+ (grain{ii}.strain.pstrain(1) < med_ps - stdlim*std_ps)
+
+ grain{ii}.strain.DT(:) = NaN;
+ grain{ii}.strain.ST(:) = NaN;
+ grain{ii}.strain.ROT(:) = NaN;
+ grain{ii}.strain.pstrain(:) = NaN;
+ grain{ii}.strain.pstraindir(:) = NaN;
+
+ end
+
+end
+
+
+toc
+
+
+end % of function
\ No newline at end of file
diff --git a/zUtil_Strain2/gtStrainPlaneNormals.m b/zUtil_Strain2/gtStrainPlaneNormals.m
new file mode 100644
index 0000000000000000000000000000000000000000..47f87ae257290a5e38f1d0d7af35821dfcc35ffe
--- /dev/null
+++ b/zUtil_Strain2/gtStrainPlaneNormals.m
@@ -0,0 +1,86 @@
+function [pldef, drel] = gtStrainPlaneNormals(pl, G)
+% FUNCTION [pldef, drel] = gtStrainPlaneNormals(pl, G)
+%
+% Returns the new orientations of plane normals 'pl' after a deformation
+% described by the deformation tensor 'G' (9 distinct components).
+% Also returns the elongations 'drel' due to the deformation in the
+% directions of the plane normals. In a crystal this equals to the ratios
+% of the d-spacings after and before the deformation.
+% The results are exact, no linear approximation is used.
+% The 'pl' and 'G' coordinates must be given in the same reference frame.
+%
+
+
+% Get an arbitrary vector e1 perpendicular to n:
+[~, maxcomp] = max(abs(pl), [], 1);
+
+
+ee = zeros(3,size(pl,2)*3);
+
+zero = zeros(1,size(pl,2));
+
+ee(:,1:3:end) = [ pl(2,:); -pl(1,:); zero];
+ee(:,2:3:end) = [ zero; pl(3,:); -pl(2,:)];
+ee(:,3:3:end) = [-pl(3,:); zero; pl(1,:)];
+
+ind2 = (0:size(pl,2)-1)*3 + maxcomp;
+e1 = ee(:,ind2);
+
+
+% Cross product of e1 and pl: e2=cross(e1,pl)
+e2 = [pl(2,:).*e1(3,:) - pl(3,:).*e1(2,:);
+ pl(3,:).*e1(1,:) - pl(1,:).*e1(3,:);
+ pl(1,:).*e1(2,:) - pl(2,:).*e1(1,:)];
+
+
+% Check
+% ne2 = sqrt(sum(e2.*e2, 1));
+% e2n = e2./ne2([1 1 1],:);
+%
+% ne1 = sqrt(sum(e1.*e1, 1));
+% e1n = e1./ne1([1 1 1],:);
+%
+% cross(e1n,e2n) - pl
+
+
+% check flip e2 if needed to get them right handed:
+% ch = [e1(2,:).*e2(3,:) - e1(3,:).*e2(2,:);
+% e1(3,:).*e2(1,:) - e1(1,:).*e2(3,:);
+% e1(1,:).*e2(2,:) - e1(2,:).*e2(1,:)];
+% sum(ch.*pl, 1) < 0
+
+
+% e-s in deformed state:
+e1def = G*e1 + e1;
+e2def = G*e2 + e2;
+
+
+% new deformed pl
+pldef = [e1def(2,:).*e2def(3,:) - e1def(3,:).*e2def(2,:);
+ e1def(3,:).*e2def(1,:) - e1def(1,:).*e2def(3,:);
+ e1def(1,:).*e2def(2,:) - e1def(2,:).*e2def(1,:)];
+
+
+% normalise pldef
+normpldef = sqrt(sum(pldef.*pldef, 1));
+pldef = pldef./normpldef([1 1 1],:);
+
+
+% check flip
+% if sum(pldef.*pl, 1) < 0
+% pldef = ...-pldef...;
+% end
+%sum(pldef.*pl, 1) >= 0
+
+
+
+% new relative d-spacing (exact calculation, non-linear in the G components)
+fvec = pl + G*pl; % change along original plane normal
+
+% new normalised d-spacing
+drel = sum(pldef.*fvec, 1);
+
+
+
+end % of function
+
diff --git a/zUtil_Strain2/gtStrainSingleGrain.m b/zUtil_Strain2/gtStrainSingleGrain.m
new file mode 100644
index 0000000000000000000000000000000000000000..0243d67cc97e5ae2b6200a52af49ad2adde20f35
--- /dev/null
+++ b/zUtil_Strain2/gtStrainSingleGrain.m
@@ -0,0 +1,146 @@
+function DT = gtStrainSingleGrain(grain, cryst, lambda, Bmat)
+% GTSTRAINSINGLEGRAIN DT = gtStrainSingleGrain(grain, cryst, lambda, Bmat)
+%
+% Computes the deformation tensor 'DT' for a single indexed grain
+% (crystal system is irrelevant). The 9 deformation matrix components are
+% computed in the same reference as the plane normals.
+% The deformation components are found by fitting all the interplanar
+% angles and distances observed (no weighting of the equations). No linear
+% approximation is used.
+% If 'Bmat' is empty, it will be calculated inside the function.
+%
+
+
+% All measured signed hkl-s in the grain
+shkl = [];
+for ii = 1:length(grain.shklind)
+ shkl(:,ii) = cryst.shklfam{grain.hklind(ii)}(grain.shklind(ii),:)';
+end
+
+
+% Measured d-spacings
+dsp = 0.5*lambda./sind(grain.theta);
+
+% D-spacings in the undeformed reference crystal
+dsp_ref = cryst.dspacing(grain.hklind);
+
+% Required relative elongations to achieve unstrained state from deformed state
+elo_ref = dsp_ref./dsp;
+
+
+% Get B matrix for hkl -> Cartesian coordinate transforms
+if isempty(Bmat)
+ Bmat = gtCrystHKL2CartesianMatrix(cryst.latticepar);
+end
+
+% Cartesian coordinates of plane normals in the unstrained reference crystal
+% (output plane normals are normalized)
+pl_ref = gtCrystHKL2Cartesian(shkl, Bmat);
+
+% All combinations of plane normal dot products in the unstrained reference crystal
+dotpro_ref = abs(pl_ref'*pl_ref);
+dotpro_ref(dotpro_ref>1) = 1; % correction for numerical errors
+
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Optimization
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Initial guess values for fitting the 9 components
+scguess = [0 0 0 0 0 0 0 0 0];
+
+% Lower and upper bounds for the components solution
+boundslow = -0.1*ones(1,9);
+boundshigh = 0.1*ones(1,9);
+
+
+% Optimization parameters
+par.tolx = 1e-6; % precision for components (1e-6 is reasonable and fast)
+par.tolfun = 1e-12; % optimization function target value (should be very low)
+par.typicalx = 1e-5*ones(1,9); % typical component values for calculating gradients
+
+
+
+% Find the "undeformation tensor" components 'sc' that would bring the
+% observed, deformed grain into the undeformed reference state.
+% Using Matlab's built-in non-linear least square solver.
+scopt = lsqnonlin( @(sc) sfDeviation(sc, grain.pl, elo_ref, dotpro_ref),...
+ scguess, boundslow, boundshigh,...
+ optimset('TolX',par.tolx, 'TolFun',par.tolfun,...
+ 'TypicalX',par.typicalx, 'FunValCheck','off', 'Display','off'));
+
+%[scopt, resnorm, residual, exitflag, optimout, lagrange, jacobian] = ...
+% lsqnonlin( @(sc) sfDeviation(sc, grain.pl, elo_ref, dotpro_ref),...
+% scguess, boundslow, boundshigh,...
+% optimset('TolX',par.tolx, 'TolFun',par.tolfun,...
+% 'TypicalX',par.typicalx, 'FunValCheck','off', 'Display','off'));
+
+
+% "Undeformation tensor"
+DTundef = reshape(scopt,3,3);
+
+I = [1 0 0; 0 1 0; 0 0 1];
+
+% The real deformation tensor is a kind of inverse of DTundef
+DT = inv(I + DTundef) - I;
+
+
+end % of main function
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Sub-function
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Gives a vector of the deviations in interplanar angles and distances
+% as a function of the deformation components 'sc'. These differences
+% should be eliminated by finding the right 'sc'.
+
+
+function dev = sfDeviation(sc, pl0, elo_ref, dotpro_ref)
+
+ DT = reshape(sc,3,3);
+
+ % 'pl0' is the as observed plane normals; 'pl' is 'pl0' after
+ % deformation 'DT'
+ [pl, elo] = gtStrainPlaneNormals(pl0, DT);
+
+
+ % Difference in elongation between computed and reference
+ elodiff = elo - elo_ref;
+
+
+ % Difference in dot products between computed and reference
+ dotprodiff = abs(pl'*pl) - dotpro_ref;
+
+ % Take upper triangle of matrix
+ inds = triu(true(size(dotprodiff)),1);
+ dotprodiff = dotprodiff(inds);
+
+
+ % All differences in one vector
+ dev = [elodiff(:); dotprodiff];
+
+end
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+