Fix issues with calculate_qw_pixels2 for recent horace

_test/test_sqw_class/test_migrated_apis.m

@@ -87,6 +87,34 @@ function test_calculate_qw_pixels(obj)
             assertEqualToTolWithSave(obj,qw,'tol',2.e-7);
         end
 
+        function test_calculate_qw_pixels2_unchanged(obj)
+            sqw_obj = read_sqw(obj.test_sqw_2d_fullpath);
+            qw = calculate_qw_pixels2(sqw_obj);
+
+            % Pixels unchanged so should be the same.
+            assertEqualToTol(qw, sqw_obj.pix.coordinates,'tol',1.e-3);
+        end
+
+        function test_calculate_qw_pixels2_symmetrised(obj)
+            sqw_obj = read_sqw(obj.test_sqw_2d_fullpath);
+
+            sym = SymopReflection([-1, 1, 0], [0, 0, 1], [-0.5, -0.5, 0]);
+
+            sqw_ref = sqw_obj.symmetrise_sqw(sym);
+            qw = calculate_qw_pixels2(sqw_ref);
+
+            % Resort because symmetrise reorders data.
+            [~, index] = sortrows(sqw_ref.pix.all_indexes');
+            [~, index_o] = sortrows(sqw_obj.pix.all_indexes');
+            recomp = qw(:, index);
+            orig = sqw_obj.pix.coordinates(:, index_o);
+            symm = sqw_ref.pix.coordinates(:, index);
+
+            % Despite reflection should regenerate original pixel locations
+            assertFalse(equal_to_tol(recomp, symm,'tol',1.e-3));
+            assertEqualToTol(recomp, orig,'tol',1.e-3);
+        end
+
         %% Cut
         function test_cut_sym(obj)
             skipTest('Incorrect test data for cut_sym');

horace_core/sqw/@sqw/calculate_qw_pixels2.m

@@ -36,8 +36,11 @@
 idet = idx(2,:)';
 ien = idx(3,:)';
 
+remap = containers.Map(unique(irun), 1:numel(win.header));
+irun = arrayfun(@(x) remap(x), irun);
+
 if ~iscell(win.header)
-    header = {win.header};
+    header = num2cell(win.header)';
 else
     header = win.header;
 end
@@ -56,15 +59,17 @@
 eps_hi = cellfun(@(x) 0.5*(x.en(end-1)+x.en(end)), header);
 n_en = cellfun(@(x) numel(x.en)-1, header);
 
+
 [~, ~, spec_to_rlu] = cellfun(...
     @(h) calc_proj_matrix(h.alatt, h.angdeg, ...
                           h.cu, h.cv, h.psi, ...
                           h.omega, h.dpsi, h.gl, h.gs), header, 'UniformOutput', false);
+
 % Join in 3rd rank leading to n x n x nhead
 spec_to_rlu = cat(3, spec_to_rlu{:});
 
 eps_diff = (eps_lo(irun) .* (n_en(irun) - ien) + eps_hi(irun) .* (ien - 1)) ./ (n_en(irun) - 1);
-[~, detdcn] = spec_coords_to_det(win.detpar);
+detdcn = spec_coords_to_det(win.detpar);
 kfix = sqrt(efix/k_to_e);
 
 switch emode
@@ -79,13 +84,15 @@
     kf = ki;
 end
 
-qw = cell(1, 4);
-qw(1:3) = calculate_q (ki, kf, detdcn(:, idet), spec_to_rlu(:, :, irun));
+qw = cell(4, 1);
+qw(1:3) = calculate_q(ki, kf, detdcn(:, idet), spec_to_rlu(:, :, irun));
 qw{4} = eps_diff;
+qw = cat(2, qw{:})';
+qw = win.data.proj.transform_hkl_to_pix(qw);
 
 end
 
-function [d_mat, detdcn] = spec_coords_to_det (detpar)
+function detdcn = spec_coords_to_det (detpar)
 % Matrix to convert coordinates in spectrometer (or laboratory) frame into detector frame
 %
 %   >> d_mat = spec_coords_to_det (detpar)
@@ -107,50 +114,17 @@
 
 %% TODO: Investigate use of transform_pix_to_hkl
 
-ndet = numel(detpar.x2);
-
-cp = reshape(cosd(detpar.phi), [1, 1, ndet]);
-sp = reshape(sind(detpar.phi), [1, 1, ndet]);
-cb = reshape(cosd(detpar.azim), [1, 1, ndet]);
-sb = reshape(sind(detpar.azim), [1, 1, ndet]);
+detdcn = {};
+for i=1:detpar.n_objects
+    ndet = numel(detpar{i}.x2);
 
-d_mat = [cp, cb.*sp, sb.*sp;...
-         -sp, cb.*cp, sb.*cp;...
-         zeros(1, 1, ndet), -sb, cb];
-
-detdcn = [cp; cb.*sp; sb.*sp];
+    cp = cosd(detpar{i}.phi);
+    sp = sind(detpar{i}.phi);
+    cb = cosd(detpar{i}.azim);
+    sb = sind(detpar{i}.azim);
 
+    detdcn{i} = [cp, cb.*sp, sb.*sp];
 end
 
-function q = calculate_q (ki, kf, detdcn, spec_to_rlu)
-% Calculate qh, qk, ql for direct geometry instrument
-%
-%   >> q = calculate_q (ki, kf, detdcn, spec_to_rlu)
-%
-% Input:
-% ------
-%   ki          Incident wavevectors for each point [Column vector]
-%   kf          Final wavevectors for each point    [Column vector]
-%   detdcn      Array of unit vectors in the direction of the detectors
-%               Size is [3, npnt]
-%   spec_to_rlu Matrix to convert momentum in spectrometer coordinates to
-%               components in r.l.u.:
-%                   v_rlu = spec_to_rlu * v_spec
-%               Size is [3, 3, npnt]
-%
-% Output:
-% -------
-%   q           Components of momentum (in rlu) for each point
-%               [Cell array of column vectors]
-%               i.e. q{1}=qh, q{2}=qk, q{3}=ql
-
-% Use in-place working to save memory (note: bsxfun not needed from 2016b an onwards)
-qtmp = -kf' .* detdcn;
-qtmp(1, :) = ki' + qtmp(1, :);            % qspec proper now
-qtmp = mtimesx_horace (spec_to_rlu, reshape(qtmp, [3, 1, numel(ki)]));
-qtmp = squeeze(qtmp);
-
-% Package output
-q = num2cell(qtmp', 1);
-
+detdcn = cat(1, detdcn{:})';
 end

horace_core/sqw/@sqw/private/calculate_q.m

@@ -21,13 +21,12 @@
 %               i.e. q{1}=qh, q{2}=qk, q{3}=ql
 
 % Use in-place working to save memory (note: bsxfun not needed from 2016b an onwards)
-qtmp = bsxfun(@times,-kf',detdcn);      % -kf in spectrometer axes
-qtmp(1,:) = ki' + qtmp(1,:);            % qspec proper now
-qtmp = mtimesx_horace (spec_to_rlu,reshape(qtmp,[3,1,numel(ki)]));
-qtmp = squeeze(qtmp);
+q = -kf' .* detdcn;
+q(1,:) = ki' + q(1,:);            % qspec proper now
+q = mtimesx_horace (spec_to_rlu, reshape(q, [3, 1, numel(ki)]));
+q = squeeze(q);
 
 % Package output
-q = cell(1,3);
-q{1} = qtmp(1,:)';
-q{2} = qtmp(2,:)';
-q{3} = qtmp(3,:)';
+q = num2cell(q', 1);
+
+end