diff --git a/zUtil_Topotomo/calc_pole_tilts_id11_index2.m b/zUtil_Topotomo/calc_pole_tilts_id11_index2.m
index 1bb3e1f2c4e684005a8435947dcc20919e3f0500..c7881f2efe7b33206e303edc3c3a4bf50158dd82 100644
--- a/zUtil_Topotomo/calc_pole_tilts_id11_index2.m
+++ b/zUtil_Topotomo/calc_pole_tilts_id11_index2.m
@@ -1,73 +1,166 @@
-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