Skip to content

Commit

Permalink
Re #1776 Merge branch '1776_kf_sphere' of https://github.com/pace-neu…
Browse files Browse the repository at this point in the history
…trons/Horace into 1776_kf_sphere
  • Loading branch information
abuts committed Jan 24, 2025
2 parents 5c8717d + 318ed2b commit b4666c4
Show file tree
Hide file tree
Showing 5 changed files with 237 additions and 16 deletions.
8 changes: 6 additions & 2 deletions herbert_core/utilities/scientific/neutron_constants.m
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
function s=neutron_constants
function s=neutron_constants(constant_name)
% Return a structure with constants for neutron units manipulation
%
% >> s=neutron_constants % return structure
Expand All @@ -24,7 +24,11 @@
end

if nargout>0
s=structure;
if narin == 0
s=structure;
else
s = structure.(constant_name);
end
else
nam=fieldnames(structure);
disp('Available constants:')
Expand Down
17 changes: 5 additions & 12 deletions horace_core/sqw/@SQWDnDBase/private/cut_parse_inputs_.m
Original file line number Diff line number Diff line change
Expand Up @@ -226,18 +226,11 @@
end

if opt.proj_given
% There are currently no situations where we want to define lattice in
% projection. so always take lattice from the source object
source_proj = obj.proj;
proj.do_check_combo_arg = false;

proj.alatt = source_proj.alatt;
proj.angdeg = source_proj.angdeg;

proj.do_check_combo_arg = true;
proj = proj.check_combo_arg();
%end

%
% Take source parameters, which projection may need to
% use in its transformation from source to target coordinate
% system.
proj = proj.copy_proj_param_from_source(obj);
else % it may be fewer parameters then actual dimensions and
% if no projection is given, we would like to append missing binning
% parameters with their default values.
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -178,7 +178,7 @@
obj = obj.check_combo_arg();
end
end
function obj = set_img_scales(obj,varargin)
function obj = set_img_scales(~,varargin)
error('HORACE:CurveProjBase:invalid_argument', ...
'You can not set image scales directly. Use projection type instead')
end
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -171,6 +171,24 @@
end
[obj,par] = init(obj,varargin{:});
end
function obj = copy_proj_param_from_source(obj,cut_source)
% a projection may need some information from the object,
% it used to process. For example, all projections used for cut
% need to know lattice parameters of sqw object used as source
%
% this is generic method, which may be overloaded by specific
% projections, which may want more information from the source
% sqw
source_proj = cut_source.proj;
%
%Retrived lattice from the source object
obj.do_check_combo_arg = false;
obj.alatt = source_proj.alatt;
obj.angdeg = source_proj.angdeg;

obj.do_check_combo_arg = true;
obj = obj.check_combo_arg();
end
%
function [obj,remains] = init(obj,varargin)
% Method normally used to initialize an empty object.
Expand Down Expand Up @@ -954,7 +972,7 @@
% for indexes to include
% en_inside -- 1D logical array, containing true, for
% contributing cells (n_cells = n_edges-1) for
% orthogonal 1D indexes on dE lattice
% orthogonal 1D indexes on dE lattice
%
% Uses knowledge about linear arrangement of 4-D array of indexes
% in memory and on disk
Expand Down
206 changes: 206 additions & 0 deletions horace_core/sqw/coord_transform/@kf_sphere_proj/kf_sphere_proj.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,206 @@
classdef kf_sphere_proj<sphere_proj
% Class defines spherical coordinate projection, used by cut_sqw
% to make spherical cuts.
%
% Unlike sphere_proj, which calculates spherical coordinates of scattrgint
% vector Q, kf_sphere_proj, calculates spherical coordinates of
% scattering vector kf, where Q = ki-kf;
%
% Usage (with positional parameters):
%
% >>sp = kf_sphere_proj(); %default construction
% >>sp = kf_sphere_proj(Ei);
% >>sp = kf_sphere_proj(Ei,u,v);
% >>sp = kf_sphere_proj(Ei,u,v,type);
% >>sp = kf_sphere_proj(Ei,u,v,type,alatt,angdeg);
% >>sp = kf_sphere_proj(Ei,u,v,type,alatt,angdeg,offset,label,title);
%
% Where:
% Ei -- Incident energy in direct spectrometer. Will work for indirect
% too though its meaning a bit unclear.
% u -- [1,3] vector of hkl direction of z-axis of the spherical
% coordinate system this projection defines. It defines the
% direction of ki vector (incident beam)
% The axis to calculate theta angle from.
% v -- [1,3] vector of hkl direction of x-axis of the spherical
% coordinate system, the axis to calculate Phi angle from.
%
% type-- 3-letter character array, defining the spherical
% coordinate system units (see type property below)
% alatt-- 3-vector of lattice parameters. Value will be ignored by cut.
% angdeg- 3-vector of lattice angles. Value will be ignored by cut.
% offset- 4-vector, defining hkldE value of centre of
% coordinates of the spherical coordinate
% system.
% label - 4-element cellarray, which defines axes labels
% title - character string to title the plots of cuts, obtained
% using this projection.
%
% all parameters may be provided as 'key',value pairs appearing in
% arbitrary order after positional parameters
% e.g.:
% >>sp = kf_sphere_proj(10,[1,0,0],[0,1,0],'arr','offset',[1,1,0]);
% >>sp = kf_sphere_proj(10,[1,0,0],'type','arr','v',[0,1,0],'offset',[1,1,0]);
%
% Default angular coordinates names and meaning of the coordinate system,
% defined by sphere_proj are chosen as follows:
% |kf| -- coordinate 1 is the modulus of the scattering momentum.
% theta -- coordinate 2, the angle between axis u
% and the direction of the kf.
% phi -- coordinate 3 is the angle between the projection of the
% scattering vector to the plane defined by vector v and
% perpendicular to u.
% dE -- coordinate 4 the energy transfer direction
%
% parent's class "type" property describes which scales are avaliable for
% each direction:
% for |Q|:
% 'a' -- Angstrom,
% 'r' -- scale = max(\vec{u}*\vec{e_h,e_k,e_l}) -- projection of u to
% unit vectors in hkl directions
% 'p' -- |u| = 1 -- i.e. scale = |u|
% 'h','k' or 'l' -- i.e. scale = (a*,b* or c*);
% for angular units theta, phi:
% 'd' - degree, 'r' -- radians
% For energy transfer:
% 'e'-energy transfer in meV (no other scaling so may be missing)
%
properties(Dependent)
Ei
end
properties(Access=protected)
Ei_ = [];
% cache for ki -- incident beam wave vector
ki_ = [1;0;0];
end
methods
function obj=kf_sphere_proj(varargin)
% Constrtuctor for spherical projection
% See init for the list of input parameters
%
obj = obj@sphere_proj();
obj.label = {'|kf|','\theta','\phi','En'};
if nargin>0
obj = obj.init(varargin{:});
end
end
%------------------------------------------------------------------
function ei = get.Ei(obj)
ei = obj.Ei_;
end
function obj = set.Ei(obj,val)
if ~isnumeric(val)||numel(val)
error('HORACE:kf_sphere_proj:invalid_arguments', ...
'Incident beam must be poisitive. Got %d',val);
end
if val<=0
error('HORACE:kf_sphere_proj:invalid_arguments', ...
'Incident beam must be poisitive. Got %d',val);
end

obj.Ei_ = val;
if obj.do_check_combo_arg_
obj = obj.check_combo_arg();
end
end
%------------------------------------------------------------------
% Particular implementation of aProjectionBase abstract interface
%------------------------------------------------------------------
function pix_transformed = transform_pix_to_img(obj,pix_data,varargin)
% Transform pixels expressed in crystal Cartesian coordinate systems
% into spherical coordinate system defined by the object
% properties
%
% Input:
% pix_data -- [3xNpix] or [4xNpix] array of pix coordinates
% expressed in crystal Cartesian coordinate system
% or instance of PixelDatBase class containing this
% information.
% Returns:
% pix_out -- [3xNpix or [4xNpix]Array the pixels coordinates
% transformed into spherical coordinate system
% defined by object properties
%
if isa(pix_data,'PixelDataBase')
pix_cc = pix_data.q_coordinates;
shift_ei = obj.offset(4) ~=0; % its ginored for the time being

ndim = 3;
input_is_obj = true;
else % if pix_input is 4-d, this will use 4-D matrix and shift
% if its 3-d -- matrix is 3-dimensional and energy is not shifted
% anyway
ndim = size(pix_input,1);
pix_cc = pix_input;
input_is_obj = false;
end
if ndim ==3
kf = obj.ki_-pix_cc;
else
kf = obj.ki_-pix_cc(1:3,:);
end
pix_transformed = transform_pix_to_img@sphere_proj(obj,kf,varargin{:});
if ndim > 3
if input_is_obj
pix_transformed = [pix_transformed;pix_data.dE];
else
pix_transformed = [pix_transformed;pix_data(4,:)];
end
end
end
function pix_cc = transform_img_to_pix(obj,pix_transformed,varargin)
% Transform pixels in image (spherical) coordinate system
% into crystal Cartesian system of pixels
kf = transform_img_to_pix@sphere_proj(obj,pix_transformed,varargin{:});
if size(kf,2) == 3
pix_cc = obj.ki_-kf;
else
pix_cc = obj.ki_-kf(1:3,:);
pix_cc = [pix_cc;kf(4,:)];
end
end
end

%=====================================================================
% SERIALIZABLE INTERFACE
%----------------------------------------------------------------------
methods
function flds = saveableFields(obj)
flds = saveableFields@sphere_proj(obj);
flds = ['Ei';flds(:)];
end
function ver = classVersion(~)
ver = 1;
end
function obj = check_combo_arg (obj)
% Check validity of interdependent fields
%
% >> obj = check_combo_arg(w)
%
% Throws HORACE:CurveProjBase:invalid_argument with the message
% suggesting the reason for failure if the inputs are incorrect
% w.r.t. each other.
%
% Normalizes input vectors to unity and constructs the
% transformation to new coordinate system when operation is
% successful
obj = check_combo_arg@sphere_proj(obj);
if obj.alatt_defined && obj.angdeg_defined
bm = obj.bmatrix();
else
bm = eye(3);
end
kf_mod=sqrt(obj.Ei_/neutron_constants('c_k_to_emev'));
obj.ki_ = bm*obj.u(:)*kf_mod;
end

end
methods(Static)
function obj = loadobj(S)
% boilerplate loadobj method, calling generic method of
% savable class. Useful for recovering class from a structure
obj = kf_sphere_proj();
obj = loadobj@serializable(S,obj);
end
end
end

0 comments on commit b4666c4

Please sign in to comment.