gtStrainPlaneNormals.m

function [pldef, drel] = gtStrainPlaneNormals(pl, defT)
% GTSTRAINPLANENORMALS Given a deformation tensor, calculates new plane 
% normals with no approximations.
% 
%   [pldef, drel] = gtStrainPlaneNormals(pl, deft)
%
% Returns the new orientations of plane normals 'pl' after a deformation 
% described by the deformation tensor 'defT' (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 approximation is used.
% The 'pl' and 'defT' coordinates must be given in the same reference frame.
%
% INPUT
%   pl    - normalised coordinates of plane normals (3xn)
%   defT  - deformation tensor(s) (3x3 or 3x3xn)
%
% OUTPUT
%   pldef - plane normals in deformed state
%   drel  - relative elongations along the (deformed) plane normals
%


% 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


% e1 and e2 vectors in deformed state:
if (size(defT,3) == 1)
    e1def = defT*e1 + e1;
    e2def = defT*e2 + e2;
else
    % Expand vectors for multiplication element-wise
    e1r   = reshape(e1,1,[]);
    defTe = reshape(defT,3,[],1).*[e1r; e1r; e1r];
    defTe = defTe(:,1:3:end) + defTe(:,2:3:end) + defTe(:,3:3:end);
    e1def = defTe + e1;

    e2r   = reshape(e2,1,[]);
    defTe = reshape(defT,3,[],1).*[e2r; e2r; e2r];
    defTe = defTe(:,1:3:end) + defTe(:,2:3:end) + defTe(:,3:3:end);
    e2def = defTe + e2;
end


% 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 length along original plane normal
% (exact calculation, non-linear in the defT components)

if (size(defT,3) == 1)
    fvec = defT*pl + pl;
else
    plr   = reshape(pl,1,[]);
    defTe = reshape(defT,3,[],1).*[plr; plr; plr];
    defTe = defTe(:,1:3:end) + defTe(:,2:3:end) + defTe(:,3:3:end);
    
    fvec = defTe + pl;   
end


% new normalised d-spacing
drel = sum(pldef.*fvec, 1);


end % of function