gtCalculateGrain: restored old fashioned ordering of reflections!

It gives the exact same order of the omegas as the old non parallel gtCalculateGrain.

Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>

4_grains/gtCalculateGrain.m

@@ -158,15 +158,23 @@ pls4(:, todel, :)   = [];
 St4(:, :, todel, :) = [];
 pl(:, todel)        = [];
 
+plorig4 = repmat(pl, [1, 1, 4]);
+
 thetatypesp = cryst.thetatypesp(~todel);
 thetatypesp = reshape(thetatypesp, [], 1);
 hkl = cryst.hkl(:, thetatypesp)';
-hklsp  = cryst.hklsp(:, ~todel)';
+hklsp = cryst.hklsp(:, ~todel)';
+
+thetatypesp4 = repmat(thetatypesp, [1, 4]);
+hkl4 = repmat(hkl, [1, 1, 4]);
+hklsp4 = repmat(hklsp, [1, 1, 4]);
+hklsp4(:, :, [3 4]) = -hklsp4(:, :, [3 4]);
+sinth4 = repmat(sinth', [1 4])';
+theta4 = repmat(theta', [1 4])';
 
 %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % % Following the convention of the matching output, the omega value
 % % smaller than 180deg (spot "a") is used in the pair data.
-% % The two Friedel pairs are 1a-1b and 2a-2b.
 
 chind4 = om4(1, :) > om4(3, :);
 
@@ -185,14 +193,32 @@ omind4([2, 4], chind4)    = omind4([4, 2], chind4);
 St4(:, :, chind4, [2, 4]) = St4(:, :, chind4, [4, 2]);
 
 %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Four omegas of the plane normal (1-3 and 2-4 are the two Friedel pairs):
+% % But the two Friedel pairs are going to be 1a-1b and 2a-2b.
+
+pllab4(:, :, [2, 3]) = pllab4(:, :, [3, 2]);
+pls4(:, :, [2, 3])   = pls4(:, :, [3, 2]);
+om4([2, 3], :)       = om4([3, 2], :);
+omind4([2, 3], :)    = omind4([3, 2], :);
+St4(:, :, :, [2, 3]) = St4(:, :, :, [3, 2]);
+hklsp4(:, :, [2, 3]) = hklsp4(:, :, [3, 2]);
+
+%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+om = reshape(om4, [], 1)';
+omind = reshape(omind4, [], 1)';
+sinth = reshape(sinth4, [], 1);
+theta = reshape(theta4, [], 1);
+
+pllab = reshape(permute(pllab4, [1 3 2]), 3, []);
+pls = reshape(permute(pls4, [1 3 2]), 3, []);
+plorig = reshape(permute(plorig4, [1 3 2]), 3, []);
+
+St = reshape(permute(St4, [1 2 4 3]), 3, 3, []);
 
-om = reshape(om4', [], 1)';
-omind = reshape(omind4', [], 1)';
-sinth = reshape(sinth, [], 1);
-theta = reshape(theta, [], 1);
-pllab = reshape(pllab4, 3, []);
-pls = reshape(pls4, 3, []);
-St = reshape(St4, 3, 3, []);
+thetatypesp = reshape(thetatypesp4', [], 1);
+hkl = reshape(permute(hkl4, [3 1 2]), [], 3);
+hklsp = reshape(permute(hklsp4, [3 1 2]), [], 3);
 
 % Diffraction vectors % takes coloumn vectors
 dveclab = gtFedPredictDiffVecMultiple(pllab, beamdir);
@@ -214,14 +240,14 @@ uvw = gtFedPredictUVWMultiple(St, dveclab, csam_v, detpos, detnorm, Qdet, uvorig
 % Initialse output variables
 grain.allblobs.uvw = uvw';
 grain.allblobs.uv = uvw(1:2, :)';
-grain.allblobs.plorig = [pl'; pl'; pl'; pl'];
+grain.allblobs.plorig = plorig';
 grain.allblobs.pl = pls';
 grain.allblobs.pllab = pllab';
-grain.allblobs.hkl = [hkl; hkl; hkl; hkl];
-grain.allblobs.hklsp = [hklsp; -hklsp; hklsp; -hklsp];
-grain.allblobs.sintheta = [sinth; sinth; sinth; sinth];
-grain.allblobs.theta = [theta; theta; theta; theta];
-grain.allblobs.thetatype = [thetatypesp; thetatypesp; thetatypesp; thetatypesp];
+grain.allblobs.hkl = hkl;
+grain.allblobs.hklsp = hklsp;
+grain.allblobs.sintheta = sinth;
+grain.allblobs.theta = theta;
+grain.allblobs.thetatype = thetatypesp;
 grain.allblobs.omega = om';
 grain.allblobs.omind = omind';
 grain.allblobs.dvec = dveclab';

7_fed/geometry/gtFedPredictOmega.m

-function [om, pllab, plsigned, Srot] = gtFedPredictOmega(pl, sinth, beamdir, rotdir, rotcomp)
+function [om, pllab, plsigned, Srot, omind] = gtFedPredictOmega(pl, sinth, beamdir, rotdir, rotcomp)
 %    FUNCTION  [om,pllab,plsigned, Srot] = gtFedPredictOmega(pl,sinth,beamdir,rotdir,rotcomp)
 %
 % Predicts the four omega angles where diffraction occurs for a plane normal pl
@@ -75,14 +75,17 @@ if all(isnan(om))
     pllab = [];
     plsigned = [];
     Srot = [];
+    omind = [];
     return
 end
 
 om = mod(om, 360);
+omind = [1; 2; 3; 4];
 
 % Exchange om1 and om2 if om1 > om2
 if (om(1) > om(2))
     om(1:2) = [om(2) om(1)];
+    omind(1:2) = [omind(2) omind(1)];
 end
 
 % Rotation matrix for each omega
@@ -104,6 +107,8 @@ or3 = tempCrossProd * pllab(:, 3);
 if ((or1 * or3) > 0) % orientation is not opposite, exchange om3 to om4
     om(3:4) = [om(4) om(3)];
 
+    omind(3:4) = [omind(4) omind(3)];
+
     Srot(:, :, [3 4]) = Srot(:, :, [4 3]);
 
     pllab(:, 3:4) = [pllab(:, 4) pllab(:, 3)];

7_fed2/gtFedPredictOmegaMultiple.m

@@ -133,111 +133,102 @@ if ~isempty(omind)
     plsigned = ss3.*pl;
     
 else
-    
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     %% All four omega indices
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    
+
     om = NaN(4,nn);
-    
+
     % Diffraction condition 1
     % Pl_rot and beamdir dotproduct equals -sin(th)
     CC = C + sinth;
     DD = D - CC.^2;
     ok = DD > 0;
     E  = sqrt(DD(ok));
-    
-    om(1,ok) = 2*atand((-B(ok)-E)./(CC(ok)-A(ok)));
-    om(2,ok) = 2*atand((-B(ok)+E)./(CC(ok)-A(ok)));
-    
-    
+
+    om(1, ok) = 2*atand((-B(ok)-E)./(CC(ok)-A(ok)));
+    om(2, ok) = 2*atand((-B(ok)+E)./(CC(ok)-A(ok)));
+
     % Diffraction condition 2
     % Pl_rot and beamdir dotproduct equals +sin(th)
     CC = C - sinth;
     DD = D - CC.^2;
     ok = DD > 0;
     E  = sqrt(DD(ok));
-    
-    om(3,ok) = 2*atand((-B(ok)+E)./(CC(ok)-A(ok)));
-    om(4,ok) = 2*atand((-B(ok)-E)./(CC(ok)-A(ok)));
-    
-    
+
+    om(3, ok) = 2*atand((-B(ok)+E)./(CC(ok)-A(ok)));
+    om(4, ok) = 2*atand((-B(ok)-E)./(CC(ok)-A(ok)));
+
     % Limit range
-    om = mod(om,360);
-    
+    om = mod(om, 360);
+
     % Omegea indices
     omind = [1; 2; 3; 4];
-    omind = omind(:,ones(1,nn));
+    omind = omind(:, ones(1, nn));
 
-    
     % For output conventions:
-    % make sure om1 is smaller then om2
-    chom = om(1,:) > om(2,:);
-    om(1:2,chom)    = [om(2,chom); om(1,chom)];
-    omind(1:2,chom) = [omind(2,chom); omind(1,chom)];
-   
-    
+    % make sure om1 is smaller than om2
+    chom = om(1, :) > om(2, :);
+
+    om([1 2], chom)    = om([2 1], chom);
+    omind([1 2], chom) = omind([2 1], chom);
+
     % ROTATION MATRICES AND PLANE NORMALS IN LAB
 
-    plsigned        = zeros(3,nn,4);
-    plsigned(:,:,1) = pl;
-    plsigned(:,:,2) = pl;
-    plsigned(:,:,3) = -pl;
-    plsigned(:,:,4) = -pl;
-    
+    plsigned = zeros(3, nn, 4);
+    plsigned(:, :, 1) = pl;
+    plsigned(:, :, 2) = pl;
+    plsigned(:, :, 3) = -pl;
+    plsigned(:, :, 4) = -pl;
+
     % Get rotation matrices and multpily the input plane normals to get
     % them in the diffracting position (avoid looping)
-    rot(:,:,:,4) = gtMathsRotationTensor(om(4,:), rotcomp);
-    rot(:,:,:,3) = gtMathsRotationTensor(om(3,:), rotcomp);
-    rot(:,:,:,2) = gtMathsRotationTensor(om(2,:), rotcomp);
-    rot(:,:,:,1) = gtMathsRotationTensor(om(1,:), rotcomp);
-    
-    
-    plt = reshape(pl,1,[]);
-    plt = plt([1 1 1],:);
-    
-    plrot = reshape(rot,3,[]).*[plt, plt, plt, plt];
-    
-    plrot = plrot(:,1:3:end) + plrot(:,2:3:end) + plrot(:,3:3:end);
-    
-    plrot = reshape(plrot,3,nn,[]);
-      
-    
+    rot(:, :, :, 4) = gtMathsRotationTensor(om(4, :), rotcomp);
+    rot(:, :, :, 3) = gtMathsRotationTensor(om(3, :), rotcomp);
+    rot(:, :, :, 2) = gtMathsRotationTensor(om(2, :), rotcomp);
+    rot(:, :, :, 1) = gtMathsRotationTensor(om(1, :), rotcomp);
+
+    plt = reshape(pl, 1, []);
+    plt = plt([1 1 1], :);
+
+    plrot = reshape(rot, 3, []) .* [plt, plt, plt, plt];
+
+    plrot = plrot(:, 1:3:end) + plrot(:, 2:3:end) + plrot(:, 3:3:end);
+
+    plrot = reshape(plrot, 3, nn, []);
+
     % Exchange om3 and om4 if needed, so that om1-om3 and om2-om4 are the
     % opposite reflections (Friedel pairs in case the rotation axis is
     % perpendicular to the beam).
     % Midplane of setup is defined by beam direction and rotation axis.
     % This is needed to consider consistently all possible setups where 
     % the beam and rotation axis are not perpendicular.
-    
+
     % Cross product beamdir and rotdir
-    br = [beamdir(2,:)*rotdir(3) - beamdir(3,:)*rotdir(2); ...
-          beamdir(3,:)*rotdir(1) - beamdir(1,:)*rotdir(3); ...
-          beamdir(1,:)*rotdir(2) - beamdir(2,:)*rotdir(1)];
-    
+    br = [ beamdir(2, :) * rotdir(3) - beamdir(3, :) * rotdir(2); ...
+           beamdir(3, :) * rotdir(1) - beamdir(1, :) * rotdir(3); ...
+           beamdir(1, :) * rotdir(2) - beamdir(2, :) * rotdir(1) ];
+
     % Dot product of br and pllab
-    if size(beamdir,1)==3 && size(beamdir,2)==1
-        dot1 = br'*plrot(:,:,1);
-        dot3 = br'*plrot(:,:,3);
+    if (size(beamdir, 1) == 3) && (size(beamdir, 2) == 1)
+        dot1 = br' * plrot(:, :, 1);
+        dot3 = br' * plrot(:, :, 3);
     else
-        dot1 = sum(br.*plrot(:,:,1),1);
-        dot3 = sum(br.*plrot(:,:,3),1);
+        dot1 = sum(br .* plrot(:, :, 1), 1);
+        dot3 = sum(br .* plrot(:, :, 3), 1);
     end
-    
+
     % Indices to be swapped
-    chom = (((dot1.*dot3) > 0) & ~isnan(dot1) & ~isnan(dot3)) ;
-    
+    chom = (((dot1 .* dot3) > 0) & ~isnan(dot1) & ~isnan(dot3)) ;
+
     % Swap 3rd and 4th
-    om([3,4],chom)      = om([4,3],chom);
-    omind([3,4],chom)   = omind([4,3],chom);
-    plrot(:,chom,[3,4]) = plrot(:,chom,[4,3]);
-    rot(:,:,chom,[3,4]) = rot(:,:,chom,[4,3]);
-    
-    
+    om([3 4], chom)        = om([4 3], chom);
+    omind([3 4], chom)     = omind([4 3], chom);
+    plrot(:, chom, [3 4])  = plrot(:, chom, [4 3]);
+    rot(:, :, chom, [3 4]) = rot(:, :, chom, [4 3]);
+
     % Change all pllab to opposite direction
-    plrot(:,:,[3,4]) = -plrot(:,:,[3,4]);
-    
-    
+    plrot(:, :, [3 4]) = -plrot(:, :, [3 4]);
 end