diff --git a/+sw_tests/+unit_tests/unittest_spinw_spec2MDHisto.m b/+sw_tests/+unit_tests/unittest_spinw_spec2MDHisto.m new file mode 100644 index 000000000..1511d2730 --- /dev/null +++ b/+sw_tests/+unit_tests/unittest_spinw_spec2MDHisto.m @@ -0,0 +1,55 @@ +classdef unittest_spinw_spec2MDHisto < sw_tests.unit_tests.unittest_super + properties + swModel = []; + tmpdir = ''; + testfilename = ''; + nsteps = {100}; + end + + properties (TestParameter) + testpars = struct(... + 'test_1_0_0', struct('q0', [0 0 0], 'qdir', [1 0 0], 'proj', [[1 0 0]' [0 1 0]' [0 0 1]'], 'nxs', 'test100mdh.nxs'), ... + 'test_1_1_0', struct('q0', [0 0 0], 'qdir', [1 1 0], 'proj', [[1 1 0]' [1 -1 0]' [0 0 1]'], 'nxs', 'test110mdh.nxs'), ... + 'test_1_1_1', struct('q0', [0 0 0], 'qdir', [1 1 1], 'proj', [[1 1 1]' [1 -1 0]' [1 1 -2]'], 'nxs', 'test111mdh.nxs'), ... + 'test_1_1_2', struct('q0', [0 0 2], 'qdir', [1 1 0], 'proj', [[1 1 0]' [1 -1 0]' [0 0 1]'], 'nxs', 'test112mdh.nxs'), ... + 'test_1_1_2_2', struct('q0', [0 0 2], 'qdir', [1 1 0], 'proj', [[1 -1 0]' [1 1 0]' [0 0 1]'], 'nxs', 'test112_2mdh.nxs'), ... + 'test_2_2_2', struct('q0', [2 2 2], 'qdir', [1 1 0], 'proj', [[1 1 0]' [1 -1 0]' [0 0 1]'], 'nxs', 'test222mdh.nxs')); + end + + methods (TestClassSetup) + function setup_model(testCase) + testCase.swModel = sw_model('triAF', 1); + end + function setup_tempdir(testCase) + testCase.tmpdir = tempdir; + end + end + + methods (TestMethodTeardown) + function remove_tmpdir(testCase) + delete(testCase.testfilename); + end + end + + methods (Test) + function test_qdirs(testCase, testpars) + q0 = testpars.q0; + qdir = testpars.qdir; + proj = testpars.proj; + testCase.disable_warnings('spinw:spinwave:NonPosDefHamiltonian'); + spec = sw_egrid(spinwave(testCase.swModel, {q0 q0+qdir testCase.nsteps{1}})); + dproj = [norm((qdir-q0))/testCase.nsteps{1}, 1e-6, 1e-6]; + testCase.testfilename = fullfile(testCase.tmpdir, testpars.nxs); + sw_spec2MDHisto(spec, proj, dproj, testCase.testfilename); + end + + function test_non_ortho(testCase) + q0 = [0 0 2]; + qdir = [1 1 0]; + spec = sw_egrid(spinwave(testCase.swModel, {q0 q0+qdir testCase.nsteps{1}})); + proj = [qdir(:) [1 0 0]' [0 0 1]']; + dproj = [1e-6, norm((qdir-q0))/testCase.nsteps{1}, 1e-6]; + verifyError(testCase,@() sw_spec2MDHisto(spec, proj, dproj, 'tmp/test_blank.nxs'), "read_struct:nonorthogonal") + end + end +end \ No newline at end of file diff --git a/swfiles/sw_spec2MDHisto.m b/swfiles/sw_spec2MDHisto.m new file mode 100644 index 000000000..8c8eea9f2 --- /dev/null +++ b/swfiles/sw_spec2MDHisto.m @@ -0,0 +1,286 @@ +function sw_spec2MDHisto(spectra,proj, dproj, filename) +% saves spectrum to MDHisto +% +% ### Syntax +% +% sw_spec2MDHisto(spectra,proj,dproj,filename)` +% +% ### Description +% +% `sw_spec2MDHisto(spectra,proj,dproj,filename)` saves a +% spectrum that is calculated by sw_egrid +% +% ### Input Arguments +% spectra: a structure calculated by sw_egrid +% +% proj: a 3x3 matrix defining an orthogonal coordinate system +% where each column is a vector defining the orientation +% of the view. One of the vectors must be identical to the Q axis +% defined by the direction of the calculation. +% +% dproj: is a 3 vector that is the bin size in each of the +% directions defined in proj. For the direction of the +% calculation, the value used is internally calcualted from the spectrum. +% It is wise to enter the step size for clarity. +% +% filename: is the name of the nexus file. It will overwrite the existing +% file if one already exists +% +% Example: +% q0 = [0 0 0]; +% qdir = [1 0 0]; +% nsteps = 100; +% spec = sw_egrid(spinwave(sw_model('triAF', 1), {q0 q0+qdir nsteps})) +% proj = [qdir(:) [0 1 0]' [0 0 1]']; +% dproj = [(qdir(1)-q0(1)/steps, 1e-6, 1e-6]; +% sw_spec2MDHisto(spec, proj, dproj, 'testmdh.nxs'); +% Note that: +% (1) In the call to `spinwave`, only one q-direction may be specified +% e.g. the HKL specifier must be of the form {q0 q0+qdir nsteps} +% (2) one column in the `proj` matrix must be the q-direction used in +% `spinwave` (e.g. `qdir`). + + +if nargin==0 + swhelp sw_spec2MDHisto + return +end +[unit_cell,Bmat,proj_out,D,dat,proj,name] = read_struct(spectra,proj,dproj); +%check if hdf file exists and delete if it does. +if exist(filename,'file') + delete(filename) +end + +h5createnwrite(filename,'/MDHistoWorkspace/coordinate_system',3); % None = 0, QLab = 1, QSample = 2, HKL = 3 +h5createnwrite(filename,'/MDHistoWorkspace/visual_normalization',0); +h5writeatt(filename,'/MDHistoWorkspace','NX_class','NXentry'); +h5writeatt(filename,'/MDHistoWorkspace','Qconvention','Inelastic'); +h5writeatt(filename,'/MDHistoWorkspace','SaveMDVersion', 2); +% write data +rtpth = NXScreategroup(filename,'/MDHistoWorkspace','data','NXdata'); +% write D dimensions +Dszs=zeros(1,4); +Dstrcell={}; +for idx=1:length(D) + szd=size(D{idx}); + Dszs(idx)=szd(2)-1; + Dstrcell{idx} = strcat('D',num2str(idx-1)); + Dpth = strcat(rtpth,'/',Dstrcell{idx}); + h5createnwritevec(filename,rtpth,Dstrcell{idx},D{idx}); + if idx '); +h5createnwritevec(filename,smplpth,'num_oriented_lattice',int32(1)) +h5createnwritevec(filename,smplpth,'num_other_samples',int32(0)) +h5createnwritevec(filename,smplpth,'geom_height', 0) +h5createnwritevec(filename,smplpth,'geom_id', int32(0)) +h5createnwritevec(filename,smplpth,'geom_thickness', 0) +h5createnwritevec(filename,smplpth,'geom_width' ,0) +%write material +mtl_pth = NXScreategroup(filename,smplpth,'material','NXdata'); +h5writeatt(filename,mtl_pth,'formulaStyle','empty') +h5writeatt(filename,mtl_pth,'name',' ') +h5writeatt(filename,mtl_pth,'version',int32(2)) +h5createnwritevec(filename,mtl_pth,'packing_fraction',1) +h5createnwritevec(filename,mtl_pth,'number_density',0) +h5createnwritevec(filename,mtl_pth,'pressure',0) +h5createnwritevec(filename,mtl_pth,'temperature',0) +%write crystal lattice +OL_pth = NXScreategroup(filename,smplpth,'oriented_lattice','NXcrystal'); +%write lattice parameters +u_parm_names={'a','b','c','alpha','beta','gamma'}; +for idx=1:length(u_parm_names) + parm_name = strcat('unit_cell_',u_parm_names{idx}); + h5createnwritevec(filename,OL_pth,parm_name,unit_cell(idx)) +end +%write orientation matrix +om_path =strcat(OL_pth,'/orientation_matrix'); +h5createnwrite(filename,om_path,Bmat); +%write instrument +instr_pth = NXScreategroup(filename,exppth,'instrument','NXinstrument'); +h5writeatt(filename,instr_pth,'version',int32(1)) +h5createnwritevec(filename,instr_pth,'name','SEQUOIA'); +end + +function h5createnwrite(filename,path,val) +dtyp = class(val); +h5create(filename,path,size(val),'Datatype',dtyp); +h5write(filename,path,val); +end + +function h5createnwritevec(filename,group,ds,val) +fh = H5F.open(filename,'H5F_ACC_RDWR','H5P_DEFAULT'); +sph = H5S.create_simple(1,length(val),length(val)); +gsh = H5G.open(fh,group); +memtype = 'H5ML_DEFAULT'; +switch class(val) + case 'int32' + H5type = 'H5T_NATIVE_INT'; + + case 'int8' + H5type = 'H5T_NATIVE_CHAR'; + case 'char' + dims= size(val); + sph = H5S.create_simple(1,fliplr(dims(1)),[]); + H5type = H5T.copy('H5T_FORTRAN_S1'); + H5T.set_size(H5type,(dims(2)+1)) + memtype = H5T.copy ('H5T_C_S1'); + H5T.set_size(memtype,dims(2)); + otherwise + H5type = 'H5T_NATIVE_DOUBLE'; +end +dsh = H5D.create(gsh,ds,H5type,sph,'H5P_DEFAULT'); +H5D.write(dsh,memtype,'H5S_ALL','H5S_ALL','H5P_DEFAULT',val) +H5S.close(sph) +H5D.close(dsh) +H5G.close(gsh) +H5F.close(fh) +end + +function pthout = NXScreategroup(filename,pth,group,NX_class) +% ### Syntax + +% pthout = NXScreategroup(filename,pth,group,NX_class) + +% ### Description +% +% create a group with attributes of a nexus class +% +% ### Input Arguments +% filename, name of hdf5 file +% pth, path to where the group should be created +% group the group name +% NX_class a string containing a valid NX_class definition +% +% ### Output Arguments +% returns a path to the group + +fh = H5F.open(filename,'H5F_ACC_RDWR','H5P_DEFAULT'); +pthh = H5G.open(fh,pth); +gh = H5G.create(pthh,group,'H5P_DEFAULT','H5P_DEFAULT','H5P_DEFAULT'); +pthout =strcat([pth,'/',group]); +H5G.close(gh) +H5G.close(pthh) +H5F.close(fh) +h5writeatt(filename,pthout,'NX_class',NX_class) +end + +function writeNXlog(filename,log_pth,log_nm,value,units) +% ### Description +% +% create and write a Nexus log +% +% ### Input Arguments +%filename, name of hdf5 file +%log_pth, path to where the log should be created +%log_nm the name of the log +% value the vlaue of the log +% the units if any +pth = NXScreategroup(filename,log_pth,log_nm,'NXlog' ); +h5createnwritevec(filename,pth,'value',value) +h5writeatt(filename,pth,'units',units) + +end +function [latt_parms,Bmat,proj_out,D,signal,proj,name] = read_struct(dstruct,proj,dproj) + +% ### Input Arguments +% dstruct: is the spinw structure +% proj: is the viewing projection matrix each column is a different axis. +% One axis must be parallel to the drection of propogation in the +% dstruct +% drproj: is the size of the bin in each direction (ignored for the direction +% of the cut). +% +% ### Output Arguments +% latt_parms: the lattice parameters from the spinw structure +% Bmat: a matrix for converting from inverse angstroms to rlu +% proj_out: +% D: a cell array of the number of steps in each direction +% signal: the signal array from the spinw spec strcuture +% name : the chemical formula from the spinW file + check_ortho(proj) + fnames=fieldnames(dstruct); + objnum = find(strcmp(fnames,'obj')); + swobj = dstruct.(subsref(fnames,substruct('{}',{objnum}))); + latt_parms = abc(swobj); + M = basisvector(swobj); + Bmat = inv(M); + name =formula(swobj).chemform; + %proj_out = proj(:); + hkls = dstruct.hkl; + hkls_sz = size(hkls); + dir_vec = hkls(:,hkls_sz(2))-hkls(:,1); + %qout = hkls'/dir_vec'; + D={}; + for idx=1:3 + procjv = proj(:,idx)/norm(proj(:,idx)); + + if abs(norm(cross(dir_vec,procjv)))> 1e-6 + dtmp = dot(hkls(:,2),procjv); + D{idx} = dtmp+dproj(idx)/2.*[-1 1]; + else + hkl_proj = hkls'/dir_vec'; + dhkl = hkl_proj(2)-hkl_proj(1); % get the spacing along the q axis + %hkl_proj = dot(hkls(:,1),procjv); + D{idx} = zeros([1,length(hkl_proj)+1]); + D{idx}(1:length(hkl_proj)) = hkl_proj-dhkl/2; + D{idx}(length(D{idx})) = hkl_proj(length(hkl_proj))+dhkl/2;% change to bin boundaries + proj(:,idx) = dir_vec; %set varying projection vector to spectra object + end + end + D{4}=dstruct.Evect; + + signal = dstruct.swConv; + proj_out = proj(:); +end + +function check_ortho(mat) +% check if the three column vectors in proj are orthogonal to each other +sm = size(mat); + for idx = 1:sm(2) + for idx2 = 1:sm(2) + if idx~=idx2 + if norm(dot(mat(:,idx),mat(:,idx2))) >1e-6 + error("read_struct:nonorthogonal","the 3 vectors in proj must form an orthogonal basis set") + end + end + end + end +end \ No newline at end of file