Skip to content
Snippets Groups Projects
Commit 41905317 authored by Nicola Vigano's avatar Nicola Vigano
Browse files

Def/Shape-functions: small optimization

parent 32e73c8a
No related branches found
No related tags found
No related merge requests found
......@@ -71,19 +71,19 @@ function [om, plrot, plsigned, rot, omind] = gtFedPredictOmegaMultiple(...
%%% Quadratic equation for omega angles
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nn = size(pl,2);
nn = size(pl, 2);
if size(beamdir,1)==3 && size(beamdir,2)==1
A = beamdir'*(rotcomp.cos*pl);
B = beamdir'*(rotcomp.sin*pl);
C = beamdir'*(rotcomp.const*pl);
if ((size(beamdir, 1) == 3) && (size(beamdir, 2) == 1))
A = beamdir' * (rotcomp.cos * pl);
B = beamdir' * (rotcomp.sin * pl);
C = beamdir' * (rotcomp.const * pl);
else
A = sum(beamdir.*(rotcomp.cos*pl),1);
B = sum(beamdir.*(rotcomp.sin*pl),1);
C = sum(beamdir.*(rotcomp.const*pl),1);
A = sum(beamdir .* (rotcomp.cos * pl), 1);
B = sum(beamdir .* (rotcomp.sin * pl), 1);
C = sum(beamdir .* (rotcomp.const * pl), 1);
end
D = A.^2 + B.^2;
D = A .^ 2 + B .^ 2;
if (~isempty(omind))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
......@@ -102,7 +102,7 @@ if (~isempty(omind))
ss = ssp - ~ssp;
ss = ss .* ok;
om = 2*atand((-B + ss .* E) ./ (CC - A));
om = 2 * atand((-B + ss .* E) ./ (CC - A));
om(~ok) = NaN;
% Limit range
......@@ -138,8 +138,8 @@ else
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)
......@@ -148,8 +148,8 @@ else
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);
......@@ -179,9 +179,7 @@ else
rot = reshape(rot, 3, 3, [], 4);
plt = reshape(pl, 1, 3, []);
plt = plt([1 1 1], :, :);
plrot = rot .* cat(4, plt, plt, plt, plt);
plrot = bsxfun(@times, rot, plt);
plrot = sum(plrot, 2);
plrot = reshape(plrot, 3, nn, []);
......@@ -220,5 +218,4 @@ else
plrot(:, :, [3 4]) = -plrot(:, :, [3 4]);
end
end % of function
\ No newline at end of file
......@@ -29,7 +29,6 @@ num_def = size(defT, 3);
% Get an arbitrary vector e1 perpendicular to n:
[~, maxcomp] = max(abs(pl), [], 1);
ee = zeros(3, num_pls * 3);
zero = zeros(1, num_pls);
......@@ -41,13 +40,11 @@ ee(:, 3:3:end) = [-pl(3, :); zero; pl(1, :)];
ind2 = (0:num_pls-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], :);
......@@ -57,48 +54,47 @@ e2 = [pl(2, :) .* e1(3, :) - pl(3, :) .* e1(2, :);
%
% 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 (num_def == 1)
e1def = defT*e1 + e1;
e2def = defT*e2 + e2;
e1def = defT * e1 + e1;
e2def = defT * e2 + e2;
elseif (num_pls == 1)
defT_t = permute(defT, [2 1 3]);
defT_t = reshape(defT_t, 3, [])';
e1def = reshape(defT_t * e1, 3, []) + e1(:, ones(num_def, 1));
e2def = reshape(defT_t * e2, 3, []) + e2(:, ones(num_def, 1));
defT_v = reshape(defT_t, 3, []);
e1def = e1' * defT_v;
e1def = reshape(e1def, 3, []);
e1def = bsxfun(@plus, e1def, e1);
e2def = e2' * defT_v;
e2def = reshape(e2def, 3, []);
e2def = bsxfun(@plus, e2def, e2);
else
% Expand vectors for multiplication element-wise
e1r = reshape(e1, 1, []);
defTe = reshape(defT, 3, [], 1).*[e1r; e1r; e1r];
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 = 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, :)];
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], :);
normpldef = sqrt(sum(pldef .^ 2, 1));
pldef = bsxfun(@times, pldef, 1 ./ normpldef);
% check flip
% if sum(pldef.*pl, 1) < 0
......@@ -106,26 +102,25 @@ pldef = pldef./normpldef([1 1 1], :);
% 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;
fvec = defT * pl + pl;
elseif (num_pls == 1)
fvec = reshape(defT_t * pl, 3, []) + pl(:, ones(num_def, 1));
fvec = pl' * defT_v;
fvec = reshape(fvec, 3, []);
fvec = bsxfun(@plus, fvec, pl);
else
plr = reshape(pl, 1, []);
defTe = reshape(defT, 3, [], 1).*[plr; plr; plr];
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);
drel = sum(pldef .* fvec, 1);
end % of function
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment