diff --git a/pathdef.m b/pathdef.m
index 50a86c112d3ca13309955f9032724efbcbea9ede..1beb498a6d76e89d574e87b76b24b0f22b4dd56b 100755
--- a/pathdef.m
+++ b/pathdef.m
@@ -122,6 +122,7 @@ p = [...
      'S:\soft\matlab\easyspin-6.0.0\documentation\highlight;', ...
      'S:\soft\matlab\easyspin-6.0.0\documentation\img;', ...
      'S:\soft\matlab\easyspin-6.0.0\easyspin;', ...
+     'S:\soft\matlab\easyspin-6.0.0\easyspin\util;', ...
      'S:\soft\matlab\easyspin-6.0.0\examples;', ...
      'S:\soft\matlab\easyspin-6.0.0\examples\analysis;', ...
      'S:\soft\matlab\easyspin-6.0.0\examples\endor;', ...
@@ -146,7 +147,7 @@ p = [...
      'S:\soft\matlab\my_util\cw;', ...
      'S:\soft\matlab\my_util\pulse;', ...
      'S:\soft\matlab\my_util\scr;', ...
-     'C:\Users\gianlum33\AppData\Local\Temp\111\Editor_lcxdh;', ...
+     'C:\Users\gianlum33\AppData\Local\Temp\27\Editor_cwocb;', ...
      matlabroot,'\toolbox\matlab\capabilities;', ...
      matlabroot,'\toolbox\matlab\datafun;', ...
      matlabroot,'\toolbox\matlab\datatypes;', ...
diff --git a/reports/oop_ciss_powder_average_zech.tex b/reports/oop_ciss_powder_average_zech.tex
index 0defddfd1a56e8431ca3dc7e00bf4f5c967ce638..08fd71bc58b56cf41b7339d96085b8fe567e5d46 100644
--- a/reports/oop_ciss_powder_average_zech.tex
+++ b/reports/oop_ciss_powder_average_zech.tex
@@ -30,6 +30,6 @@
 		-\frac{1}{2}\sin[2(d - J)\tau]\left[\frac{1}{2}\sin{(\beta)}\sin^2{(2\xi)} + \sin(2\beta) \cos^4(\xi)\right],
 	\label{eq:oop-eseem_singlet}
 	\end{equation}
-	where $d$ and $J$ are the dipolar and exchange coupling respectively, $\beta$ is the turning angle of the first microwave pulse of the sequence and $\xi$ is the mixing angle. In particular, we define $\xi = \arctan{[(d + 2J)/(\Omega_A - \Omega_B)]}$, where $\Omega_{j} = g_{j}\frac{\mu_B B_0}{\hbar} - \omega_{MW}$. This definition implies that in the limit of weakly interacting spins $\xi \rightarrow 0$ and $E_x \rightarrow \sin(2\beta)$, while in the limit of strongly coupled spins $\xi \rightarrow \pi/2$ and $E_x \rightarrow \sin(\beta)4\xi^2$, therefore approaches $\sin(\beta)$ but vanishes at the limit. 
+	where $d$ and $J$ are the dipolar and exchange coupling respectively, $\beta$ is the turning angle of the first microwave pulse of the sequence and $\xi$ is the mixing angle. In particular, we define $\xi = \arctan{[(d + 2J)/(\Omega_A - \Omega_B)]}$, where $\Omega_{j} = g_{j}\frac{\mu_B B_0}{\hbar} - \omega_{MW}$. This definition implies that in the limit of weakly interacting spins $\xi \rightarrow 0$ and $E_x \rightarrow \sin(2\beta)$, while in the limit of strongly coupled spins $\xi \rightarrow \pi/2$ and $E_x \rightarrow \sin(\beta)\xi^2$, therefore approaches $\sin(\beta)$ but vanishes at the limit. 
 	\printbibliography
 \end{document}
\ No newline at end of file
diff --git a/zech_psiOopEseem_powderAverage_inhBroadening.asv b/zech_psiOopEseem_powderAverage_inhBroadening.asv
deleted file mode 100755
index 8cdef813785bae3b8650002400787344e73ef286..0000000000000000000000000000000000000000
--- a/zech_psiOopEseem_powderAverage_inhBroadening.asv
+++ /dev/null
@@ -1,544 +0,0 @@
-%
-clearvars
-
-% Import
-% cd('/net/storage/gianlum33/projects/oop_ciss_calculations/')
-expName = 'data_extracted/zech_p46_oopEseem_expData.csv';
-fitName = 'data_extracted/zech_p46_oopEseem_fit.csv';
-
-importedData = readtable(expName);
-x = importedData{:, 1}*pi/180;
-y = importedData{:, 2};
-importedData = readtable(fitName);
-xx = importedData{:, 1}*pi/180;
-yy = importedData{:, 2};
-
-crystalEseem = @(p, beta) -( ...
-    0.5*sin(2*beta) * sin(2*p)^4 + ...
-    sin(beta) * cos(2*p)^2 * sin(2*p)^2);  % p = alpha
-
-p0Z = 31.5*pi/180;
-p0Timmel = 41*pi/180;
-initFit = rescaledata(crystalEseem(p0Z, xx), 'maxabs');
-expectFit = rescaledata(crystalEseem(p0Timmel, xx), 'maxabs');
-
-crystalEseem = @(p, beta) -(...
-    1/2 * sin(2*beta) .* cos(p).^4 + ...
-    1/4 * sin(beta) .* sin(2*p).^2);  % p = 2alpha
-
-p0B = pi/2 - 2*p0Z;  % xi_Bittl = pi/2 - 2*alpha_Zech
-initFit2 = rescaledata(crystalEseem(p0B, xx), 'maxabs');
-
-clf
-plot(x*180/pi, rescaledata(y, 'maxabs'), 'o-')
-hold on
-plot(xx*180/pi, rescaledata(yy, 'maxabs'))
-plot(xx*180/pi, initFit)
-plot(xx*180/pi, initFit2)
-plot(xx*180/pi, expectFit)
-% esfit()
-
-%% Dipolar interaction matrix, g1 and g2
-
-% Dipolar interaction, p = [dip, theta, phi]
-dFunc = @(p) p(1)*((sin(p(2))*cos(p(3)))^2 - 1/3);
-
-% Mixing angle, p = [dd, JJ, B0, g1, g2] in [MHz, MHz, T, adim, adim]
-xiFunc = @(p) atan( (p(1) + 2*p(2))*1e6 * planck/(2*pi*(bmagn*p(3))) ./ (p(4) - p(5)) );
-
-% Rotation matrices
-Rz = @(x) [cos(x) sin(x)    0; 
-            -sin(x) cos(x)  0;
-            0       0       1];
-Ry = @(x) [cos(x)   0   -sin(x); 
-           0        1   0;
-           sin(x)   0   cos(x)];
-rotMatrixZyz = @(euAng) Rz(euAng(3))*Ry(euAng(2))*Rz(euAng(1));
-
-% Parameters of the system
-eulerAnglesZech = [81, 126, 182]*pi/180;  % Zech, 'A structural model ...'
-eulerAnglesZech = [-182, 126, 81]*pi/180;  % Zech, 'A structural model ...'
-
-Sys.S = [1/2 1/2];  % Taken from Easyspin example
-Sys.g = [2.0030 2.0026 2.0023; ... % P700+
-         2.0062 2.0051 2.0022];    % A1-
-Sys.gFrame = [-10 -128 -83; ...
-                0    0   0]*pi/180;
-Sys.eeFrame = [0 90 0]*pi/180;  % zD direction is along -x of A1-
-eulerAnglesEasyspin = Sys.gFrame(1, :);
-% Sys.J = -unitconvert(1e-3,'mT->MHz'); % MHz
-% Sys.dip = unitconvert(+0.177,'mT->MHz'); % MHz
-
-dip = unitconvert(0.177,'mT->MHz');
-JJ = unitconvert(-1e-3,'mT->MHz');
-mwFreq = 9.7e9;  % Hz
-B0 = 0.3459;  % T, static magnetic field (max of trEPR spectrum)
-
-thetas = 0:0.5:180;
-nTheta = numel(thetas);
-phis = 0:0.5:360;
-nPhi = numel(phis);
-[g1, g2, dd, cosThetaD, xiNoLw] = deal(zeros(nTheta, nPhi));
-[nVers, nVersInFrame1] = deal(repmat({[0, 0, 0]}, nTheta, nPhi));
-signalPowderNoLw = zeros(numel(xx), 1);
-signalNoLw_ = repmat({signalPowderNoLw}, nTheta, nPhi);
-
-tic
-for ith = 1:nTheta
-    % Update waitbar and message
-    if ~exist('wbarTh', 'var')
-        wbarTh = waitbar(0, '1', 'Name', 'Theta');
-    end
-    waitbar(ith/nTheta, wbarTh, sprintf('%d out of %d', ith, nTheta))
-
-    theta_ = thetas(ith)*pi/180;
-    for iph = 1:nPhi
-        phi_ = phis(iph)*pi/180;
-        
-        % Direction of B0 in the frame of spin 2 and g_eff of spin 2
-        nVers{ith, iph} = ...
-            [cos(phi_)*sin(theta_), sin(phi_)*sin(theta_), cos(theta_)];
-        g2(ith, iph) = sqrt( sum( (Sys.g(2, :).*nVers{ith, iph}).^2));
-
-        % Angle between zD and B0, the sign is not important
-        cosThetaD(ith, iph) = -sin(theta_)*cos(phi_);
-        dd(ith, iph) = dFunc([dip, theta_, phi_]);
-        
-        % Direction of B0 in the frame of spin 1 and g_eff of spin 1
-        nVersInFrame1{ith, iph} = ...
-            rotMatrixZyz(eulerAnglesEasyspin)*nVers{ith, iph}';
-        g1(ith, iph) = sqrt( sum( (Sys.g(1, :).*nVersInFrame1{ith, iph}').^2));
-
-        xiNoLw(ith, iph) = ...
-            xiFunc([dd(ith, iph), JJ, B0, g1(ith, iph), g2(ith, iph)]);
-        
-        % Calculate signal
-        signalNoLw_{ith, iph} = crystalEseem(xiNoLw(ith, iph), xx);
-        theta_ = thetas(ith)*pi/180;  % Solid angle normalization
-        signalNoLw_{ith, iph} = signalNoLw_{ith, iph}*sin(theta_);
-        signalPowderNoLw = signalPowderNoLw + signalNoLw_{ith, iph};
-    end
-end
-toc
-
-delete(wbarTh)
-clear('wbarTh')        
-        
-%%
-
-% TODO maybe set same colorbar extremes for both
-figure(1)
-clf
-tL = tiledlayout(1, 2);
-nexttile
-% imagesc(phis, thetas, g1)
-imagesc(g1)
-colorbar
-title('g1')
-nexttile
-% imagesc(phis, thetas, g2)
-imagesc(g2)
-colorbar
-labelaxesfig(tL, 'Phi', 'Theta')
-title('g2')
-
-figure(2)
-clf
-% imagesc(dip*(valueOfCosThetaD.^2 - 1/3))
-% imagesc(thetas, phis, dd)
-% imagesc(phis, thetas, dd)
-imagesc(dd)
-labelaxesfig(gca, 'Phi', 'Theta')
-% imagesc(thetas, phis, xiNoLw)
-cbar = colorbar;
-% set(cbar, 'YDir', 'reverse');
-title('Dipolar interaction dd')
-
-
-figure(3)
-clf
-% imagesc(phis, thetas, g1 - g2)
-imagesc(g1 - g2)
-labelaxesfig(gca, 'Phi', 'Theta')
-cbar = colorbar;
-title('g1 minus g2')
-
-figure(4)
-clf
-% imagesc(dip*(valueOfCosThetaD.^2 - 1/3))
-% imagesc(phis, thetas, xiNoLw)
-imagesc(xiNoLw)
-cbar = colorbar;
-labelaxesfig(gca, 'Phi', 'Theta')
-% set(cbar, 'YDir', 'reverse');
-title('xiNoLw')
-
-%
-figure(5)
-clf
-plot(xx, signalNoLw_{22, 200}, 'DisplayName', string(xiNoLw(22, 200)*180/pi) + " deg")
-hold on
-plot(xx, signalNoLw_{45, 100}, 'DisplayName', string(xiNoLw(45, 100)*180/pi) + " deg")
-plot(xx, signalNoLw_{46, 200}, 'DisplayName', string(xiNoLw(46, 200)*180/pi) + " deg")
-plot(xx, signalNoLw_{43, 200}, 'DisplayName', string(xiNoLw(43, 200)*180/pi) + " deg")
-plot(xx, signalNoLw_{335, 540}, 'DisplayName', string(xiNoLw(335, 540)*180/pi) + " deg")
-legend()
-
-%% Optimized for loop
-
-% Inhomogenous line broadening due to hfi (isotropic)
-% Prisner et al., Time-resolved W-band..., 1995
-lw1 = 10e6;  % Hz
-lw2 = 10e6;  % Hz
-% The g-factor linewidth
-glw1 = planck/bmagn/B0*lw1;
-glw2 = planck/bmagn/B0*lw2;
-% Spacing between g-values
-dgax = 1e-5;
-
-% pulseBw = 0.0005;  % T
-% gPulseBw = planck*mwFreq/bmagn*(pulseBw/B0/B0);
-% gPulsePosition = planck*mwFreq/bmagn/B0;
-% excitRange = [-1/2, +1/2]*gPulseBw + gPulsePosition;
-
-signalPowder = zeros(numel(xx), 1);
-signal_ = repmat({signalPowder}, nTheta, nPhi);
-xi = zeros(nTheta, nPhi);
-
-tic
-for ith = 30
-    % Update waitbar and message
-    if ~exist('wbarTh', 'var')
-        wbarTh = waitbar(0, '1', 'Name', 'Theta');
-    end
-    waitbar(ith/nTheta, wbarTh, sprintf('%d out of %d', ith, nTheta))
-
-    theta_ = thetas(ith)*pi/180;
-    % theta_ = theta_/multTheta;
-    % theta_ = theta_ * 10;
-    for iph = 70
-        phi_ = phis(iph)*pi/180;
-        
-        % Finite pulse bandwidth: no signal from spins outside excitRange
-        %{
-        if g1(ith, iph) < excitRange(1) || ...
-            g1(ith, iph) > excitRange(2) || ... 
-            g2(ith, iph) < excitRange(1) || ...
-            g2(ith, iph) > excitRange(2)
-            continue
-        end
-        %}
-        
-        % Consider inhomogenous broadening
-        % Create g-axis in order to 'sample' the inh. broadened g-values
-        gaxLarge = creategaxis(g1(ith, iph), g2(ith, iph), glw1, glw2, dgax);
-        % Inhomogeneously broadened gaussian distributions
-        gInh1Large = gaussian(gaxLarge, g1(ith, iph), glw1);
-        gInh2Large = gaussian(gaxLarge, g2(ith, iph), glw2);
-        % Normalization
-        gInh1Large = gInh1Large/sum(gInh1Large);
-        gInh2Large = gInh2Large/sum(gInh2Large);
-        
-        [gInh1, gInh2, gax] = ...
-            cutinhdistributions(gInh1Large, gInh2Large, gaxLarge);
-        
-        ngax = numel(gax);
-        
-        xi_ = zeros(1, ngax);
-        dd_ = dd(ith, iph);
-        signalInh = zeros(numel(xx), ngax);
-        gaxFirstValue = gax(1);
-
-        gaussianFactors = gInh1'*gInh2;
-        signalInhFactors = zeros(1, ngax);
-        
-        xi_ = atan( (dd_ + 2*JJ)*1e6 * ...
-                planck/(2*pi*(bmagn*0.35)) ./ (gax - gaxFirstValue) );
-        signalInh = crystalEseem(xi_, xx);
-        parfor ig = 1:ngax
-            % Calculate xi at fixed g2 = gax(1) value
-%             xi_(ig) = atan( (dd_ + 2*JJ)*1e6 * ...
-%                 planck/(2*pi*(bmagn*0.35)) ./ (gax(ig) - gaxFirstValue) );
-%             signalInh(:, ig) = crystalEseem(xi_(ig), xx);
-            
-            % Calculate multiplicity
-%             if ig ~= 1
-%                 iDiag = ig - 1;
-%                 signalInhFactors(ig) = sum(spdiags(gaussianFactors, -iDiag)) + ...
-%                 sum(spdiags(gaussianFactors, iDiag));
-%             else  % ig == 1
-%                 iDiag = ig - 1;  % = 0
-%                 signalInhFactors(ig) = sum(spdiags(gaussianFactors, iDiag));
-%             end
-            iDiag = ig - 1;
-            signalInhFactors(ig) = sum(spdiags(gaussianFactors, -iDiag)) + ...
-                sum(spdiags(gaussianFactors, iDiag));
-        end
-        
-        xi(ith, iph) = mean(xi_.*signalInhFactors);
-        signalInh = signalInh.*signalInhFactors;
-        signal_{ith, iph} = sum(signalInh, 2)*sin(theta_);
-    end
-
-    signalPowder = signalPowder + signal_{ith, iph};
-end
-toc
-
-% save('powderAverage_inhBroadening_EasyspinAngles.mat', 'xi', 'signal_')
-
-delete(wbarTh)
-clear('wbarTh')
-
-clf
-plot(xx, signal_{ith, iph})
-hold on
-plot(xx, signalNoLw_{ith, iph})
-
-max(signal_{ith, iph})
-%% Test
-
-ith = 30;
-iph = 400;
-
-gax = creategaxis(g1(ith, iph), g2(ith, iph), glw1, glw2, dgax);
-ngax = numel(gax);
-% Inhomogeneously broadened gaussian distributions
-gInh1 = gaussian(gax, g1(ith, iph), glw1);
-gInh2 = gaussian(gax, g2(ith, iph), glw2);
-% Normalization (?)
-gInh1 = gInh1/sum(gInh1);
-gInh2 = gInh2/sum(gInh2);
-
-diff = zeros(ngax, ngax);
-for i1 = 1:ngax
-    for i2 = 1:ngax
-        diff(i1, i2) = gInh1(i1) - gInh2(i2);
-    end
-end
-
-g1(ith, iph) - g2(ith, iph)
-mean(abs(diff(:)))
-
-
-clf
-imagesc(diff)
-colorbar
-% clf
-% plot(gax, gInh1, gax, gInh2)
-% 
-% nStep = round(2*max(diff(:))/dgax);
-% diffDistr = zeros(nStep, 1);
-% for i1 = 1:ngax
-%     for i2 = 1:ngax
-%         aa = round(diff(i1, i2)/dgax);
-%     end
-% end
-
-%%
-
-size(crystalEseem(xi_(1:10)', xx))
-
-%% Some plots
-
-ith = 1;
-iph = 1;
-fig = figure(6);
-tiledlayout(2, 2, 'TileSpacing', 'compact', 'Padding', 'tight')
-nexttile
-%contour(xi{ith, iph})
-colorbar
-
-nexttile
-% plot(1:numel(xi{ith, iph}(:, 100)), xi{ith, iph}(:, 100))
-% textStr = sprintf('Max = %.3f, Min = %.3f', ...
-%     max(xi{ith, iph}(:, 100)), min(xi{ith, iph}(:, 100)));
-% text(10, -1.5, textStr)
-
-nexttile
-%%
-ith = 30;
-iph = 59;
-gax = creategaxis(g1(ith, iph), g2(ith, iph), glw1, glw2, dgax);
-gInh1 = gaussian(gax, g1(ith, iph), glw1);
-gInh2 = gaussian(gax, g2(ith, iph), glw2);
-gInh1 = gInh1/sum(gInh1);
-gInh2 = gInh2/sum(gInh2);
-contour(gInh1'*gInh2)
-colorbar
-
-clf
-plot(gax, gInh1, '-o')
-hold on
-plot(gax, gInh2, '-o')
-yyaxis right
-plot(gax, gInh1.*gInh2, '-o')
-
-max(gInh1)
-max(gInh2)
-max(gInh1.*gInh2)
-%{
-nexttile(5, [2, 2])
-plot(x*180/pi, rescaledata(y, 'maxabs'), 'o-')
-hold on
-plot(xx*180/pi, rescaledata(yy, 'maxabs'))
-plot(xx*180/pi, initFit)
-plot(xx*180/pi, initFit2)
-plot(xx*180/pi, rescaledata(signalPowder, 'maxabs'), 'ko-')
-textStr = sprintf('Max = %.3f, Min = %.3f', ...
-    max(rescaledata(signalPowder, 'maxabs')), ...
-    min(rescaledata(signalPowder, 'maxabs')));
-text(10, 0.5, textStr)
-
-textStr = "thetas = 0:0.5:180; phis = 0:0.5:360;" + newline + ...
-    "dgax = 0.00005; computational time 5.16 h";
-text(110, -0.5, textStr)
-%}
-
-%% 
-
-%{
-% [nTheta, nPhi] = size(xi);
-signalPowder = zeros(numel(xx), 1);
-
-wbarTh = waitbar(0, '1', 'Name', 'signalPowder');
-% wbarPh = waitbar(0, '1', 'Name', 'In the for cycle (Phi)');
-% 
-for ith = 1:nTheta
-    % Update waitbar and message
-    waitbar(ith/nTheta, wbarTh, sprintf('%d out of %d', ith, nTheta))
-    
-    for iph = 1:nPhi
-%         waitbar(iph/nPhi, wbarPh, sprintf('%d out of %d', iph, nPhi))
-        xi_ = xi{ith, iph};
-        [ngax, ~] = size(xi_);
-        signalSummed = zeros(numel(xx), ngax);
-        for ig1 = 1:ngax
-            for ig2 = 1:ngax
-                signalSummed(:, ig1) = signalSummed(:, ig1) + ...
-                    crystalEseem(xi_(ig1, ig2), xx);
-            end
-            signal_{ith, iph} = sum(signalSummed, 2);
-        end
-        % Solid angle normalization
-%         theta_ = thetas(ith)*pi/180;
-%         signal_{ith, iph} = sum(signalSingleSumg2, 2)*sin(theta_);
-%         signalPowder = signalPowder + signal_{ith, iph};
-    end
-    theta_ = thetas(ith)*pi/180;
-    signalPowder = signalPowder + signal_{ith, iph}*sin(theta_);
-end
-
-delete(wbarTh)
-% delete(wbarPh)
-%}
-
-% signalPowder = zeros(numel(xx), 1);
-% for ith = 1:nTheta
-%     theta_ = thetas(ith)*pi/180;
-%     for iph = 1:nPhi
-%         signalPowder = signalPowder + signal_{ith, iph}*sin(theta_);
-%     end
-% end
-
-figure(6)
-clf
-plot(x*180/pi, rescaledata(y, 'maxabs'), 'o-')
-hold on
-plot(xx*180/pi, rescaledata(yy, 'maxabs'))
-plot(xx*180/pi, initFit)
-plot(xx*180/pi, initFit2)
-plot(xx*180/pi, expectFit)
-plot(xx*180/pi, rescaledata(signalPowder, 'maxabs'), 'k-')
-legend('Exp. data', 'Fit Zech', 'Fit Gianluca using alpha', ...
-    'Fit Gianluca using xi', 'Sum of ESEEM signals', ...
-    'Location', 'northwest')
-labelaxesfig(gca, 'Flip angle beta', 'ESEEM intensity')
-% plot(xx*180/pi, -rescaledata(sin(2*xx), y))
-
-%%
-aa.signalPowder = zeros(numel(xx), 1);
-
-for ith = 1:nTheta
-    for iph = 1:nPhi
-        aa.signalPowder = aa.signalPowder + aa.signal_{ith, iph};
-    end
-end
-%% Mean xi
-
-% for ith = 1:nTheta
-%     theta_ = thetas(ith)*pi/180;
-%     xitest(ith, :) = repmat(cos(theta_)^4, 1, nPhi);
-%     xitest(ith, :) = repmat(1, 1, nPhi);
-% end
-
-for ith = 1:nTheta
-    theta_ = thetas(ith)*pi/180;
-    xiMean = xiMean + sum(xi(ith, :))*sin(theta_);
-end
-normTh = pi/2/nTheta;
-normPh = 1/nPhi;
-normTot = normPh*normTh;
-xiMean = xiMean*normTot*180/pi % /nTheta*pi/2
-% xiMean = xiMean/normTh
-
-%%
-
-g1_ = g1(1, 1);
-g2_ = g2(1, 1);
-lw1 = 15e6;  % Hz
-lw2 = 15e6;  % Hz
-% The g-factor linewidth for each one is (delta_nu_i = lw_i):
-% delta_g_i mu_B B0 = h delta_nu_i
-glw1 = planck/bmagn/0.35*lw1;
-glw2 = planck/bmagn/0.35*lw2;
-
-xx = creategaxis(g1_, g2_, glw1, glw2, 0.0001);
-% xx = 1.995:0.0001:2.015;
-
-gd1 = gaussian(xx, g1_, glw1);
-gd2 = gaussian(xx, g2_, glw2);
-
-figure()
-plot(xx, gd1, xx, gd2)
-
-%%
-
-clf
-aa = creategaxis(g1(140, 1), g2(140, 1), glw1, glw2, dgax);
-plot(aa, gaussian(aa, g1(140, 1), glw1), ...
-    aa, gaussian(aa, g2(140, 1), glw2))
-yyaxis right
-plot(aa, gaussian(aa, g1(140, 1), glw1) .* gaussian(aa, g2(140, 1), glw2))
-
-%%
-
-function xx = creategaxis(g1, g2, glw1, glw2, dg)
-    xxFactor = 1;
-    xmin = min([g1 - xxFactor*glw1, g2 - xxFactor*glw2]);
-    xmax = max([g1 + xxFactor*glw1, g2 + xxFactor*glw2]);
-    xx = xmin:dg:xmax;
-end
-
-function [gInh1, gInh2, gax] = cutinhdistributions(gInh1Large, gInh2Large, gaxLarge)
-    xxFactor = 2;
-    multipl = gInh1Large.*gInh2Large;
-    [maxMult, idxMax] = max(multipl);
-    [~, idxFwhm] = min(abs(multipl - maxMult/2));
-    leftIdx = idxMax - xxFactor*abs(idxMax - idxFwhm);
-    rightIdx = idxMax + xxFactor*abs(idxMax - idxFwhm);
-    gInh1 = gInh1Large(leftIdx:rightIdx);
-    gInh2 = gInh2Large(leftIdx:rightIdx);
-    gax = gaxLarge(leftIdx:rightIdx);
-end
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/zech_psiOopEseem_powderAverage_inhBroadening.m b/zech_psiOopEseem_powderAverage_inhBroadening.m
index e7fefc9fffc8d9ac086962e0ed433aaad83b7917..45bb8a3450cab45612b84b239aceeaac547abca6 100755
--- a/zech_psiOopEseem_powderAverage_inhBroadening.m
+++ b/zech_psiOopEseem_powderAverage_inhBroadening.m
@@ -42,9 +42,11 @@ plot(xx*180/pi, expectFit)
 
 % Dipolar interaction, p = [dip, theta, phi]
 dFunc = @(p) p(1)*((sin(p(2))*cos(p(3)))^2 - 1/3);
+dipFunc = @(dip, theta, phi) dip*((sin(theta)*cos(phi)).^2 - 1/3);
 
 % Mixing angle, p = [dd, JJ, B0, g1, g2] in [MHz, MHz, T, adim, adim]
-xiFunc = @(p) atan( (p(1) + 2*p(2))*1e6 * planck/(2*pi*(bmagn*p(3))) ./ (p(4) - p(5)) );
+xiFunc = @(dd, JJ, B0, g1, g2) atan( (dd + JJ).*1e6 .* ...
+    planck./(2*pi*(bmagn*B0)) ./ (g1 - g2) );
 
 % Rotation matrices
 Rz = @(x) [cos(x) sin(x)    0; 
@@ -73,57 +75,72 @@ dip = unitconvert(0.177,'mT->MHz');
 JJ = unitconvert(-1e-3,'mT->MHz');
 mwFreq = 9.7e9;  % Hz
 B0 = 0.3459;  % T, static magnetic field (max of trEPR spectrum)
+% B0 = 0.3461;
 
-thetas = 0:0.5:180;
+pulseBw = 0.5e-3;  % T
+gPulseBw = planck*mwFreq/bmagn*(pulseBw/B0/B0);
+gPulsePosition = planck*mwFreq/bmagn/B0;
+excitRange = [-1/2, +1/2]*gPulseBw + gPulsePosition;
+
+thetas = (0:0.05:180)*pi/180;
 nTheta = numel(thetas);
-phis = 0:0.5:360;
+phis = (0:0.05:360)*pi/180;
 nPhi = numel(phis);
-[g1, g2, dd, cosThetaD, xiNoLw] = deal(zeros(nTheta, nPhi));
-[nVers, nVersInFrame1] = deal(repmat({[0, 0, 0]}, nTheta, nPhi));
-signalPowderNoLw = zeros(numel(xx), 1);
-signalNoLw_ = repmat({signalPowderNoLw}, nTheta, nPhi);
 
-tic
-for ith = 1:nTheta
-    % Update waitbar and message
-    if ~exist('wbarTh', 'var')
-        wbarTh = waitbar(0, '1', 'Name', 'Theta');
-    end
-    waitbar(ith/nTheta, wbarTh, sprintf('%d out of %d', ith, nTheta))
+clear('dd', 'nVers', 'nVersInFrame1', 'g1', 'g2', 'xiNoLw', 'xiNoLw3D', ...
+    'signalNoLw', 'signalPowderNoLw')
+dd = dipFunc(dip, thetas', phis);
+nVers(1, :, :) = sin(thetas')*cos(phis);
+nVers(2, :, :) = sin(thetas')*sin(phis);
+nVers(3, :, :) = cos(thetas')*ones(1, nPhi);
+g2 = squeeze(sqrt( sum( (Sys.g(2, :)'.*nVers).^2)));
+eulerRotMatrix = rotMatrixZyz(eulerAnglesEasyspin);
+nVersInFrame1 = pagemtimes(eulerRotMatrix, nVers);
+g1 = squeeze(sqrt( sum( (Sys.g(1, :)'.*nVersInFrame1).^2)));
 
-    theta_ = thetas(ith)*pi/180;
-    for iph = 1:nPhi
-        phi_ = phis(iph)*pi/180;
-        
-        % Direction of B0 in the frame of spin 2 and g_eff of spin 2
-        nVers{ith, iph} = ...
-            [cos(phi_)*sin(theta_), sin(phi_)*sin(theta_), cos(theta_)];
-        g2(ith, iph) = sqrt( sum( (Sys.g(2, :).*nVers{ith, iph}).^2));
-
-        % Angle between zD and B0, the sign is not important
-        cosThetaD(ith, iph) = -sin(theta_)*cos(phi_);
-        dd(ith, iph) = dFunc([dip, theta_, phi_]);
+tic
+% Mixing angles
+xiNoLw = xiFunc(dd, JJ, B0, g1, g2);
+% signal
+xiNoLw3D(1, :, :) = xiNoLw;
+signalNoLw = crystalEseem(xiNoLw3D, xx);
+% Solid angle normalization
+xiNoLwMean = ...
+    sum(sin(thetas').*abs(xiNoLw), 'all')/sum(sin(thetas')*ones(1, nPhi), 'all');
+signalNoLw = sin(thetas).*signalNoLw;
+% Sum
+signalPowderNoLw = sum(squeeze(sum(signalNoLw, 3)), 2);
+toc     
         
-        % Direction of B0 in the frame of spin 1 and g_eff of spin 1
-        nVersInFrame1{ith, iph} = ...
-            rotMatrixZyz(eulerAnglesEasyspin)*nVers{ith, iph}';
-        g1(ith, iph) = sqrt( sum( (Sys.g(1, :).*nVersInFrame1{ith, iph}').^2));
+%%
+figure(6)
+mySignal = rescaledata(signalPowderNoLw, 'maxabs');
+testManualFit = rescaledata(crystalEseem(6.18*pi/180, xx), 'maxabs');
+crystalMean = rescaledata(crystalEseem(xiNoLwMean, xx), 'maxabs');
 
-        xiNoLw(ith, iph) = ...
-            xiFunc([dd(ith, iph), JJ, B0, g1(ith, iph), g2(ith, iph)]);
-        
-        % Calculate signal
-        signalNoLw_{ith, iph} = crystalEseem(xiNoLw(ith, iph), xx);
-        theta_ = thetas(ith)*pi/180;  % Solid angle normalization
-        signalNoLw_{ith, iph} = signalNoLw_{ith, iph}*sin(theta_);
-        signalPowderNoLw = signalPowderNoLw + signalNoLw_{ith, iph};
-    end
-end
-toc
+clf
+% tL = tiledlayout(2, 1, 'TileSpacing', 'compact', 'Padding', 'Compact');
+% nexttile
+plot(x*180/pi, rescaledata(y, 'maxabs'), 'o-', 'DisplayName', 'Zech exp. data')
+hold on
+plot(xx*180/pi, rescaledata(yy, 'maxabs'), 'DisplayName', 'Zech Fit')
+% plot(xx*180/pi, initFit)
+% plot(xx*180/pi, initFit2)
+plot(xx*180/pi, mySignal, 'DisplayName', 'Gianluca powder')
+plot(xx*180/pi, expectFit, 'k', 'DisplayName', '8 degrees')
+plot(xx*180/pi, testManualFit, '--', 'DisplayName', '6.18 deg test')
+plot(xx*180/pi, crystalMean, 'DisplayName', 'X-tal of avg xi (4deg)')
+xlim(setaxlim(xx*180/pi, 1))
+
+% nexttile
+% plot(xx*180/pi, mySignal - expectFit)
+legend('Location', 'northwest')
+labelaxesfig(gca, 'Flip angle beta', 'ESEEM intensity')
 
-delete(wbarTh)
-clear('wbarTh')        
-        
+% [temp01, temp1] = min(abs(expectFit(5:end - 1)));
+% temp1 = 4 + temp1;
+% xx(temp1)*180/pi
+% mySignal(temp1)
 %%
 
 % TODO maybe set same colorbar extremes for both
@@ -188,21 +205,16 @@ legend()
 
 % Inhomogenous line broadening due to hfi (isotropic)
 % Prisner et al., Time-resolved W-band..., 1995
-lw1 = 10e6;  % Hz
-lw2 = 10e6;  % Hz
+lw1 = 8e6;  % Hz
+lw2 = 8e6;  % Hz
 % The g-factor linewidth
 glw1 = planck/bmagn/B0*lw1;
 glw2 = planck/bmagn/B0*lw2;
 % Spacing between g-values
-dgax = 1e-5;
-
-% pulseBw = 0.0005;  % T
-% gPulseBw = planck*mwFreq/bmagn*(pulseBw/B0/B0);
-% gPulsePosition = planck*mwFreq/bmagn/B0;
-% excitRange = [-1/2, +1/2]*gPulseBw + gPulsePosition;
+dgax = 1e-6;
 
 signalPowder = zeros(numel(xx), 1);
-signal_ = repmat({signalPowder}, nTheta, nPhi);
+signal_ = zeros(nTheta, nPhi, numel(xx));
 xi = zeros(nTheta, nPhi);
 
 tic
@@ -213,13 +225,19 @@ for ith = 30
     end
     waitbar(ith/nTheta, wbarTh, sprintf('%d out of %d', ith, nTheta))
 
-    theta_ = thetas(ith)*pi/180;
+    theta_ = thetas(ith);
     % theta_ = theta_/multTheta;
     % theta_ = theta_ * 10;
-    for iph = 70
-        phi_ = phis(iph)*pi/180;
+    
+    signalPhi_ = zeros(nPhi, numel(xx));
+    ddPhi_ = dd(ith, :);
+    g1Phi_ = g1(ith, :);
+    g2Phi_ = g2(ith, :);
+    disp('hi')
+    parfor iph = 70
+        phi_ = phis(iph);
         
-        % Finite pulse bandwidth: no signal from spins outside excitRange
+%         Finite pulse bandwidth: no signal from spins outside excitRange
         %{
         if g1(ith, iph) < excitRange(1) || ...
             g1(ith, iph) > excitRange(2) || ... 
@@ -231,10 +249,10 @@ for ith = 30
         
         % Consider inhomogenous broadening
         % Create g-axis in order to 'sample' the inh. broadened g-values
-        gaxLarge = creategaxis(g1(ith, iph), g2(ith, iph), glw1, glw2, dgax);
+        gaxLarge = creategaxis(g1Phi_(iph), g2Phi_(iph), glw1, glw2, dgax);
         % Inhomogeneously broadened gaussian distributions
-        gInh1Large = gaussian(gaxLarge, g1(ith, iph), glw1);
-        gInh2Large = gaussian(gaxLarge, g2(ith, iph), glw2);
+        gInh1Large = gaussian(gaxLarge, g1Phi_(iph), glw1);
+        gInh2Large = gaussian(gaxLarge, g2Phi_(iph), glw2);
         % Normalization
         gInh1Large = gInh1Large/sum(gInh1Large);
         gInh2Large = gInh2Large/sum(gInh2Large);
@@ -242,45 +260,29 @@ for ith = 30
         [gInh1, gInh2, gax] = ...
             cutinhdistributions(gInh1Large, gInh2Large, gaxLarge);
         
-        ngax = numel(gax);
-        
-        xi_ = zeros(1, ngax);
-        dd_ = dd(ith, iph);
-        signalInh = zeros(numel(xx), ngax);
         gaxFirstValue = gax(1);
-
         gaussianFactors = gInh1'*gInh2;
-        signalInhFactors = zeros(1, ngax);
         
-%         xi_ = atan( (dd_ + 2*JJ)*1e6 * ...
-%                 planck/(2*pi*(bmagn*0.35)) ./ (gax - gaxFirstValue) );
-%         signalInh = crystalEseem(xi_, xx);
-        parfor ig = 1:ngax
-            % Calculate xi at fixed g2 = gax(1) value
-%             xi_(ig) = atan( (dd_ + 2*JJ)*1e6 * ...
-%                 planck/(2*pi*(bmagn*0.35)) ./ (gax(ig) - gaxFirstValue) );
-%             signalInh(:, ig) = crystalEseem(xi_(ig), xx);
-            
-            % Calculate multiplicity
-%             if ig ~= 1
-%                 iDiag = ig - 1;
-%                 signalInhFactors(ig) = sum(spdiags(gaussianFactors, -iDiag)) + ...
-%                 sum(spdiags(gaussianFactors, iDiag));
-%             else  % ig == 1
-%                 iDiag = ig - 1;  % = 0
-%                 signalInhFactors(ig) = sum(spdiags(gaussianFactors, iDiag));
-%              end
-            iDiag = ig - 1;
-            signalInhFactors(ig) = sum(spdiags(gaussianFactors, -iDiag)) + ...
-                sum(spdiags(gaussianFactors, iDiag));
-        end
+        xi_ = atan( (ddPhi_(iph) + 2*JJ)*1e6 * ...
+                planck/(2*pi*(bmagn*B0)) ./ (gax - gaxFirstValue) );
+        signalInh = crystalEseem(xi_, xx);
         
-        xi(ith, iph) = mean(xi_.*signalInhFactors);
-        signalInh = signalInh.*signalInhFactors;
-        signal_{ith, iph} = sum(signalInh, 2)*sin(theta_);
+        xiComplete_ = [flip(-xi_(2:end)), xi_];
+        xiMultiplicity = sum(spdiags(gaussianFactors), 1);
+        
+        % Calculate inh signal multiplicity using the fact that the signal
+        % is symmetric for xi -> -xi (hence sum the values of the
+        % non-central (meaning g1 - g2 != 0) multiplicity values)
+        inhMultiplicity = xiMultiplicity(end/2 + 0.5:end) + ...
+            xiMultiplicity(end/2 + 0.5:-1:1);
+        inhMultiplicity(1) = inhMultiplicity(1)/2;
+
+        xiPhi_ = mean(xiComplete_.*xiMultiplicity);
+        signalInh = signalInh.*inhMultiplicity;
+        signalPhi_(iph, :) = sum(signalInh, 2)*sin(theta_);
     end
-
-    signalPowder = signalPowder + signal_{ith, iph};
+    toc
+    signal_(ith, :, :) = signalPhi_;
 end
 toc
 
@@ -289,23 +291,35 @@ toc
 delete(wbarTh)
 clear('wbarTh')
 
-clf
-plot(xx, signal_{ith, iph})
-hold on
-plot(xx, signalNoLw_{ith, iph})
+% clf
+% plot(xx, signal_{ith, iph})
+% hold on
+% plot(xx, signalNoLw_{ith, iph})
 
-max(signal_{ith, iph})
+% max(signal_{ith, iph})
 
 %% Test
 
-aa0 = spdiags(gaussianFactors);
-aa1 = sum(aa0, 1);
-aa2 = aa1(end/2 + 0.5:end) + aa1(end/2 + 0.5:-1:1);
+% tempInhMultiplicity = sum(spdiags(gaussianFactors), 1);
+% InhMultiplicity = aa1(end/2 + 0.5:end) + aa1(end/2 + 0.5:-1:1);
 clf
-plot(1:293, aa2, 1:293, signalInhFactors)
+ngax = numel(gax);
+plot(1:2*ngax - 1, xiMultiplicity, 1:ngax, inhMultiplicity)
+hold on
+plot(1:ngax, gInh1, 1:ngax, gInh2)
+plot(1:2*ngax - 1, xiComplete_*1e-4)
 
 %% Test
 
+deltagDistr = linspace(0, 1e-4, 2*ngax - 1);
+xiDistr = atan( (dd(30, 70) + 2*JJ)*1e6 * ...
+        planck/(2*pi*(bmagn*B0)) ./ deltagDistr );
+mean(xiDistr)*180/pi
+atan( (dd(30, 70) + 2*JJ)*1e6 * ...
+        planck/(2*pi*(bmagn*B0)) ./ mean(deltagDistr))*180/pi
+    
+%% Test
+
 ith = 30;
 iph = 400;
 
@@ -345,7 +359,7 @@ colorbar
 
 %%
 
-size(crystalEseem(xi_(1:10)', xx))
+% size(crystalEseem(xi_(1:10)', xx))
 
 %% Some plots
 
@@ -429,11 +443,11 @@ for ith = 1:nTheta
             signal_{ith, iph} = sum(signalSummed, 2);
         end
         % Solid angle normalization
-%         theta_ = thetas(ith)*pi/180;
+%         theta_ = thetas(ith);
 %         signal_{ith, iph} = sum(signalSingleSumg2, 2)*sin(theta_);
 %         signalPowder = signalPowder + signal_{ith, iph};
     end
-    theta_ = thetas(ith)*pi/180;
+    theta_ = thetas(ith);
     signalPowder = signalPowder + signal_{ith, iph}*sin(theta_);
 end
 
@@ -443,25 +457,13 @@ delete(wbarTh)
 
 % signalPowder = zeros(numel(xx), 1);
 % for ith = 1:nTheta
-%     theta_ = thetas(ith)*pi/180;
+%     theta_ = thetas(ith);
 %     for iph = 1:nPhi
 %         signalPowder = signalPowder + signal_{ith, iph}*sin(theta_);
 %     end
 % end
 
-figure(6)
-clf
-plot(x*180/pi, rescaledata(y, 'maxabs'), 'o-')
-hold on
-plot(xx*180/pi, rescaledata(yy, 'maxabs'))
-plot(xx*180/pi, initFit)
-plot(xx*180/pi, initFit2)
-plot(xx*180/pi, expectFit)
-plot(xx*180/pi, rescaledata(signalPowder, 'maxabs'), 'k-')
-legend('Exp. data', 'Fit Zech', 'Fit Gianluca using alpha', ...
-    'Fit Gianluca using xi', 'Sum of ESEEM signals', ...
-    'Location', 'northwest')
-labelaxesfig(gca, 'Flip angle beta', 'ESEEM intensity')
+
 % plot(xx*180/pi, -rescaledata(sin(2*xx), y))
 
 %%
@@ -475,13 +477,13 @@ end
 %% Mean xi
 
 % for ith = 1:nTheta
-%     theta_ = thetas(ith)*pi/180;
+%     theta_ = thetas(ith);
 %     xitest(ith, :) = repmat(cos(theta_)^4, 1, nPhi);
 %     xitest(ith, :) = repmat(1, 1, nPhi);
 % end
 
 for ith = 1:nTheta
-    theta_ = thetas(ith)*pi/180;
+    theta_ = thetas(ith);
     xiMean = xiMean + sum(xi(ith, :))*sin(theta_);
 end
 normTh = pi/2/nTheta;
@@ -522,14 +524,14 @@ plot(aa, gaussian(aa, g1(140, 1), glw1) .* gaussian(aa, g2(140, 1), glw2))
 %%
 
 function xx = creategaxis(g1, g2, glw1, glw2, dg)
-    xxFactor = 1;
+    xxFactor = 10;
     xmin = min([g1 - xxFactor*glw1, g2 - xxFactor*glw2]);
     xmax = max([g1 + xxFactor*glw1, g2 + xxFactor*glw2]);
     xx = xmin:dg:xmax;
 end
 
 function [gInh1, gInh2, gax] = cutinhdistributions(gInh1Large, gInh2Large, gaxLarge)
-    xxFactor = 2;
+    xxFactor = 5;
     multipl = gInh1Large.*gInh2Large;
     [maxMult, idxMax] = max(multipl);
     [~, idxFwhm] = min(abs(multipl - maxMult/2));