diff --git a/+sw_tests/+system_tests/systemtest_spinwave_incommensurate_and_supercell_consistency.m b/+sw_tests/+system_tests/systemtest_spinwave_incommensurate_and_supercell_consistency.m index 2ce7ec5b6..7eed5dce7 100644 --- a/+sw_tests/+system_tests/systemtest_spinwave_incommensurate_and_supercell_consistency.m +++ b/+sw_tests/+system_tests/systemtest_spinwave_incommensurate_and_supercell_consistency.m @@ -4,14 +4,19 @@ reference_data_file = []; tol = 1e-5; end - methods - function om = remove_ghosts(testCase, spec) - om = spec.omega(find(abs(spec.Sperp) > testCase.tol)); - om = sort(unique(round(om / testCase.tol) * testCase.tol)); + methods (Static) + function om = remove_ghosts(spec, tol) + om = spec.omega(find(abs(spec.Sperp) > tol)); + om = sort(unique(round(om / tol) * tol)); end - + end + methods function assert_super_and_incom_consistency(testCase, swobj, ... - spec_super, spec_incom) + spec_super, spec_incom, ... + ghost_tol) + if nargin < 5 + ghost_tol = testCase.tol; + end % test cross-section in q,En bins testCase.verify_test_data(spec_incom.swConv, ... spec_super.swConv) @@ -22,8 +27,8 @@ function assert_super_and_incom_consistency(testCase, swobj, ... 6*n_matom); testCase.assertEqual(size(spec_super.Sperp, 1), ... 2*prod(nExt)*n_matom); - testCase.assertEqual(testCase.remove_ghosts(spec_super),... - testCase.remove_ghosts(spec_incom)); + testCase.assertEqual(testCase.remove_ghosts(spec_super, ghost_tol), ... + testCase.remove_ghosts(spec_incom, ghost_tol)); end end @@ -35,31 +40,33 @@ function test_AFM_kagome(testCase) 'angled',[90 90 120], 'sym','P -3') AF33kagome.addatom('r',[1/2 0 0],'S', 1, 'label','MCu1') AF33kagome.gencoupling('maxDistance',7); - AF33kagome.addmatrix('label','J1','value',1.00) + AF33kagome.addmatrix('label','J1','value',1) AF33kagome.addcoupling('mat','J1','bond',1); % sqrt3 x sqrt(3) magnetic structure k = [-1/3 -1/3 0]; n = [0, 0, 1]; S = [0 0 -1; 1 1 -1; 0 0 0]; % binning for spinwave spectrum - qarg = {[-1/2 0 0] [0 0 0] [1/2 1/2 0] 3}; - evec = 0:0.5:2; + qarg = {[-1/2 0 0] [0 0 0] [1/2 1/2 0] 50}; + evec = 0:0.1:1.5; % use structural unit cell with incommensurate k AF33kagome.genmagstr('mode','helical','unit','lu', 'k', k,... 'n',n, 'S', S, 'nExt',[1 1 1]); testCase.disable_warnings('spinw:spinwave:NonPosDefHamiltonian'); spec_incom = AF33kagome.spinwave(qarg, 'hermit', true); - spec_incom = sw_egrid(spec_incom, 'component','Sperp', 'Evect',evec); + spec_incom = sw_egrid(spec_incom, 'component','Sperp', 'Evect', evec, ... + 'zeroEnergyTol', 1e-2); % use supercell k=0 structure AF33kagome.genmagstr('mode','helical','unit','lu', 'k', k,... 'n',n, 'S', S, 'nExt', [3,3,1]); + spec_super = AF33kagome.spinwave(qarg, 'hermit', true); - spec_super = sw_egrid(spec_super, 'component','Sperp', 'Evect',evec); + spec_super = sw_egrid(spec_super, 'component','Sperp', 'Evect', evec); testCase.assert_super_and_incom_consistency(AF33kagome, ... spec_super, ... - spec_incom); + spec_incom, 5e-2); end function test_two_matom_per_unit_cell(testCase) @@ -81,8 +88,8 @@ function test_two_matom_per_unit_cell(testCase) k = [1/2, 0, 0]; S = [0 0;1 1;0 0]; % binning for spinwave spectrum - qarg = {[0 0 0] [1/2 0 0] 5}; - evec = 0:1.5:5; + qarg = {[0 0 0] [0, 0.5, 0] 5}; + evec = 0:0.5:5; % use structural unit cell with incommensurate k testCase.disable_warnings('spinw:spinwave:NonPosDefHamiltonian', ... diff --git a/+sw_tests/+unit_tests/unittest_sw_egrid.m b/+sw_tests/+unit_tests/unittest_sw_egrid.m new file mode 100644 index 000000000..e45266c13 --- /dev/null +++ b/+sw_tests/+unit_tests/unittest_sw_egrid.m @@ -0,0 +1,373 @@ +classdef unittest_sw_egrid < sw_tests.unit_tests.unittest_super + % Runs through unit test for @spinw/spinwave.m + + properties + swobj = []; + swobj_tri = []; + spectrum = struct(); + sw_egrid_out = struct(); + sw_egrid_out_pol = struct(); + sw_egrid_out_sperp = struct(); + qh5 = [0:0.25:1; zeros(2,5)]; + end + + properties (TestParameter) + % Components that require sw_neutron be called first + pol_components = {'Mxx', 'Pxy', 'Pz'}; + end + + methods (TestClassSetup) + function setup_chain_model(testCase) + % Just create a very simple FM 1D chain model + testCase.swobj = spinw; + testCase.swobj.genlattice('lat_const', [3 8 8], 'angled', [90 90 90]); + testCase.swobj.addatom('r', [0 0 0],'S', 1, 'label', 'MNi2'); + testCase.swobj.gencoupling('maxDistance', 7); + testCase.swobj.addmatrix('value', -eye(3), 'label', 'Ja'); + testCase.swobj.addcoupling('mat', 'Ja', 'bond', 1); + testCase.swobj.genmagstr('mode', 'direct', 'k', [0 0 0], ... + 'S', [0; 1; 0]); + end + end + + methods (TestMethodSetup) + function setup_default_spectrum_and_egrid(testCase) + % Spectrum input to sw_egrid + testCase.spectrum.obj = testCase.swobj; + testCase.spectrum.formfact = false; + testCase.spectrum.incomm = false; + testCase.spectrum.helical = false; + testCase.spectrum.norm = false; + testCase.spectrum.nformula = int32(0); + testCase.spectrum.param = struct('notwin', true, ... + 'sortMode', true, ... + 'tol', 1e-4, ... + 'omega_tol', 1e-5, ... + 'hermit', true); + testCase.spectrum.title = 'Numerical LSWT spectrum'; + testCase.spectrum.gtensor = false; + testCase.spectrum.datestart = '01-Jan-2023 00:00:01'; + testCase.spectrum.dateend = '01-Jan-2023 00:00:02'; + testCase.spectrum.hkl = [0:0.25:1; zeros(2,5)]; + testCase.spectrum.hklA = testCase.spectrum.hkl*2/3*pi; + testCase.spectrum.omega = [ 1e-5 2. 4. 2. -1e-5; ... + -1e-5 -2. -4. -2. 1e-5]; + testCase.spectrum.Sab = zeros(3, 3, 2, 5); + Sab1 = [0.5 0 0.5j; 0 0 0; -0.5j 0 0.5]; + Sab2 = [0.5 0 -0.5j; 0 0 0; 0.5j 0 0.5]; + testCase.spectrum.Sab(:, :, 1, 1:4) = repmat(Sab1, 1, 1, 1, 4); + testCase.spectrum.Sab(:, :, 2, 5) = Sab1; + testCase.spectrum.Sab(:, :, 2, 1:4) = repmat(Sab2, 1, 1, 1, 4); + testCase.spectrum.Sab(:, :, 1, 5) = Sab2; + + % Output from sw_egrid when 'component' is specified + testCase.sw_egrid_out = testCase.spectrum; + testCase.sw_egrid_out.param.sumtwin = true; + testCase.sw_egrid_out.T = 0; + testCase.sw_egrid_out.Evect = linspace(0, 4.4, 501); + testCase.sw_egrid_out.swConv = zeros(500, 5); + testCase.sw_egrid_out.swConv([728, 1455, 1728]) = 0.5; + testCase.sw_egrid_out.swInt = 0.5*ones(2, 5); + testCase.sw_egrid_out.component = 'Sperp'; + + % Default output from sw_egrid - when Sperp is used there are + % extra fields + testCase.sw_egrid_out_sperp = testCase.sw_egrid_out; + testCase.sw_egrid_out_sperp.intP = []; + testCase.sw_egrid_out_sperp.Mab = []; + testCase.sw_egrid_out_sperp.Pab = []; + testCase.sw_egrid_out_sperp.Sperp = 0.5*ones(2, 5); + testCase.sw_egrid_out_sperp.param.n = [0 0 1]; + testCase.sw_egrid_out_sperp.param.pol = false; + testCase.sw_egrid_out_sperp.param.uv = {}; + testCase.sw_egrid_out_sperp.param.sumtwin = true; + + % Output from sw_egrid when polarisation required - when + % sw_neutron has to be called first + testCase.sw_egrid_out_pol = testCase.sw_egrid_out_sperp; + testCase.sw_egrid_out_pol.param.pol = true; + testCase.sw_egrid_out_pol.intP = 0.5*ones(3, 2, 5); + testCase.sw_egrid_out_pol.Pab = repmat(diag([-0.5 -0.5 0.5]), 1, 1, 2, 5); + Mab_mat = [0.5 0 1i*0.5; 0 0 0; -1i*0.5 0 0.5]; + testCase.sw_egrid_out_pol.Mab = repmat(Mab_mat, 1, 1, 2, 5); + testCase.sw_egrid_out_pol.Mab(:, :, 1, 5) = conj(Mab_mat); + testCase.sw_egrid_out_pol.Mab(:, :, 2, :) = conj(testCase.sw_egrid_out_pol.Mab(:, :, 1, :)); + + end + end + + methods (Test) + function test_noInput(testCase) + % Tests that if sw_egrid is called with no input, it calls the help + % First mock the help call + help_function = sw_tests.utilities.mock_function('swhelp'); + sw_egrid(); + testCase.assertEqual(help_function.n_calls, 1); + testCase.assertEqual(help_function.arguments, {{'sw_egrid'}}); + end + function test_pol_component_no_sw_neutron_causes_error(testCase, pol_components) + testCase.verifyError(... + @() sw_egrid(testCase.spectrum, 'component', pol_components), ... + 'sw_egrid:WrongInput'); + end + function test_invalid_component(testCase) + testCase.verifyError(... + @() sw_egrid(testCase.spectrum, 'component', 'Sab'), ... + 'sw_parstr:WrongString'); + end + function test_epsilon_deprecated_warning(testCase) + testCase.verifyWarning(... + @() sw_egrid(testCase.spectrum, 'epsilon', 1e-5), ... + 'sw_egrid:DeprecationWarning'); + end + function test_defaults(testCase) + out = sw_egrid(testCase.spectrum); + testCase.verify_obj(out, testCase.sw_egrid_out_sperp); + end + function test_Sperp(testCase) + out = sw_egrid(testCase.spectrum, 'component', 'Sperp'); + testCase.verify_obj(out, testCase.sw_egrid_out_sperp); + end + function test_Sxx(testCase) + component = 'Sxx'; + expected_out = testCase.sw_egrid_out; + expected_out.component = component; + out = sw_egrid(testCase.spectrum, 'component', component); + testCase.verify_obj(out, expected_out); + end + function test_Sxy(testCase) + component = 'Sxy'; + expected_out = testCase.sw_egrid_out; + expected_out.component = component; + expected_out.swInt = zeros(2, 5); + expected_out.swConv = zeros(500, 5); + out = sw_egrid(testCase.spectrum, 'component', component); + testCase.verify_obj(out, expected_out); + end + function test_diff_Sxx_Szz(testCase) + component = 'Sxx-Szz'; + expected_out = testCase.sw_egrid_out; + expected_out.component = component; + expected_out.swInt = zeros(2, 5); + expected_out.swConv = zeros(500, 5); + out = sw_egrid(testCase.spectrum, 'component', component); + testCase.verify_obj(out, expected_out); + end + function test_sum_Sxx_Szz_Sxy(testCase) + component = 'Sxx+Szz+Sxy'; + expected_out = testCase.sw_egrid_out; + expected_out.component = component; + expected_out.swInt = 2*expected_out.swInt; + expected_out.swConv = 2*expected_out.swConv; + out = sw_egrid(testCase.spectrum, 'component', component); + testCase.verify_obj(out, expected_out); + end + function test_cell_array_component(testCase) + component = {'Sxx', 'Sxy'}; + expected_out = testCase.sw_egrid_out; + expected_out.component = component; + expected_out.swInt = {expected_out.swInt; zeros(2, 5)}; + expected_out.swConv = {expected_out.swConv; zeros(500, 5)}; + out = sw_egrid(testCase.spectrum, 'component', component); + testCase.verify_obj(out, expected_out); + end + function test_Mzz(testCase) + component = 'Mzz'; + expected_out = testCase.sw_egrid_out_pol; + expected_out.component = component; + + neutron_out = sw_neutron(testCase.spectrum, 'pol', true); + out = sw_egrid(neutron_out, 'component', component); + testCase.verify_obj(out, expected_out); + end + function test_sum_Pzz_Mxx(testCase) + component = 'Pzz+Mxx'; + expected_out = testCase.sw_egrid_out_pol; + expected_out.component = component; + expected_out.swInt = 2*expected_out.swInt; + expected_out.swConv = 2*expected_out.swConv; + + neutron_out = sw_neutron(testCase.spectrum, 'pol', true); + out = sw_egrid(neutron_out, 'component', component); + testCase.verify_obj(out, expected_out); + end + + function test_Px(testCase) + component = 'Px'; + expected_out = testCase.sw_egrid_out_pol; + expected_out.component = component; + neutron_out = sw_neutron(testCase.spectrum, 'pol', true); + out = sw_egrid(neutron_out, 'component', component); + testCase.verify_obj(out, expected_out); + end + + function test_fName_component(testCase) + % add field to spectrum + component = 'Sperp'; + spectrum = testCase.spectrum; + spectrum.(component) = ones(2,5); % double actual + + % note do not add fields normally added to output when Sperp + expected_out = testCase.sw_egrid_out; + expected_out.component = component; + expected_out.(component) = spectrum.(component); + expected_out.swInt = 2*expected_out.swInt; + expected_out.swConv = 2*expected_out.swConv; + + out = sw_egrid(spectrum, 'component', component); + + testCase.verify_obj(out, expected_out); + end + function test_Evect(testCase) + Evect = linspace(1, 3, 201); + expected_out = testCase.sw_egrid_out_sperp; + expected_out.swConv = zeros(200, 5); + expected_out.swConv([301, 701]) = 0.5; + expected_out.Evect = Evect; + out = sw_egrid(testCase.spectrum, 'Evect', Evect); + testCase.verify_obj(out, expected_out); + end + function test_Evect_cbin(testCase) + Evect_in = linspace(1.005, 2.995, 200); + expected_out = testCase.sw_egrid_out_sperp; + expected_out.swConv = zeros(200, 5); + expected_out.swConv([301, 701]) = 0.5; + expected_out.Evect = linspace(1, 3, 201); + out = sw_egrid(testCase.spectrum, 'Evect', Evect_in, 'binType', 'cbin'); + testCase.verify_obj(out, expected_out); + end + function test_temp(testCase) + temp = 300; + expected_out = testCase.sw_egrid_out_sperp; + expected_out.swConv([728, 1728]) = 6.70976913583173; + expected_out.swConv(1455) = 3.48826657260066; + expected_out.T = temp; + out = sw_egrid(testCase.spectrum, 'T', temp); + testCase.verify_obj(out, expected_out, 'rel_tol', 1e-10); + end + function test_single_ion_temp(testCase) + temp = 300; + spectrum = testCase.spectrum; + % Copy swobj here so don't interfere with other tests + spectrum.obj = copy(testCase.swobj); + spectrum.obj.single_ion.T = temp; + expected_out = testCase.sw_egrid_out_sperp; + expected_out.obj = spectrum.obj; + expected_out.swConv([728, 1728]) = 6.70976913583173; + expected_out.swConv(1455) = 3.48826657260066; + expected_out.T = temp; + out = sw_egrid(spectrum); + testCase.verify_obj(out, expected_out, 'rel_tol', 1e-10); + end + function test_twin(testCase) + swobj_twin = copy(testCase.swobj); + swobj_twin.addtwin('axis', [0 0 1], 'phid', [60 120], 'vol', [1 2]); + spectrum = swobj_twin.spinwave(testCase.spectrum.hkl); + spectrum.datestart = testCase.spectrum.datestart; + spectrum.dateend = testCase.spectrum.dateend; + out = sw_egrid(spectrum); + + expected_out = testCase.sw_egrid_out_sperp; + expected_out.obj = spectrum.obj; + expected_out.omega = spectrum.omega; + expected_out.Sab = spectrum.Sab; + expected_out.param.notwin = false; + expected_out.swConv = zeros(500, 5); + expected_out.swConv([728, 1455, 1728]) = 0.125; + expected_out.swConv([567, 1228, 1888, 2455]) = 0.65625; + expected_out.swInt = 0.78125*ones(2, 5); + expected_out.intP = cell(1,3); + expected_out.Pab = cell(1,3); + expected_out.Mab = cell(1,3); + expected_out.Sperp = {0.5*ones(2, 5) 0.875*ones(2, 5) 0.875*ones(2, 5)}; + + testCase.verify_obj(out, expected_out); + end + function test_twin_nosum(testCase) + swobj_twin = copy(testCase.swobj); + swobj_twin.addtwin('axis', [0 0 1], 'phid', [60 120], 'vol', [1 2]); + spectrum = swobj_twin.spinwave(testCase.spectrum.hkl); + spectrum.datestart = testCase.spectrum.datestart; + spectrum.dateend = testCase.spectrum.dateend; + out = sw_egrid(spectrum, 'sumtwin', false); + + expected_out = testCase.sw_egrid_out_sperp; + expected_out.obj = spectrum.obj; + expected_out.omega = spectrum.omega; + expected_out.Sab = spectrum.Sab; + expected_out.param.notwin = false; + expected_out.param.sumtwin = false; + expected_out.component = {'Sperp'}; + expected_out.swConv = cell(1, 3); + expected_out.swConv{1} = testCase.sw_egrid_out_sperp.swConv; + expected_out.swInt = cell(1, 3); + expected_out.swInt{1} = testCase.sw_egrid_out_sperp.swInt; + expected_out.intP = cell(1,3); + expected_out.Pab = cell(1,3); + expected_out.Mab = cell(1,3); + expected_out.Sperp = {0.5*ones(2, 5) 0.875*ones(2, 5) 0.875*ones(2, 5)}; + + testCase.verify_obj(out, expected_out); + end + function test_imagChk(testCase) + dE = 1; + testCase.spectrum.omega(1) = 0 + 2i*dE; % imag > dE bins + testCase.verifyError(... + @() sw_egrid(testCase.spectrum, ... + 'Evect', 0:dE:4, 'imagChk', true), ... + 'egrid:BadSolution'); + end + function test_autoEmin(testCase) + eps_imag = 1e-8; + imag_omega = 0 + 1i*eps_imag; + testCase.spectrum.omega(1) = imag_omega; + component = 'Sxx'; + expected_out = testCase.sw_egrid_out; + expected_out.component = component; + expected_out.Evect(1) = expected_out.Evect(1) + eps_imag; + expected_out.omega(1) = imag_omega; + expected_out.swConv(1) = 0; + out = sw_egrid(testCase.spectrum, 'component', component, 'autoEmin', true); + testCase.verify_obj(out, expected_out); + end + function test_modeIdx(testCase) + % only consisder -ve mode (magnon anhillation/ energy gain) + component = 'Sxx'; + expected_out = testCase.sw_egrid_out; + expected_out.component = component; + for modeIdx = 1:2 + out = sw_egrid(testCase.spectrum, 'component', component, ... + 'modeIdx', modeIdx); + if modeIdx == 2 + % only consisder -ve mode (magnon anhillation/energy gain) + % which is not in range of default energy bins + expected_out.swConv = zeros(500, 5); + end + testCase.verify_obj(out, expected_out); + end + end + function test_zeroEnergyTol(testCase) + % set zeroEnergyTol > max energy to prodcuce zero intensity + out = sw_egrid(testCase.spectrum, ... + 'binType', 'ebin' , 'zeroEnergyTol', 5); + expected_out = testCase.sw_egrid_out_sperp; + expected_out.swConv = zeros(size(expected_out.swConv)); + testCase.verify_obj(out, expected_out); + end + function test_negative_zeroEnergyTol(testCase) + out = sw_egrid(testCase.spectrum, 'zeroEnergyTol', -1); + expected_out = testCase.sw_egrid_out_sperp; + expected_out.swConv([1,2001]) = 0.5; % intensity at zero energy + testCase.verify_obj(out, expected_out); + end + function test_maxDSF(testCase) + % set maxDSF low to zero all intensity + out = sw_egrid(testCase.spectrum, ... + 'binType', 'ebin' , 'maxDSF', 1e-2); + expected_out = testCase.sw_egrid_out_sperp; + expected_out.swConv = zeros(size(expected_out.swConv)); + testCase.verify_obj(out, expected_out); + end + end + +end diff --git a/swfiles/sw_egrid.m b/swfiles/sw_egrid.m index b8e6d7743..dfea70596 100644 --- a/swfiles/sw_egrid.m +++ b/swfiles/sw_egrid.m @@ -113,8 +113,9 @@ % include all modes. % % `'epsilon'` -% : Error limit, used to determine whether a given energy bin is -% uniform or not. Default value is $10^{-5}$. +% : DEPRECATED (previously the error limit, used to determine whether a +% given energy bin is uniform or not. Default value is $10^{-5}$). This +% parameter is no longer relevant and is ignored. % % `'autoEmin'` % : Due to the finite numerical precision, the spin wave energies @@ -158,6 +159,10 @@ % `Evect` % : Input energy bin vector, defines the energy bin **edge** positions % (converted from the given bin centers if necessary). +% +% `zeroEnergyTol` +% : Eigenvalues with magnitude of the real component less than zeroEnergyTol +% will be not be included in the structure factor binning % % `param` % : All the input parameters. @@ -195,17 +200,28 @@ end inpForm.fname = {'Evect' 'T' 'component' 'sumtwin' 'modeIdx' 'epsilon'}; -inpForm.defval = {E0 T0 'Sperp' true zeros(1,0) 1e-5 }; +inpForm.defval = {E0 T0 'Sperp' true zeros(1,0) NaN }; inpForm.size = {[1 -1] [1 1] [1 -2] [1 1] [1 -4] [1 1] }; inpForm.soft = {true false false false false false }; -inpForm.fname = [inpForm.fname {'autoEmin' 'imagChk' 'binType'}]; -inpForm.defval = [inpForm.defval {false true 'ebin' }]; -inpForm.size = [inpForm.size {[1 1] [1 1] [1 -5] }]; -inpForm.soft = [inpForm.soft {false false false }]; +inpForm.fname = [inpForm.fname {'autoEmin' 'imagChk' 'binType' 'zeroEnergyTol'}]; +inpForm.defval = [inpForm.defval {false true 'ebin' 5e-4}]; +inpForm.size = [inpForm.size {[1 1] [1 1] [1 -5] [1 1]}]; +inpForm.soft = [inpForm.soft {false false false false}]; + +inpForm.fname = [inpForm.fname {'maxDSF'}]; +inpForm.defval = [inpForm.defval {1e6}]; +inpForm.size = [inpForm.size {[1 1]}]; +inpForm.soft = [inpForm.soft {false}]; param = sw_readparam(inpForm, varargin{:}); +if ~isnan(param.epsilon) + % non default value of epsilon + warning('sw_egrid:DeprecationWarning', ... + 'epsilon is deprecated - it is not relevant and will be ignored.'); +end + switch param.binType case 'ebin' eBin = true; @@ -385,7 +401,7 @@ end if ~isempty(intP{1}) - intP{ii} = reshape(intP{ii},3,3,nMode,[]); + intP{ii} = reshape(intP{ii},3,nMode,[]); end if ~isempty(Pab{1}) Pab{ii} = reshape(Pab{ii},3,3,nMode,[]); @@ -430,86 +446,40 @@ end % save the edge bins -Evect = sort(param.Evect); +param.Evect = sort(param.Evect); if eBin - eEvect = Evect; + ebin_edges = param.Evect; + ebin_cens = (ebin_edges(2:end)+ebin_edges(1:(end-1)))/2; + dE = diff(ebin_cens); else - dE = diff(Evect); - eEvect = [Evect-[dE(1) dE]/2 Evect(end)+dE(end)/2]; + ebin_cens = param.Evect; + dE = diff(ebin_cens); + ebin_edges = [ebin_cens-0.5*[dE(1) dE], ebin_cens(end) + dE(end)/2]; end +nE = numel(ebin_cens); -if isfield(spectra,'omega') - % Create vector for energy values, and put extra value below minimum and - % above maximum for easy indexing swConv. - epsilon = 1e-8; - - if eBin - Evect = (Evect(2:end)+Evect(1:(end-1)))/2; - end - % energy bin parameters - nE = numel(Evect); - dE = diff(Evect); - +if isfield(spectra,'omega') if param.imagChk % find the maximum of the imaginary part of the spin wave energies % checks only the first twin! ioMax = max(abs(imag(omega{1}(:)))); - if ioMax > max(abs(dE(:))) error('egrid:BadSolution',['The imaginary part of the spin '... 'wave energes is larger than the bin size! Improve '... 'your calculation or disable imagChk option!']); + elseif param.autoEmin + if abs(ebin_edges(1)) < ioMax + ebin_edges(1) = ebin_edges(1)+ioMax; + ebin_cens(1) = ebin_cens(1)+ioMax/2; + end end - end - - if param.autoEmin && abs(Evect(1)-dE(1)/2)=0); else - Evect = [-epsilon; epsilon]; - end - - dE = dE(1); - - % Create indices in the matrix by searching for the closest value, size - % nMode x nHkl. Put all the modes to the positive side for magnon creation. - % The negative side will be the same, however with different Bose factor - % for non-zero temperature. - idxE = cell(1,nTwin); - - - for tt = 1:nTwin - % put the modes that are not in the modeIdx parameter above the - % energy bin vector - omega{tt}(~ismember(1:nMode,param.modeIdx),:) = Evect(end); - - % faster binning - if isequalE - idxE{tt} = floor((real(omega{tt})-E0)/dE)+2; - idxE{tt}(idxE{tt}<2) = 1; - idxE{tt}(idxE{tt}>nE+2) = nE+2; - else - % memory intensive binnin, bad for large numel(omega) - [~, idxE{tt}] = min(abs(repmat(real(omega{tt}),[1 1 nE+2])-repmat(permute(Evect,[2 3 1]),[nMode nHkl 1])),[],3); - end - - % Creates indices in the swConv matrix. - idxE{tt} = idxE{tt} + repmat((0:nHkl-1).*(nE+2),[nMode 1]); - idxE{tt} = idxE{tt}(:); + nBose = 1./(exp(abs(ebin_cens)./(spectra.obj.unit.kB*param.T))-1)+double(ebin_cens>=0); end % Sums up the intensities in DSF into swConv. @@ -517,22 +487,31 @@ for tt = 1:nTwin for ii = 1:nConv - swConv{ii,tt} = reshape(accumarray(idxE{tt},DSF{ii,tt}(:),[(nE+2)*nHkl 1]),[(nE+2) nHkl]); - end - end - - % Calculate Bose temperature factor for magnons - if param.T==0 - nBose = double(Evect(2:(end-1))>=0); - else - nBose = 1./(exp(abs(Evect(2:(end-1)))./(spectra.obj.unit.kB*param.T))-1)+double(Evect(2:(end-1))>=0); - end - - % Multiply the intensities with the Bose factor. - for tt = 1:nTwin - for ii = 1:nConv - swConv{ii,tt} = bsxfun(@times,swConv{ii,tt}(2:(end-1),:),nBose); - swConv{ii,tt}(isnan(swConv{ii,tt})) = 0; + real_eigvals = real(omega{tt}(param.modeIdx, :)); + % find energy bin (cen) index coinciding with evals in omega + ien = discretize(real_eigvals, ebin_edges); + % set eigvals < zeroEnergyTol to naN (will be ignored) + izero_eigval = abs(real_eigvals(:)) < param.zeroEnergyTol; + ien(izero_eigval) = NaN; + % NaN in ien implies eigvals not in extent of Evect + ien_valid = ~isnan(ien(:)); + % get hkl index of each ien bin (column index in real_eigvals) + [~, ihkl] = ind2sub(size(ien), 1:numel(ien)); + % get index of bins in final sxConv field + % normally ien(ien_valid) is a colulmn vector but not in case + % of one modeIDx specified so need to reshape + sw_conv_idx = [reshape(ien(ien_valid), [], 1), ihkl(ien_valid)']; + if ~isempty(sw_conv_idx) + % sum intensities and pad energies above max eigval with 0 + DSF_valid = DSF{ii,tt}(param.modeIdx, :); + DSF_valid(DSF_valid > param.maxDSF) = 0; + swConv{ii,tt} = accumarray(sw_conv_idx, DSF_valid(ien_valid), [nE, nHkl]); + % Multiply the intensities with the Bose factor. + swConv{ii,tt} = bsxfun(@times,swConv{ii,tt},nBose'); + swConv{ii,tt}(isnan(swConv{ii,tt})) = 0; + else + swConv{ii,tt} = zeros(numel(ebin_cens), nHkl); + end end end @@ -565,7 +544,7 @@ spectra.T = param.T; -spectra.Evect = eEvect; +spectra.Evect = ebin_edges; if isfield(spectra,'swRaw') spectra = rmfield(spectra,'swRaw'); diff --git a/test_data/system_tests/systemstest_spinwave_KCu3As2O7.mat b/test_data/system_tests/systemstest_spinwave_KCu3As2O7.mat index d67911814..f8d3b84f7 100644 Binary files a/test_data/system_tests/systemstest_spinwave_KCu3As2O7.mat and b/test_data/system_tests/systemstest_spinwave_KCu3As2O7.mat differ diff --git a/test_data/system_tests/systemstest_spinwave_af33kagome.mat b/test_data/system_tests/systemstest_spinwave_af33kagome.mat index 52bfcfec7..489010db1 100644 Binary files a/test_data/system_tests/systemstest_spinwave_af33kagome.mat and b/test_data/system_tests/systemstest_spinwave_af33kagome.mat differ diff --git a/test_data/system_tests/systemstest_spinwave_biquadratic.mat b/test_data/system_tests/systemstest_spinwave_biquadratic.mat index 3469854b0..df8de7cc8 100644 Binary files a/test_data/system_tests/systemstest_spinwave_biquadratic.mat and b/test_data/system_tests/systemstest_spinwave_biquadratic.mat differ diff --git a/test_data/system_tests/systemstest_spinwave_pcsmo.mat b/test_data/system_tests/systemstest_spinwave_pcsmo.mat index c2a73b32d..23e0cef45 100644 Binary files a/test_data/system_tests/systemstest_spinwave_pcsmo.mat and b/test_data/system_tests/systemstest_spinwave_pcsmo.mat differ diff --git a/test_data/system_tests/systemstest_spinwave_yb2ti2o7.mat b/test_data/system_tests/systemstest_spinwave_yb2ti2o7.mat index 1e9efc863..cd51f6a8c 100644 Binary files a/test_data/system_tests/systemstest_spinwave_yb2ti2o7.mat and b/test_data/system_tests/systemstest_spinwave_yb2ti2o7.mat differ