Skip to content
Snippets Groups Projects
Commit 94af1a64 authored by Wolfgang Ludwig's avatar Wolfgang Ludwig
Browse files

Take into account initial gonio offset positions 


Signed-off-by: default avatarWolfgang Ludwig <wolfgang.ludwig@esrf.fr>
parent ceeb3be9
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
zUtil_Topotomo/calc_pole_tilts_id11_index2.m
+ 166
− 73
View file @ 94af1a64
function [phx, phy] = calc_pole_tilts_id11_index(parameters, gid, thetatypes)
% assumes initial tilts of the goniometer to be at 0, 0 - uses ID11
% goniometer geometry with maximum tilt of 20 deg for samrx and 10 deg for samry
load(sprintf('%s/4_grains/phase_01/index.mat',parameters.acq.dir));
maxvalue = tand(20);
% Sample reference system
LabX = [1 0 0]'; samgeo.dirx = [1 0 0];
LabY = [0 1 0]'; samgeo.diry = [0 1 0];
LabZ = [0 0 1]'; samgeo.dirz = [0 0 1];
% Directions of instrument axis at sample rotation omega = 0 (...due to sample omega = 0 corresponds to instrument rotation = -90 deg)
instrgeo.dirx = [ 0 -1 0];
instrgeo.diry = [ 1 0 0];
instrgeo.dirz = [ 0 0 1];
% Rotation Matrix (based on angle & axis -> Rodrigues Rotation Formula) v_rotated = R(angle, axis) * v
R = @(angle,axis)gtMathsRotationTensor(angle, gtMathsRotationMatrixComp(axis, 'col'));
% Loop through grains and their reflections...
for i = 1:length(gid)
if exist(sprintf('%s/4_grains/phase_01/grain_%04d.mat',parameters.acq.dir, gid(i)), 'file')
g=load(sprintf('%s/4_grains/phase_01/grain_%04d.mat',parameters.acq.dir, gid(i)));
else
g=gtCalculateGrain(grain{gid(i)},parameters);
end
ind = find(g.allblobs.pl(:,3) > maxvalue & ismember(g.allblobs.thetatype, thetatypes));
if ~isempty(ind)
[pl, ia, ic] = unique(g.allblobs.pl(ind,:), 'rows', 'stable');
ind = ind(ia);
for j = 1:length(ind)
hkl = g.allblobs.hkl(ind(j), :);
theta = g.allblobs.theta(ind(j),:);
G_sam = g.allblobs.pl(ind(j),:);
diffrz = g.allblobs.omega(ind(j))-90;
uv = g.allblobs.detector.uvw(ind(j),1:2);
% GInstr = Plane normal in Instrument reference system:
% samrx is upper tilt and rotates around LabX at diffrz = 0 (omega = 90)
% ramry is lower tile and rotates around LabY at diffrz = 0 (omega = 90)
G_instr = gtGeoSam2Sam(G_sam, samgeo, instrgeo, 1, 1);
% first solve for upper tilt: RX * G_sam = [r 0 s] -> G_sam(2)*cos(phx) - G_sam(3)*sin(phx) = 0
phx = atand(G_instr(2)/G_instr(3));
tmp = R(phx, LabX) * G_instr';
% now solve for lower tilt: RY * tmp = [0 0 1], RY = [cos(phy) 0 sin(phy); 0 1 0; -sin(phy) 0 cos(phy)]; ->
% phy = atand(-tmp(1)/tmp(3));
phy = atand(-tmp(1)/tmp(3));
% testing
% G_aligned = R(phy, LabY) * tmp;
if (abs(phy) < 10 && abs(phx) <20)
disp(sprintf('Grain %d: Found %d %d %d reflection: diffry = %f samrx = %f samry = %f\n', gid(i), hkl(1), hkl(2), hkl(3), theta, phx, phy));
disp(sprintf('Positions in DCT scan: diffrz = %f, u = %f v = %f\n', diffrz, uv(1), uv(2)));
end
end
end
end
function out = calc_pole_tilts_id11_index2(parameters, gid, thetatypes)
% assumes initial tilts of the goniometer to be at 0, 0 - uses ID11
% goniometer geometry with maximum tilt of 20 deg for samrx and 10 deg for samry
% To do: - automatically determine offsets
% - use Diffractometer class
cd(parameters.acq.dir)
conf = gtGetOffsets(parameters);
conf.samrx_max = 18;
conf.samry_max = 14;
conf.thetatype_slipplane = 1;
conf.test_results = 0;
conf.nfdtx_offset = 95.77;
conf.d3tz_offset = 0;
conf.distance = 16.43;
conf.omegas = [0 : 4: 356];
conf.loaddir = [0 0 1];
load(sprintf('%s/4_grains/phase_01/index.mat',parameters.acq.dir));
maxvalue = tand(20);
%% Rotation Matrix (based on angle & axis -> Rodrigues Rotation Formula) v_rotated = R(angle, axis) * v
R = @(angle,axis)gtMathsRotationTensor(angle, gtMathsRotationMatrixComp(axis, 'col'));
%% Sample reference system
samgeo = parameters.samgeo;
LabX = [1 0 0]';
LabY = [0 1 0]';
LabZ = [0 0 1]';
%% Directions of instrument axis at "dct" rotation angle omega = 0 (diffrz rotation angle of first image in scan)
T = R(conf.diffrz_offset, LabZ) * R(conf.samry_offset, LabY) * R(conf.samrx_offset, LabX);
instrgeo.orig = [0, 0, 0];
instrgeo.dirx = (T * LabX)';
instrgeo.diry = (T * LabY)';
instrgeo.dirz = (T * LabZ)';
instrgeo.voxsize = [1 1 1];
RZ0 = R(conf.diffrz_offset, LabZ);
%% Loop through grains and their reflections...
for i = 1:length(gid)
tt_id = 1;
tt_pars = {};
gr_cen_instr = gtGeoSam2Sam(grain{gid(i)}.center, samgeo, instrgeo, false, false);
if exist(sprintf('%s/4_grains/phase_01/grain_%04d.mat',parameters.acq.dir, gid(i)), 'file')
g=load(sprintf('%s/4_grains/phase_01/grain_%04d.mat',parameters.acq.dir, gid(i)));
else
g=gtCalculateGrain(grain{gid(i)},parameters);
end
ind = find(g.allblobs.pl(:,3) > maxvalue & ismember(g.allblobs.thetatype, thetatypes));
if ~isempty(conf.thetatype_slipplane)
slip_ind = find(g.allblobs.thetatype == conf.thetatype_slipplane & g.allblobs.omind == 1);
slip_planes = g.allblobs.pl(slip_ind, :);
end
if ~isempty(ind)
[pl, ia, ic] = unique(g.allblobs.pl(ind, :), 'rows', 'stable');
ind = ind(ia);
for j = 1:length(ind)
hklsp = g.allblobs.hklsp(ind(j), :);
hkl = g.allblobs.hkl(ind(j), :);
theta = g.allblobs.theta(ind(j), :);
G_sam = g.allblobs.pl(ind(j), :);
omega = g.allblobs.omega(ind(j));
uv = g.allblobs.detector.uvw(ind(j), 1:2);
%% GInstr = Plane normal in Instrument reference system:
% samrx is upper tilt and rotates around LabX at diffrz = 0
% samrx range: [-20 20]
% samry is lower tilt and rotates around LabY at diffrz = 0
% samry range: [-15 15]
%G_instr = gtGeoSam2Sam(G_sam, samgeo, instrgeo, 1, 1);
%G_instr = R(-conf.offset, LabZ) * G_sam'; G_instr = G_instr';
G_instr = T' * G_sam';
% first solve for upper tilt: RX * G_sam = [r 0 s] -> G_sam(2)*cos(phx) - G_sam(3)*sin(phx) = 0
phx = atand(G_instr(2)/G_instr(3));
tmp = R(phx, LabX) * G_instr;
% now solve for lower tilt: RY * tmp = [0 0 1], RY = [cos(phy) 0 sin(phy); 0 1 0; -sin(phy) 0 cos(phy)]; ->
% phy = atand(-tmp(1)/tmp(3));
phy = atand(-tmp(1)/tmp(3));
if conf.test_results
fprintf('samrx = %f samry = %f for %s\n',phx, phy, num2str(hklsp))
G_aligned = R(phy, LabY) * tmp;
end
if (abs(phy) < conf.samry_max & abs(phx) < conf.samrx_max)
fprintf('Grain %d: Found %s reflection: diffry = %f samrx = %f samry = %f\n', gid(i), num2str(hklsp), theta, phx, phy)
fprintf('Positions in DCT scan: omega = %f (diffrz = %f), u = %f v = %f\n', omega, omega + conf.diffrz_offset, uv(1), uv(2))
% Calculate edge_on positions for slip-planes
tt_pars{tt_id}.nfdtx = conf.distance + conf.nfdtx_offset;
tt_pars{tt_id}.d3tz = conf.distance * tand(2 *theta) + conf.d3tz_offset;
tt_pars{tt_id}.diffry = -theta; % the id11 diffractometer can only rotate this way round...
tt_pars{tt_id}.samrx = phx;
tt_pars{tt_id}.samry = phy;
tt_pars{tt_id}.samtx = conf.samtx_offset - gr_cen_instr(1);
tt_pars{tt_id}.samty = conf.samty_offset - gr_cen_instr(2);
tt_pars{tt_id}.samtz = conf.samtz_offset - gr_cen_instr(3);
tt_pars{tt_id}.samrx_offset = conf.samrx_offset;
tt_pars{tt_id}.samry_offset = conf.samry_offset;
tt_pars{tt_id}.samtx_offset = conf.samtx_offset;
tt_pars{tt_id}.samty_offset = conf.samty_offset;
tt_pars{tt_id}.samtz_offset = conf.samtz_offset;
tt_pars{tt_id}.diffrz_offset= conf.diffrz_offset;
tt_pars{tt_id}.omegas = conf.omegas;
fprintf('\nDetector positions for %f mm distance: mv nfdtx %f, d3tz %f\n', conf.distance, tt_pars{tt_id}.nfdtx, tt_pars{tt_id}.d3tz)
fprintf('mv diffry %f samrx %f samry %f samtx %f samty %f samtz %f\n\n', -theta, phx, phy, tt_pars{tt_id}.samtx, ...
tt_pars{tt_id}.samty, tt_pars{tt_id}.samtz)
if ~isempty(conf.thetatype_slipplane)
for ii = 1 : size(slip_planes, 1)
RUT = R(phx, LabX);
RLT = R(phy, LabY);
slip_n = RZ0 * RLT * RUT * T' * slip_planes(ii, :)'; % tilted slip plane normal direction in Lab coordinates
% Here we try to solve an equation of type a * sin(om) + b * cos(om) = c
% <=> sqrt(a?? + b??) sin(om + beta) -> x = 2 atand(...)
% the following terms result from dot (beam_dir,
% pl_rot) = 0 ! (slip plane normal perpendicular to
% diffracted beam direction)
a = -cosd(theta) * slip_n(2);
b = cosd(theta) * slip_n(1);
c = -sind(theta) * slip_n(3);
D = a*a + b*b - c*c;
if D > 0
d = sqrt(D);
om(1) = 2 * atand((a + d) / (b + c));
om(2) = 2 * atand((a - d) / (b + c));
om = mod(om + 360, 360);
fprintf('The %s slipplane will be observed at omega = %f and %f (diffrz = %f and %f)\n', ...
num2str(g.allblobs.hklsp(slip_ind(ii), :)), om(1), om(2), om(1) + conf.diffrz_offset, om(2) + conf.diffrz_offset);
if conf.test_results
ROM1 = R(om(1), LabZ);
ROM2 = R(om(2), LabZ);
test1 = ROM1 * slip_n;
test2 = ROM2 * slip_n;
dotprod1 = test1' * [cosd(theta), 0 , sind(theta)]';
fprintf('[%f %f %f] * [cosd(theta), 0, sind(theta)] = %f\n', test1(1), test1(2), test1(3), dotprod1)
dotprod2 = test2' * [cosd(theta), 0 , sind(theta)]';
fprintf('[%f %f %f] * [cosd(theta), 0, sind(theta)] = %f\n', test2(1), test2(2), test2(3), dotprod2)
end
end
end
fprintf('\n\n')
end
tt_id = tt_id + 1;
end
end
end
out{i}.tt_pars = tt_pars;
out{i}.T = T;
end
end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment