Merge qspec funcs

horace_core/sqw/@Experiment/Experiment.m

@@ -304,6 +304,18 @@
             % compatibility with sqw interface
             obj.samples = val;
         end
+
+        function [qspec, en] = calc_qspec(obj)
+            % Compute the Q and E in spectrometer coordinates from
+            % experimental info.
+            efix = obj.get_efix();
+            emode = obj.get_emode();
+            en = obj.expdata(1).en;
+            det_direction = obj.detector_arrays.det_direction();
+            spec_to_rlu = obj.detector_arrays.dmat;
+
+            [qspec, en] = calc_qspec(det_direction, efix, en, emode);
+        end
     end
     %----------------------------------------------------------------------
     methods(Static)
@@ -590,10 +602,10 @@
             %
             %   - Multiple arguments can be passed, one for each run that
             %     constitutes the sqw object, by having one row per run
-            %   	i.e
-            %       	scalar      ---->   column vector (nrun elements)
+            %           i.e
+            %           scalar      ---->   column vector (nrun elements)
             %           row vector  ---->   2D array (nrun rows)
-            %        	string      ---->   cell array of strings
+            %           string      ---->   cell array of strings
             %
             % Throws if not valid form
             [args,npar] = check_and_expand_function_args_(varargin{:});
@@ -707,4 +719,3 @@
     end
     %----------------------------------------------------------------------
 end
-

horace_core/sqw/@rundatah/rundatah.m

@@ -164,7 +164,7 @@
             if size(obj.S,1)+1 == numel(en)
                 en = 0.5*(en(1:end-1)+en(2:end));
             end
-            [qspec,en]=calc_qspec_(detdcn,obj.efix,en,obj.emode);
+            [qspec,en]=calc_qspec(detdcn,obj.efix,en,obj.emode);
         end
 
         function [pix_or_data_range,pix,obj] = calc_projections(obj,detdcn)

horace_core/sqw/@sqw/calculate_qw_pixels2.m

@@ -27,130 +27,7 @@
           'Only a single sqw object is valid - cannot take an array of sqw objects')
 end
 
-c = neutron_constants;
-k_to_e = c.c_k_to_emev;
-
-% as column vectors
-idx = win.pix.all_indexes();
-irun = idx(1,:)';
-idet = idx(2,:)';
-ien = idx(3,:)';
-
-if ~iscell(win.header)
-    header = {win.header};
-else
-    header = win.header;
-end
-
-emode = cellfun(@(x) x.emode, header);
-
-if ~all(emode==emode(1))
-    error('HORACE:calculate_qw_pixels2:invalid_argument', ...
-          'Contributing runs to an sqw object must be all be direct geometry or all indirect geometry')
-end
-
-emode = emode(1);
-
-efix = cellfun(@(x) x.efix, header);
-eps_lo = cellfun(@(x) 0.5*(x.en(1)+x.en(2)), header);
-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);
-kfix = sqrt(efix/k_to_e);
-
-switch emode
-  case 1
-    ki = kfix(irun);
-    kf = sqrt((efix(irun)-eps_diff)/k_to_e);
-  case 2
-    ki = sqrt((efix(irun)+eps_diff)/k_to_e);
-    kf = kfix(irun);
-  otherwise
-    ki = kfix(irun);
-    kf = ki;
-end
-
-qw = cell(1, 4);
-qw(1:3) = calculate_q (ki, kf, detdcn(:, idet), spec_to_rlu(:, :, irun));
-qw{4} = eps_diff;
-
-end
-
-function [d_mat, 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)
-%
-% Input:
-% ------
-%   detpar      Detector parameter structure with fields as read by get_par
-%
-% Output:
-% -------
-%   d_mat       Matrix size [3, 3, ndet] to take coordinates in spectrometer
-%              frame and convert in detector frame.
-%
-%   detdcn      Direction of detector in spectrometer coordinates ([3 x ndet] array)
-%               [cos(phi); sin(phi).*cos(azim); sin(phi).sin(azim)]
-%
-% The detector frame is one with x axis along kf, y radially outwards. This is the
-% original Tobyfit detector frame.
-
-%% 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]);
-
-d_mat = [cp, cb.*sp, sb.*sp;...
-         -sp, cb.*cp, sb.*cp;...
-         zeros(1, 1, ndet), -sb, cb];
-
-detdcn = [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);
+[qspec, en] = win.experiment_info.calc_qspec();
+qw = {qspec(:, 1), qspec(:, 2), qspec(:, 3), en};
 
 end

horace_core/sqw/calc_qspec.m

@@ -0,0 +1,77 @@
+function [qspec, en]=calc_qspec(detdcn, efix, eps, emode)
+% Calculate the components of Q in reference frame fixed w.r.t. spectrometer
+%
+%   >> qspec = calc_qspec(detdcn, efix, eps, emode)
+%
+% Input:
+% ------
+%   detdcn  Direction of detector in spectrometer coordinates ([3 x ndet] array)
+%             [cos(phi); sin(phi).*cos(azim); sin(phi).sin(azim)]
+%   efix    incident energy for direct mode spectrometer or analyser(s)
+%           energy(s) for indirect.
+%   eps     bin centres for energy transfer energies
+%   emode   instrument and q-conversion type (1-direct, 2-indirect 0-
+%            elastic)
+%
+% Output:
+% -------
+%   qspec   Momentum in spectrometer coordinates
+%           (x-axis along ki, z-axis vertically upwards) ([3, ne*ndet] array)
+%   en      Energy transfer for all pixels ([1, ne*ndet] array)
+%
+%  Note: We sometimes use this routine with the energy bin boundaries replaced with
+%        bin centres i.e. have fudged the array data.en
+
+% T.G.Perring 15/6/07
+
+% Get components of Q in spectrometer frame (x || ki, z vertical)
+ne = numel(eps);
+ndet = size(detdcn, 2);
+
+c = neutron_constants;
+k_to_e = c.c_k_to_emev;  % used by calc_projections_c;
+
+switch emode
+  case 0 % Elastic
+    lambda=exp(eps);        % just pass the values as bin centres
+
+    k=(2*pi)./lambda;   % [ne x 1]
+    Q_by_k = repmat([1;0;0], [1, ndet]) - detdcn;   % [3 x ndet]
+    qspec = repmat(k', [3, ndet]).*reshape(repmat(reshape(Q_by_k, [3, 1, ndet]), [1, ne, 1]), [3, ne*ndet]);
+    en=zeros(1, ne*ndet);
+
+  case 1 % Direct Geometry
+    ki=sqrt(efix/k_to_e);
+    kf=sqrt((efix-eps)/k_to_e); % [ne x 1]
+    qspec = repmat([ki;0;0], [1, ne*ndet]) - ...
+        repmat(kf', [3, ndet]).*reshape(repmat(reshape(detdcn, [3, 1, ndet]), [1, ne, 1]), [3, ne*ndet]); %CM:detdcn shape wrong
+    en=repmat(eps(:)', 1, ndet);
+
+  case 2 % Indirect Geometry
+    %kf=sqrt(efix(:)'/k_to_e);
+    kf=sqrt(efix/k_to_e);
+    if numel(efix) > 1
+        ki = sqrt(bsxfun(@plus, efix(:)', eps(:))/k_to_e); % [ne x n_det]
+        if size(ki, 2) ~= ndet
+            error('HORACE:instr_proj:invalid_argument', ...
+                  'Number of detector''s energies in indirect mode(%d) must be equal to the number of detectors %d', ...
+                  size(ki, 2), ndet);
+        end
+        nde = ndet*ne;
+        qspec = [reshape(ki, 1, nde);zeros(2, nde)]  - ...
+            reshape(repmat(kf, [ne, 1, 3]), [nde, 3])'.*...
+            reshape(repmat(reshape(detdcn, [3, 1, ndet]), [1, ne, 1]), [3, nde]);
+    else
+        ki=sqrt((efix+eps(:))/k_to_e); % [ne x 1]
+        % building sequence Det1{dE1, dE2, ...dEN}, Det2{dE1, dE2, ...dEN}... DetN{dE1, dE2, ...dEN}
+        % note [3, _1_, ndet] vs [1, _ne_, 1]
+        qspec = repmat([ki';zeros(2, ne)], [1, ndet]) - ...
+            kf*reshape(repmat(reshape(detdcn, [3, 1, ndet]), [1, ne, 1]), [3, ne*ndet]);
+    end
+
+    en=repmat(eps(:)', 1, ndet);
+
+  otherwise
+    error('HORACE:instr_proj:invalid_argument', ...
+        'EMODE must =1 (direct geometry), =2 (indirect geometry), or =0 (elastic)')
+end