diff --git a/_test/test_sqw_class/test_migrated_apis.m b/_test/test_sqw_class/test_migrated_apis.m index a022921a1d..52f32dfb30 100644 --- a/_test/test_sqw_class/test_migrated_apis.m +++ b/_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'); diff --git a/horace_core/sqw/@sqw/calculate_qw_pixels2.m b/horace_core/sqw/@sqw/calculate_qw_pixels2.m index 3f0497171c..b0d8495549 100644 --- a/horace_core/sqw/@sqw/calculate_qw_pixels2.m +++ b/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 diff --git a/horace_core/sqw/@sqw/private/calculate_q.m b/horace_core/sqw/@sqw/private/calculate_q.m index 77ba68f52d..46989c371e 100644 --- a/horace_core/sqw/@sqw/private/calculate_q.m +++ b/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