Skip to content

Commit

Permalink
fix small bug in pepper w eigenfields
Browse files Browse the repository at this point in the history
  • Loading branch information
stestoll committed Feb 25, 2024
1 parent 1b6a329 commit 317b7e3
Show file tree
Hide file tree
Showing 7 changed files with 36 additions and 30 deletions.
2 changes: 1 addition & 1 deletion easyspin/curry.m
Original file line number Diff line number Diff line change
Expand Up @@ -286,7 +286,7 @@

Exp.R_sample = p_samplerotmatrix(Exp.SampleRotation);

[Exp,Opt] = p_symandgrid(Sys,Exp,Opt);
[Exp,Opt,nOrientations] = p_symandgrid(Sys,Exp,Opt);

% Process crystal orientations, crystal symmetry, and frame transforms
[Orientations,nOrientations,~,~] = p_crystalorientations(Exp,Opt);
Expand Down
14 changes: 10 additions & 4 deletions easyspin/pepper.m
Original file line number Diff line number Diff line change
Expand Up @@ -576,8 +576,7 @@
Opt.Intensity = anisotropicIntensities;

% Set up grid etc.
[Exp,Opt] = p_symandgrid(Sys,Exp,Opt);
nOrientations = size(Exp.SampleFrame,1);
[Exp,Opt,nOrientations] = p_symandgrid(Sys,Exp,Opt);

% Fold orientational distribution function into grid region
if partiallyOrderedSample
Expand Down Expand Up @@ -776,6 +775,7 @@
xAxis = linspace(SweepRange(1),SweepRange(2),Exp.nPoints);
Exp.deltaX = xAxis(2)-xAxis(1);


%=======================================================================
%=======================================================================
% SPECTRUM CONSTRUCTION (incl. INTERPOLATION and PROJECTION)
Expand All @@ -790,6 +790,7 @@
logmsg(1,'-absorption spectrum construction----------------------');

BruteForceSum = useEigenFields | Opt.BruteForce;

axialGrid = Opt.nOctants==0;
usingGrid = Opt.nOctants>=0;

Expand Down Expand Up @@ -885,6 +886,7 @@
msg = '(separate transitions)';
else
nSpectra = 1;
msg = '';
end
spec = zeros(nSpectra,Exp.nPoints);
logmsg(1,' spectrum array size: %dx%d %s',size(spec,1),size(spec,2),msg);
Expand Down Expand Up @@ -1121,9 +1123,13 @@

else % if Opt.ImmediateBinning elseif ~BruteForceSum ...

logmsg(1,' no interpolation',nOrientations);
logmsg(1,' constructing stick spectrum');
logmsg(1,' constructing stick spectrum (no interpolation)');
logmsg(1,' summation over %d orientations',nOrientations);

if ~isempty(Opt.separate)
error('Cannot return separate subspectra when constructing stick spectrum.');
end

spec = zeros(1,Exp.nPoints);
prefactor = (Exp.nPoints-1)/(SweepRange(2)-SweepRange(1));
for iOri = 1:nOrientations
Expand Down
2 changes: 1 addition & 1 deletion easyspin/private/p_symandgrid.m
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@
% Opt.nOctants
% Opt.GridParams (for powders)

function [Exp,Opt] = p_symandgrid(Sys,Exp,Opt)
function [Exp,Opt,nOrientations] = p_symandgrid(Sys,Exp,Opt)

logmsg(1,'-orientations------------------------------------------');

Expand Down
2 changes: 1 addition & 1 deletion easyspin/saffron.m
Original file line number Diff line number Diff line change
Expand Up @@ -301,7 +301,7 @@
% Set up orientation loop
Exp.R_sample = p_samplerotmatrix(Exp.SampleRotation);
if ~isfield(Exp,'OriWeights')
[Exp,Opt] = p_symandgrid(Sys,Exp,Opt);
[Exp,Opt,nOrienations] = p_symandgrid(Sys,Exp,Opt);
end
% Process crystal orientations, crystal symmetry, and frame transforms
[Orientations,nOrientations,nSites] = p_crystalorientations(Exp,Opt);
Expand Down
3 changes: 1 addition & 2 deletions easyspin/salt.m
Original file line number Diff line number Diff line change
Expand Up @@ -452,8 +452,7 @@
end
%====================================================================

[Exp,Opt] = p_symandgrid(Sys,Exp,Opt);
nOrientations = size(Exp.SampleFrame,1);
[Exp,Opt,nOrientations] = p_symandgrid(Sys,Exp,Opt);

% Fold orientational distribution function into grid region.
if ~isempty(Exp.Ordering)
Expand Down
25 changes: 12 additions & 13 deletions tests/resfieldseig_intensity_matrix.m → tests/pepper_eig_matrix.m
Original file line number Diff line number Diff line change
@@ -1,10 +1,9 @@
function ok = test()

%=======================================================
% resfields_eig should give the same integral intensity as matrix diagonalization
% (tested for isotropic powder and single crystal)
% resfields_eig should give the same integral intensity as resfields
% (tested using pepper simulation of isotropic powder and single crystal)
%=======================================================
clear Sys Exp

Sys.g = 2.1;
Sys.Nucs = '1H';
Expand All @@ -17,25 +16,25 @@

% isotropic powder
Opt.Method = 'matrix';
[x,y1] = pepper(Sys,Exp,Opt);
[B,spc1] = pepper(Sys,Exp,Opt);
Opt.Method = 'eig';
[x,y2] = pepper(Sys,Exp,Opt);
dx = x(2)-x(1);
[B,spc2] = pepper(Sys,Exp,Opt);
dx = B(2)-B(1);

integral_matrix = sum(y1)*dx;
integral_resfields_eig = sum(y2)*dx;
integral_matrix = sum(spc1)*dx;
integral_resfields_eig = sum(spc2)*dx;

ok(1) = areequal(integral_matrix,integral_resfields_eig,0.001,'abs');

% single crystal
Exp.SampleFrame = rand(1,3)*pi;
Opt.Method = 'matrix';
[x,y1] = pepper(Sys,Exp,Opt);
[B,spc1] = pepper(Sys,Exp,Opt);
Opt.Method = 'eig';
[x,y2] = pepper(Sys,Exp,Opt);
dx = x(2)-x(1);
[B,spc2] = pepper(Sys,Exp,Opt);
dx = B(2)-B(1);

integral_matrix = sum(y1)*dx;
integral_resfields_eig = sum(y2)*dx;
integral_matrix = sum(spc1)*dx;
integral_resfields_eig = sum(spc2)*dx;

ok(2) = areequal(integral_matrix,integral_resfields_eig,0.001,'abs');
18 changes: 10 additions & 8 deletions tests/salt_pt_threehalf_onehalf.m
Original file line number Diff line number Diff line change
@@ -1,25 +1,27 @@
function ok = test(opt)

%-------------------------------------------------------------
%S=3/2, I=1/2
% ENDOR spectrum of a S=3/2, I=1/2 system
%-------------------------------------------------------------

Sys.S = 3/2;
Sys.D = 200;
Sys.Nucs='1H';
Sys.A = [3 5];
Sys.lwEndor = 0.2;

Exp.Field = 350;
Exp.mwFreq = 9.5;
Exp.Range = [0 30];
Opt.Method='matrix';
[x,a]=salt(Sys,Exp,Opt);
Opt.Method='perturb2';
[x,b]=salt(Sys,Exp,Opt);
b=rescaledata(b,a,'minmax');

Opt.Method = 'matrix';
[nu,spc_matrix] = salt(Sys,Exp,Opt);
Opt.Method = 'perturb2';
[nu,spc_pt2] = salt(Sys,Exp,Opt);

if opt.Display
plot(x,a,'k',x,b,'r');
plot(nu,spc_matrix,nu,spc_pt2);
legend('matrix','perturb2');
end

ok = areequal(b,a,0.01,'rel');
ok = areequal(spc_pt2,spc_matrix,0.01,'rel');

0 comments on commit 317b7e3

Please sign in to comment.