diff --git a/4_grains/gtCalculateGrain.m b/4_grains/gtCalculateGrain.m
index 9acfe56eba3ad74ccfbfba3a7e57bb601e5e2946..91d6ff2ed24f34c37abdafa982cd0ad4165ea19b 100644
--- a/4_grains/gtCalculateGrain.m
+++ b/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';
diff --git a/7_fed/geometry/gtFedPredictOmega.m b/7_fed/geometry/gtFedPredictOmega.m
index 4d8861586ac24f6589df9e6f20118aecefb3e11b..ee1b2461c034966b95b295373433c25e91c26011 100644
--- a/7_fed/geometry/gtFedPredictOmega.m
+++ b/7_fed/geometry/gtFedPredictOmega.m
@@ -1,4 +1,4 @@
-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)];
diff --git a/7_fed2/gtFedPredictOmegaMultiple.m b/7_fed2/gtFedPredictOmegaMultiple.m
index 9a3f2a91c30201a4541f4a03718a3779d1cd6fbf..f7e7127035dde80f1cf5cbe22b9e58e04d30bbb7 100644
--- a/7_fed2/gtFedPredictOmegaMultiple.m
+++ b/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