diff --git a/ciss_trEPRstickSpectra.asv b/ciss_trEPRstickSpectra.asv
new file mode 100644
index 0000000000000000000000000000000000000000..919b9df5fcf9da041c1bf59a9024fcd4b0dd0118
--- /dev/null
+++ b/ciss_trEPRstickSpectra.asv
@@ -0,0 +1,128 @@
+%
+clearvars
+
+%% System definition
+
+Sys0.S = [1/2 1/2];
+Sys0.g = [2.0030 2.0026 2.0023; ... % P700+
+ 2.0062 2.0051 2.0022]; % A1-
+Sys0.gFrame = [-10 -128 -83; ...
+ 0 0 0]*pi/180;
+Sys0.eeFrame = [0 90 0]*pi/180; % zD direction is along -x of A1-
+
+Sys.g = Sys0.g;
+Sys.gFrame = Sys0.gFrame;
+Sys.dip = mt2mhz(-0.170); % MHz
+Sys.J = mt2mhz(1e-3); % MHz
+Sys.eeFrame = Sys0.eeFrame;
+Sys.nNuc = 3;
+Sys.A = [0, 0, 0, 9, 9.4, 12.8]; % MHz
+Sys.AFrame = [0 0 0 60 -90 0]*pi/180;
+Sys.trlwpp = 0.35;
+
+Sys.rho = zeros(4, 4);
+% Up-down and down-up
+Sys.rho(2, 2) = 1;
+Sys.rho(3, 3) = 1;
+Sys.rho = 1/2*Sys.rho;
+% Singlet
+% Sys.rho(2, 2) = 1;
+% Sys.rho(2, 3) = -1;
+% Sys.rho(3, 2) = -1;
+% Sys.rho(3, 3) = 1;
+% Sys.rho = 1/2*Sys.rho;
+
+Sys.mwFreq = 9.6;
+Sys.x = linspace(340.5, 344, 301); % mT
+Sys.nTheta = 50;
+Sys.nPhi = 20;
+
+%%
+
+test = treprstickspectrum(Sys);
+aa0 = averageoversolidangle(test, Sys.nTheta, Sys.nPhi, 2);
+% clf
+hold on
+box on
+plot(Sys.x, aa0/)
+
+%%
+function signal = treprstickspectrum(Sys)
+
+ fieldAxis = Sys.x;
+ nField = numel(fieldAxis);
+ dip = Sys.dip;
+ JJ = Sys.J;
+ trlwpp = Sys.trlwpp;
+ mwFreq = Sys.mwFreq;
+ nTheta = Sys.nTheta;
+ nPhi = Sys.nPhi;
+ rhoInit = Sys.rho; % This should be in the T+, S, T0, T- basis
+ if isfield(Sys, 'nNuc')
+ nNuc = Sys.nNuc;
+ else
+ nNuc = 0;
+ end
+ nHfiLine = nNuc + 1;
+
+ % Grid
+ [thetas, phis] = createthetaphigrid(nTheta, nPhi);
+ nSolidAngle = nTheta*nPhi;
+
+
+ % Direction of B0
+ nVers = ang2vec(ones(nTheta, 1)*phis, thetas*ones(1, nPhi));
+ size(nVers)
+ % Effective g-values
+ g1Tensor = rotatematrixeuangles(diag(Sys.g(1, :)), Sys.gFrame(1, :));
+ g1 = sqrt( sum( (g1Tensor*nVers).^2, 1));
+ g2Tensor = rotatematrixeuangles(diag(Sys.g(2, :)), Sys.gFrame(2, :));
+ g2 = sqrt( sum( (g2Tensor*nVers).^2, 1));
+ % Dipolar interaction
+ zD = erot(Sys.eeFrame)*[0, 0, 1]';
+ dd = dipinteraction(dip, nVers, zD);
+
+ % Hyperfine
+ if nNuc > 0
+ [APlus, AMinus] = calculateaplusaminusnohfi(...
+ Sys.A, Sys.AFrame, nVers);
+ else
+ APlus = zeros(nSolidAngle, 1);
+ AMinus = zeros(nSolidAngle, 1);
+ end
+
+ signal = zeros(nField, nSolidAngle);
+ for ii = 1:nSolidAngle
+ disp(floor(ii/nSolidAngle*100))
+ gPlus = 1/2*(g1(ii) + g2(ii));
+ gMinus = 1/2*(g1(ii) - g2(ii));
+
+ w0_ = gvalue2freq(fieldAxis, gPlus);
+ deltaw_ = gvalue2freq(fieldAxis, gMinus);
+
+ quantNumNuc = -nNuc/2:1:nNuc/2;
+ if nNuc > 0
+ pascalMatrix = pascal(nHfiLine);
+ % Antidiag
+ pascalFactor = pascalMatrix(nHfiLine:nHfiLine - 1:end - 1);
+ pascalFactor = pascalFactor'/sum(pascalFactor);
+ else
+ pascalFactor = 1;
+ end
+
+ for itrans = 1:4
+ for ihfi = 1:nHfiLine
+ w0__ = w0_ + quantNumNuc(ihfi)*APlus(ihfi);
+ deltaw__ = deltaw_ + quantNumNuc(ihfi)*AMinus(ihfi);
+ wReson = myeigenenergies(w0__, deltaw__, JJ, dd(ii), itrans);
+ intensityReson = intensityresonance(rhoInit, nVers(:, ii), deltaw__, JJ, dd(ii), itrans);
+ signal__ = gaussianresonancebsweep( ...
+ wReson*1e-3, mwFreq, mt2mhz(trlwpp)*1e-3, "lwpp");
+ signal(:, ii) = signal(:, ii) + intensityReson'.*signal__'*pascalFactor(ihfi);
+
+ end
+ end
+ end
+
+end
+
diff --git a/ciss_trEPRstickSpectra.m b/ciss_trEPRstickSpectra.m
new file mode 100755
index 0000000000000000000000000000000000000000..d4de77183b65fc196ead5b9112c849c2e5a10224
--- /dev/null
+++ b/ciss_trEPRstickSpectra.m
@@ -0,0 +1,130 @@
+%
+clearvars
+
+%% System definition
+
+Sys0.S = [1/2 1/2];
+Sys0.g = [2.0030 2.0026 2.0023; ... % P700+
+ 2.0062 2.0051 2.0022]; % A1-
+Sys0.gFrame = [-10 -128 -83; ...
+ 0 0 0]*pi/180;
+Sys0.eeFrame = [0 90 0]*pi/180; % zD direction is along -x of A1-
+
+Sys.g = Sys0.g;
+Sys.gFrame = Sys0.gFrame;
+Sys.dip = mt2mhz(-0.170); % MHz
+Sys.J = mt2mhz(1e-3); % MHz
+Sys.eeFrame = Sys0.eeFrame;
+Sys.nNuc = 3;
+Sys.A = [0, 0, 0, 9, 9.4, 12.8]; % MHz
+Sys.AFrame = [0 0 0 60 -90 0]*pi/180;
+Sys.trlwpp = 0.35;
+
+Sys.rho = zeros(4, 4);
+% Up-down and down-up
+Sys.rho(2, 2) = 1;
+Sys.rho(3, 3) = 1;
+Sys.rho = 1/2*Sys.rho;
+% Singlet
+% Sys.rho(2, 2) = 1;
+% Sys.rho(2, 3) = -1;
+% Sys.rho(3, 2) = -1;
+% Sys.rho(3, 3) = 1;
+% Sys.rho = 1/2*Sys.rho;
+
+Sys.mwFreq = 9.6;
+Sys.x = linspace(340.5, 344, 301); % mT
+Sys.nTheta = 50;
+Sys.nPhi = 20;
+
+%%
+
+test = treprstickspectrum(Sys);
+aa0 = averageoversolidangle(test, Sys.nTheta, Sys.nPhi, 2);
+aa0 = real(aa0);
+figure()
+clf
+hold on
+box on
+plot(Sys.x, aa0/max(abs(aa0)))
+
+%%
+function signal = treprstickspectrum(Sys)
+
+ fieldAxis = Sys.x;
+ nField = numel(fieldAxis);
+ dip = Sys.dip;
+ JJ = Sys.J;
+ trlwpp = Sys.trlwpp;
+ mwFreq = Sys.mwFreq;
+ nTheta = Sys.nTheta;
+ nPhi = Sys.nPhi;
+ rhoInit = Sys.rho; % This should be in the T+, S, T0, T- basis
+ if isfield(Sys, 'nNuc')
+ nNuc = Sys.nNuc;
+ else
+ nNuc = 0;
+ end
+ nHfiLine = nNuc + 1;
+
+ % Grid
+ [thetas, phis] = createthetaphigrid(nTheta, nPhi);
+ nSolidAngle = nTheta*nPhi;
+
+
+ % Direction of B0
+ nVers = ang2vec(ones(nTheta, 1)*phis, thetas*ones(1, nPhi));
+ size(nVers)
+ % Effective g-values
+ g1Tensor = rotatematrixeuangles(diag(Sys.g(1, :)), Sys.gFrame(1, :));
+ g1 = sqrt( sum( (g1Tensor*nVers).^2, 1));
+ g2Tensor = rotatematrixeuangles(diag(Sys.g(2, :)), Sys.gFrame(2, :));
+ g2 = sqrt( sum( (g2Tensor*nVers).^2, 1));
+ % Dipolar interaction
+ zD = erot(Sys.eeFrame)*[0, 0, 1]';
+ dd = dipinteraction(dip, nVers, zD);
+
+ % Hyperfine
+ if nNuc > 0
+ [APlus, AMinus] = calculateaplusaminusnohfi(...
+ Sys.A, Sys.AFrame, nVers);
+ else
+ APlus = zeros(nSolidAngle, 1);
+ AMinus = zeros(nSolidAngle, 1);
+ end
+
+ signal = zeros(nField, nSolidAngle);
+ for ii = 1:nSolidAngle
+ disp(floor(ii/nSolidAngle*100))
+ gPlus = 1/2*(g1(ii) + g2(ii));
+ gMinus = 1/2*(g1(ii) - g2(ii));
+
+ w0_ = gvalue2freq(fieldAxis, gPlus);
+ deltaw_ = gvalue2freq(fieldAxis, gMinus);
+
+ quantNumNuc = -nNuc/2:1:nNuc/2;
+ if nNuc > 0
+ pascalMatrix = pascal(nHfiLine);
+ % Antidiag
+ pascalFactor = pascalMatrix(nHfiLine:nHfiLine - 1:end - 1);
+ pascalFactor = pascalFactor'/sum(pascalFactor);
+ else
+ pascalFactor = 1;
+ end
+
+ for itrans = 1:4
+ for ihfi = 1:nHfiLine
+ w0__ = w0_ + quantNumNuc(ihfi)*APlus(ihfi);
+ deltaw__ = deltaw_ + quantNumNuc(ihfi)*AMinus(ihfi);
+ wReson = myeigenenergies(w0__, deltaw__, JJ, dd(ii), itrans);
+ intensityReson = intensityresonance(rhoInit, nVers(:, ii), deltaw__, JJ, dd(ii), itrans);
+ signal__ = gaussianresonancebsweep( ...
+ wReson*1e-3, mwFreq, mt2mhz(trlwpp)*1e-3, "lwpp");
+ signal(:, ii) = signal(:, ii) + intensityReson'.*signal__'*pascalFactor(ihfi);
+
+ end
+ end
+ end
+
+end
+
diff --git a/reports/20241018_castleMonthlyUpdates/Template.pptx b/reports/20241018_castleMonthlyUpdates/Template.pptx
new file mode 100644
index 0000000000000000000000000000000000000000..bc0c07abb19426d76512f7cf4a91ba25d76b69c7
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/Template.pptx differ
diff --git a/reports/20241018_castleMonthlyUpdates/castle_02_motivationNaamanPSIchirality.jpg b/reports/20241018_castleMonthlyUpdates/castle_02_motivationNaamanPSIchirality.jpg
new file mode 100755
index 0000000000000000000000000000000000000000..8edc1d295e639cc4dbd8a9d4e6598e59eb3cf8f8
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/castle_02_motivationNaamanPSIchirality.jpg differ
diff --git a/reports/20241018_castleMonthlyUpdates/castle_03_02_sigAmpTemperaturePSI.png b/reports/20241018_castleMonthlyUpdates/castle_03_02_sigAmpTemperaturePSI.png
new file mode 100644
index 0000000000000000000000000000000000000000..dbe41311933a1f2e5b517c622d7b9a4b15436b9a
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/castle_03_02_sigAmpTemperaturePSI.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/castle_03_and_06_sigAmpTemperaturePSI.m b/reports/20241018_castleMonthlyUpdates/castle_03_and_06_sigAmpTemperaturePSI.m
new file mode 100644
index 0000000000000000000000000000000000000000..1e4965a9a834d8aeb1a308cad48646ea2262a0af
--- /dev/null
+++ b/reports/20241018_castleMonthlyUpdates/castle_03_and_06_sigAmpTemperaturePSI.m
@@ -0,0 +1,237 @@
+% FIGURE SLIDE 3: castle_03_02_sigAmpPSI.png
+clearvars
+
+% Import
+loadFolder = "/net/storage/gianlum33/projects/zech_psi/data/processed/";
+loadFileName = loadFolder + "ZePSI-E-011";
+
+load(loadFileName)
+xx = x;
+yOld = y;
+yLoad = yOld(1, 1:8);
+
+loadFileName = loadFolder + "ZePSI-E-011-part2";
+load(loadFileName)
+
+yLoad(9:12) = y(4:-1:1); % 230 240 250 260
+
+nMeas = numel(yLoad);
+
+%%
+
+[nt, nB] = size(yLoad{1});
+
+y = yLoad;
+% for ii = 1:nMeas
+% for iB = 1:nB
+% y{ii}(:, iB) = datasmooth(yLoad{ii}(:, iB), 20, 'savgol');
+% end
+% end
+
+% Scrollable traces
+clf
+tiledlayout('flow', 'TileSpacing', 'compact', 'Padding', 'compact')
+for ii = 1:nMeas
+ nexttile
+ h = ScrollableAxes();
+ % ii = 8;
+ if ii < 9
+ plot(h, xx{1}, xx{2}, yLoad{ii});
+ else
+ plot(h, x{1}, x{2}, yLoad{ii});
+ end
+ hold on
+ title(string(60 + ii*20))
+ % plot(h, x{1}, x{2}, y{ii}');
+end
+
+%%
+
+% temperatures = [80:20:260, 265, 260:-10:220];
+temperatures = [80:20:220, 230:10:260];
+% temperatures = 260:-10:220;
+
+noiseRange = [min(x{2}), 338.7; 341.5, max(x{2})];
+[iTimeMax, snr0AtMax, snr1AtMax, ppAmpMax, sigAmpMax, ...
+ noiseLev0AtMax, noiseLev1AtMax, peaksRatio, ...
+ maxR, maxL] = deal(zeros(size(temperatures)));
+for ii = 1:nMeas
+ iTimeMax_ = 0;
+ % ppAmpMax_ = 0;
+ sigAmpMax_ = 0;
+ snrAtMax_ = 0;
+ noise1LevAtMax_ = 0;
+ for it = 1:1200
+ % [snr1, ppAmp, noiseLev1] = getSNR(x{2}, y{ii}(it, :), noiseRange);
+ sigAmp = max(y{ii}(it, :));
+ % if ppAmp > ppAmpMax_
+ if sigAmp > sigAmpMax_
+ iTimeMax_ = it;
+ sigAmpMax_ = sigAmp;
+ % [snr1AtMax_, ppAmpMax_, noiseLev1AtMax_] = getSNR( ...
+ % x{2}, y{ii}(it, :), noiseRange);
+ if ii < 9
+ [~, ~, noiseLev1AtMax_] = getSNR( ...
+ xx{2}, y{ii}(it, :), noiseRange);
+ snr1AtMax_ = sigAmpMax_/noiseLev1AtMax_;
+ [~, ~, noiseLev0AtMax_] = getSNR(...
+ xx{2}, yLoad{ii}(it, :), noiseRange);
+ snr0AtMax_ = sigAmpMax_/noiseLev0AtMax_;
+ else
+ [~, ~, noiseLev1AtMax_] = getSNR( ...
+ x{2}, y{ii}(it, :), noiseRange);
+ snr1AtMax_ = sigAmpMax_/noiseLev1AtMax_;
+ [~, ~, noiseLev0AtMax_] = getSNR(...
+ x{2}, yLoad{ii}(it, :), noiseRange);
+ snr0AtMax_ = sigAmpMax_/noiseLev0AtMax_;
+ end
+ end
+ end
+
+ iTimeMax(ii) = iTimeMax_;
+ snr1AtMax(ii) = snr1AtMax_;
+ snr0AtMax(ii) = snr0AtMax_;
+ % ppAmpMax(ii) = ppAmpMax_;
+ sigAmpMax(ii) = sigAmpMax_;
+ noiseLev1AtMax(ii) = noiseLev1AtMax_;
+ noiseLev0AtMax(ii) = noiseLev0AtMax_;
+ maxR(ii) = max(y{ii}(iTimeMax(ii), :));
+ maxL(ii) = max(y{ii}(iTimeMax(ii), 100)); % iLeftPeak = 100
+ peaksRatio(ii) = maxR/maxL;
+
+end
+% facMult = sigAmpMax(8)/0.52/sigAmpMax(end);
+
+% Adjust for factor between days
+for ii = 1:nMeas
+
+ % if ii >= 12
+ % sigAmpMax(ii) = sigAmpMax(ii)*1.1;
+ % snr1AtMax(ii) = snr1AtMax(ii)*1.1;
+ % snr0AtMax(ii) = snr0AtMax(ii)*1.1;
+ % See ZePSI_011_02_factorBetweenDays.m
+ if ii >= 5 && ii < 9 %&& strcmp(loadFileName, "../data/processed/ZePSI-E-011")
+ % ppAmpMax(ii) = ppAmpMax(ii)/0.52;
+ sigAmpMax(ii) = sigAmpMax(ii)/0.52;
+ snr1AtMax(ii) = snr1AtMax(ii)/0.52;
+ snr0AtMax(ii) = snr0AtMax(ii)/0.52;
+ end
+ if ii >= 9
+ % ppAmpMax(ii) = ppAmpMax(ii)/2.56;
+ sigAmpMax(ii) = sigAmpMax(ii)/1.1;
+ snr1AtMax(ii) = snr1AtMax(ii)/1.1;
+ snr0AtMax(ii) = snr0AtMax(ii)/1.1;
+ end
+end
+
+figure(4)
+xTicks = temperatures;
+% xTicks = flip(temperatures);
+
+clf
+tiledlayout('flow', 'TileSpacing', 'compact', 'Padding', 'compact')
+nexttile(1, [2, 4])
+cmap = viridis(nMeas);
+for ii = 1:nMeas
+ % yplot = y{ii}(4000, :);
+ if ii < 9
+ xplot = xx;
+ yplot = rescaledata(y{ii}(400, :), 'maxabs');
+ else
+ xplot = x;
+ yplot = rescaledata(y{ii}(200, :), 'maxabs');
+ end
+ plot(xplot{2}, yplot, ...
+ 'DisplayName', string(temperatures(ii)) + " K", ...
+ 'Color', cmap(ii, :))
+ hold on
+end
+legend()
+xlim(setaxlim(x{2}, 1))
+xline(x{2}(100), 'DisplayName', 'Position of shoulder')
+
+nexttile
+plot(temperatures, snr1AtMax, 'o-') % Filtered ppAmp / filtered noise
+yyaxis right
+% plot(temperatures, snr0AtMax, 'o-') % Filtered ppAmp / non-filtered noise
+xlim(setaxlim(temperatures, 1.05))
+% xticks(xTicks)
+title('SNR at max')
+legend('Smooth', 'No filter', 'Location', 'southwest')
+%
+nexttile
+plot(temperatures, noiseLev1AtMax, 'o-') % Filtered noise
+yyaxis right
+% plot(temperatures, noiseLev0AtMax, 'o-') % Non-filtered noise
+xlim(setaxlim(temperatures, 1.05))
+% xticks(xTicks)
+title('noiseLev at max')
+%
+nexttile
+errorbar(temperatures, sigAmpMax, 2*noiseLev0AtMax, 'o-')
+xline(140)
+xline(220)
+xlim(setaxlim(temperatures, 1.05))
+% xticks(xTicks)
+% title('ppAmp max')
+title('max of signal')
+%
+nexttile
+% errbar from non-filtered data. Also possible: hypot(1./maxR, 1./maxL)
+errbar = peaksRatio.*noiseLev0AtMax.*(1./maxR + 1./maxL);
+errorbar(temperatures, peaksRatio, errbar, 'o-')
+hold on
+xlim(setaxlim(temperatures, 1.05))
+ylim(setaxlim(peaksRatio, 1.5))
+xticks(xTicks)
+title('peaks ratio at max')
+
+%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% FIGURE castle_03_02_sigAmpTemperaturePSI
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% figure()
+clf
+errorbar(temperatures, sigAmpMax/max(sigAmpMax), 5*noiseLev0AtMax/max(sigAmpMax), 'o-k')
+ax0 = gca;
+
+xlim(setaxlim(temperatures, 1.1))
+ylim(setaxlim(sigAmpMax/max(sigAmpMax), 1.125))
+xticks(100:40:260)
+yticks(0:0.5:1)
+
+set(gca, 'Xdir', 'reverse')
+labelaxesfig(gca, 'Temperature / K', 'trEPR amplitude norm. / a.u.')
+
+saveFolder = "/net/storage/gianlum33/projects/oop_ciss_calculations/reports/20241018_castleMonthlyUpdates";
+
+savePath = fullfile(saveFolder, "castle_03_02_sigAmpTemperaturePSI.png");
+% exportgraphics(gcf, savePath, 'Resolution', 600)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% FIGURE castle_06_01_zoomSigAmpTemperaturePSI
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Create a new pair of axes inside current figure
+ax1 = axes('Position', [.475 .25 .4 .4]);
+box on
+
+iPlot = 6:12;
+xplot = temperatures(iPlot);
+yplot = sigAmpMax(iPlot)/max(sigAmpMax);
+errBar = 5*noiseLev0AtMax(iPlot)/max(sigAmpMax);
+
+errorbar(xplot, yplot, errBar, 'o-k')
+
+xticks(ax1, 200:10:260)
+yticks(ax1, 0:0.5:1)
+xlim(setaxlim([200, 260], 1.075))
+ylim([-0.05 1])
+% coordPos = [0.15, 0.175, 0.275, 0.6];
+% annotation('rectangle', coordPos)
+% annotation("arrow", [coordPos(1) + coordPos(3), 0.6 )
+set(gca, 'Xdir', 'reverse')
+
+savePath = fullfile(saveFolder, "castle_06_01_zoomSigAmpTemperaturePSI.png");
+% exportgraphics(gcf, savePath, 'Resolution', 600)
\ No newline at end of file
diff --git a/reports/20241018_castleMonthlyUpdates/castle_04_01_simulExample.png b/reports/20241018_castleMonthlyUpdates/castle_04_01_simulExample.png
new file mode 100644
index 0000000000000000000000000000000000000000..9436a6730775d2300716e27c919711ed34ee6c29
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/castle_04_01_simulExample.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/castle_04_02_simulExample.png b/reports/20241018_castleMonthlyUpdates/castle_04_02_simulExample.png
new file mode 100644
index 0000000000000000000000000000000000000000..661b4c2f1ec707c5ae4766ce317752f5c90df53c
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/castle_04_02_simulExample.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/castle_04_03_simulExample.png b/reports/20241018_castleMonthlyUpdates/castle_04_03_simulExample.png
new file mode 100644
index 0000000000000000000000000000000000000000..14befbfd6d9772d76b140cccf0fb121e3a39c7e8
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/castle_04_03_simulExample.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/castle_04_04_simulExample.png b/reports/20241018_castleMonthlyUpdates/castle_04_04_simulExample.png
new file mode 100644
index 0000000000000000000000000000000000000000..2f342ee15fd0bcb64d15f7460c9ad460bd8eabbf
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/castle_04_04_simulExample.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/castle_04_05_cissPowder.png b/reports/20241018_castleMonthlyUpdates/castle_04_05_cissPowder.png
new file mode 100755
index 0000000000000000000000000000000000000000..4e9ad0caa6b173b5cf2d2a714513182834153c26
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/castle_04_05_cissPowder.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/castle_04_simulExample.m b/reports/20241018_castleMonthlyUpdates/castle_04_simulExample.m
new file mode 100644
index 0000000000000000000000000000000000000000..b3b20e5c15f7a6fddc264013227108260034801e
--- /dev/null
+++ b/reports/20241018_castleMonthlyUpdates/castle_04_simulExample.m
@@ -0,0 +1,414 @@
+%
+saveFolder = "/net/storage/gianlum33/projects/oop_ciss_calculations/reports/20241018_castleMonthlyUpdates";
+load('plotColors.mat')
+
+FLAG_CALCULATE_INH_BROAD = 0;
+FLAG_SAVE_FIG = 1;
+
+g1 = 2.006;
+g2 = 2.002;
+bfieldg1 = mhz2mt(9.7, g1)*1e3;
+bfieldg2 = mhz2mt(9.7, g2)*1e3;
+
+d0 = mhz2mt(-5); % mT
+J = 0;
+
+nPoint = 1500;
+bWidth = 2.5; % mT
+
+gInhLwSmall = 0.1;
+gInhLwLarge = 1;
+trlwppSmall = gInhLwSmall - 0.01;
+trlwppLarge = gInhLwLarge - 0.2;
+
+gInhLw = gInhLwSmall;
+% gInhLw = gInhLwLarge;
+trlwpp = trlwppSmall;
+% trlwpp = trlwppLarge;
+
+bCenter = 1/2*(bfieldg2 + bfieldg1);
+bDiff = 1/2*(bfieldg2 - bfieldg1);
+
+Omega = hypot(bDiff, (J + d0/2))'; % mT
+bMin = bCenter - bWidth/2;
+bMax = bCenter + bWidth/2;
+xx = linspace(bMin, bMax, nPoint);
+
+% Find bres1 and bres2 in xx
+[~, ires1] = min(abs(xx - bfieldg1));
+[~, ires2] = min(abs(xx - bfieldg2));
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% FIGURE 1
+% castle_04_01_simulExample.png
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+fig = figure(1);
+clf
+plot([xx(ires1), xx(ires2)], [0, 0], 'xk', 'HandleVisibility', 'off')
+box on
+yline(0, 'HandleVisibility', 'off')
+xline([bfieldg1, bfieldg2], ':', 'HandleVisibility', 'off')
+
+ann1 = "g_{eff} = " + string(g1);
+xann1 = 344.8;
+xann2 = 346.25;
+yann = -0.13;
+text(xann1, yann, ann1, 'FontSize', 14)
+ann2 = "g_{eff} = " + string(g2);
+text(xann2, yann, ann2, 'FontSize', 14)
+
+xlim(setaxlim(xx))
+xticks(343:1:350)
+yticks(-1:0.5:1)
+labelaxesfig(gca, 'Magnetic field / mT', 'trEPR intensity / a.u.')
+
+if gInhLw == gInhLwSmall && FLAG_SAVE_FIG
+ savePath = fullfile(saveFolder, "castle_04_01_simulExample.png");
+ exportgraphics(fig, savePath, 'Resolution', 600)
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% FIGURE 2
+% castle_04_02_simulExample
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+bres1 = bCenter - (J - d0) - Omega;
+bres2 = bCenter + (J - d0) - Omega;
+bres3 = bCenter - (J - d0) + Omega;
+bres4 = bCenter + (J - d0) + Omega;
+
+fig = figure(2);
+clf
+plot([xx(ires1), xx(ires2)], [0, 0], 'xk', 'HandleVisibility', 'off')
+box on
+yline(0, 'HandleVisibility', 'off')
+xline([bres1, bres2, bres3, bres4], '--', 'HandleVisibility', 'off')
+
+labelaxesfig(gca, 'Magnetic field / mT', 'trEPR intensity / a.u.')
+ax = gca;
+axCenter = ax.Position(1) + ax.Position(3)/2;
+xOmega = [-1, 1]*Omega/bWidth*ax.Position(3) + axCenter;
+xDiff1 = [-1, 1]*d0/bWidth*ax.Position(3) + xOmega(1);
+xDiff2 = [-1, 1]*d0/bWidth*ax.Position(3) + xOmega(2);
+annotation('doublearrow', xOmega, [.45, .45])
+annotation('doublearrow', xDiff1, [.7, .7])
+annotation('doublearrow', xDiff2, [.7, .7])
+
+xlim(setaxlim(xx))
+xticks(343:1:350)
+yticks(-1:0.5:1)
+
+if gInhLw == gInhLwSmall && FLAG_SAVE_FIG
+ savePath = fullfile(saveFolder, "castle_04_02_simulExample.png");
+ exportgraphics(fig, savePath, 'Resolution', 600)
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% FIGURE 3
+% castle_04_03_simulExample
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+scale = 0.75;
+yspc1 = scale*gaussian(xx, bres1, trlwpp)/max(gaussian(xx, bres1, trlwpp));
+yspc2 = scale*gaussian(xx, bres2, trlwpp)/max(gaussian(xx, bres2, trlwpp));
+yspc3 = scale*gaussian(xx, bres3, trlwpp)/max(gaussian(xx, bres3, trlwpp));
+yspc4 = scale*gaussian(xx, bres4, trlwpp)/max(gaussian(xx, bres4, trlwpp));
+yspc = (yspc1 - yspc2 + yspc3 - yspc4);
+fig = figure(3);
+clf
+plot([xx(ires1), xx(ires2)], [0, 0], 'xk', 'HandleVisibility', 'off')
+box on
+hold on
+% plot(0, 0, 'HandleVisibility', 'off') % Chaotic way to skip one color
+% plot(xx, yspc1, 'DisplayName', 'Reson. 1')
+% plot(xx, -yspc2, 'DisplayName', 'Reson. 2')
+% plot(xx, yspc3, 'DisplayName', 'Reson. 3')
+% plot(xx, -yspc4, 'DisplayName', 'Reson. 4')
+plot(xx, yspc, 'Color', plotColors(2, :), 'DisplayName', 'No CISS')
+% yline(0, 'HandleVisibility', 'off')
+xline([bres1, bres2, bres3, bres4], '--', 'HandleVisibility', 'off')
+
+labelaxesfig(gca, 'Magnetic field / mT', 'trEPR intensity / a.u.')
+legend()
+
+xlim(setaxlim(xx))
+ylim([-1, 1])
+xticks(343:1:350)
+yticks(-1:0.5:1)
+
+if gInhLw == gInhLwSmall && FLAG_SAVE_FIG
+ savePath = fullfile(saveFolder, "castle_04_03_simulExample.png");
+ exportgraphics(fig, savePath, 'Resolution', 600)
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% FIGURE 4
+% castle_04_04_simulExample
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+scale = 0.75;
+yspc1 = 5/4*scale*gaussian(xx, bres1, trlwpp)/max(gaussian(xx, bres1, trlwpp));
+yspc2 = 5/4*scale*gaussian(xx, bres2, trlwpp)/max(gaussian(xx, bres2, trlwpp));
+yspc3 = 3/4*scale*gaussian(xx, bres3, trlwpp)/max(gaussian(xx, bres3, trlwpp));
+yspc4 = 3/4*scale*gaussian(xx, bres4, trlwpp)/max(gaussian(xx, bres4, trlwpp));
+yspcCISS = (yspc1 - yspc2 + yspc3 - yspc4);
+fig = figure(4);
+clf
+plot([xx(ires1), xx(ires2)], [0, 0], 'xk', 'HandleVisibility', 'off')
+box on
+hold on
+% plot(0, 0, 'HandleVisibility', 'off') % Chaotic way to skip one color
+% plot(xx, yspc1, 'DisplayName', 'Reson. 1')
+% plot(xx, -yspc2, 'DisplayName', 'Reson. 2')
+% plot(xx, yspc3, 'DisplayName', 'Reson. 3')
+% plot(xx, -yspc4, 'DisplayName', 'Reson. 4')
+plot(xx, yspc, 'Color', plotColors(2, :), 'DisplayName', 'No CISS')
+hold on
+plot(xx, yspcCISS, 'Color', plotColors(6, :), 'DisplayName', 'CISS')
+
+% yline(0, 'HandleVisibility', 'off')
+xline([bres1, bres2, bres3, bres4], '--', 'HandleVisibility', 'off')
+
+labelaxesfig(gca, 'Magnetic field / mT', 'trEPR intensity / a.u.')
+legend()
+
+xlim(setaxlim(xx))
+ylim([-1, 1])
+xticks(343:1:350)
+yticks(-1:0.5:1)
+
+if gInhLw == gInhLwSmall && FLAG_SAVE_FIG
+ savePath = fullfile(saveFolder, "castle_04_04_simulExample.png");
+ exportgraphics(fig, savePath, 'Resolution', 600)
+end
+
+%%
+%{
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% FIGURE 4
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ scale = 0.75;
+
+fig = figure(4);
+clf
+plot([xx(ires1), xx(ires2)], [0, 0], 'xk', 'HandleVisibility', 'off')
+box on
+hold on
+plot(xx, yspc, 'k', 'DisplayName', 'Method 1')
+% yline(0, 'HandleVisibility', 'off')
+
+labelaxesfig(gca, 'Magnetic field / mT', 'trEPR intensity / a.u.')
+legend()
+
+xlim(setaxlim(xx))
+ylim([-1, 1])
+xticks(343:1:350)
+yticks(-1:0.5:1)
+
+if gInhLw == gInhLwSmall && FLAG_SAVE_FIG
+ savePath = fullfile(saveFolder, "trEPRsimul_example_04.png");
+ exportgraphics(fig, savePath, 'Resolution', 600)
+end
+%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% FIGURE 5 AND FIGURE 9
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+bgInh1 = gaussian(xx, bfieldg1, gInhLw);
+bgInh1 = bgInh1/max(bgInh1)*0.5;
+bgInh2 = gaussian(xx, bfieldg2, gInhLw);
+bgInh2 = bgInh2/max(bgInh2)*0.5;
+
+fig = figure(5);
+clf
+plot([xx(ires1), xx(ires2)], [0, 0], 'xk', 'HandleVisibility', 'off')
+box on
+hold on
+plot(xx, bgInh1, '--', 'Color', plotColors(2, :), 'DisplayName', 'g-value distr.')
+plot(xx, bgInh2, '--', 'Color', plotColors(2, :), 'HandleVisibility', 'off')
+yline(0, 'HandleVisibility', 'off')
+xline([bfieldg1, bfieldg2], ':', 'HandleVisibility', 'off')
+
+text(xann1, yann, ann1, 'FontSize', 14)
+text(xann2, yann, ann2, 'FontSize', 14)
+
+xlim(setaxlim(xx))
+ylim([-1, 1])
+xticks(343:1:350)
+yticks(-1:0.5:1)
+labelaxesfig(gca, 'Magnetic field / mT', 'trEPR intensity / a.u.')
+legend()
+
+if FLAG_SAVE_FIG
+ if gInhLw == gInhLwSmall
+ savePath = fullfile(saveFolder, "trEPRsimul_newVsOld_05.png");
+ % exportgraphics(fig, savePath, 'Resolution', 600)
+ elseif gInhLw == gInhLwLarge
+ savePath = fullfile(saveFolder, "trEPRsimul_newVsOld_09.png");
+ % exportgraphics(fig, savePath, 'Resolution', 600)
+ end
+end
+
+%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% FIGURE 6
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+iCut = bgInh1 > 0.045;
+bgInh1o = bgInh1(iCut);
+xx1o = xx(iCut);
+iCut = 1:2:numel(xx1o);
+bgInh1o = bgInh1o(iCut);
+xx1o = xx1o(iCut);
+iCut = 1:5:numel(bgInh1o);
+bgInh1oPlot = bgInh1o(iCut);
+xx1oPlot = xx1o(iCut);
+
+iCut = bgInh2 > 0.045;
+bgInh2o = bgInh2(iCut);
+xx2o = xx(iCut);
+iCut = 1:2:numel(xx2o);
+bgInh2o = bgInh2o(iCut);
+xx2o = xx2o(iCut);
+iCut = 1:5:numel(bgInh2o);
+bgInh2oPlot = bgInh2o(iCut);
+xx2oPlot = xx2o(iCut);
+
+fig = figure(6);
+clf
+plot([xx(ires1), xx(ires2)], [0, 0], 'xk', 'HandleVisibility', 'off')
+box on
+hold on
+plot(xx, bgInh1, '--', 'Color', plotColors(2, :), 'DisplayName', 'g-value distr.')
+plot(xx, bgInh2, '--', 'Color', plotColors(2, :), 'HandleVisibility', 'off')
+plot(xx1oPlot, bgInh1oPlot, 'o', 'Color', plotColors(2, :), 'DisplayName', 'g-value sampling')
+plot(xx2oPlot, bgInh2oPlot, 'o', 'Color', plotColors(2, :), 'HandleVisibility', 'off')
+yline(0, 'HandleVisibility', 'off')
+xline([bfieldg1, bfieldg2], ':', 'HandleVisibility', 'off')
+
+text(xann1, yann, ann1, 'FontSize', 14)
+text(xann2, yann, ann2, 'FontSize', 14)
+
+xlim(setaxlim(xx))
+ylim([-1, 1])
+xticks(343:1:350)
+yticks(-1:0.5:1)
+labelaxesfig(gca, 'Magnetic field / mT', 'trEPR intensity / a.u.')
+legend()
+
+if gInhLw == gInhLwSmall && FLAG_SAVE_FIG
+ % savePath = fullfile(saveFolder, "trEPRsimul_newVsOld_06.png");
+ % exportgraphics(fig, savePath, 'Resolution', 600)
+end
+
+%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% FIGURE 7
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if gInhLw == gInhLwSmall
+ smallLw = gInhLw/5;
+else
+ smallLw = gInhLw/10;
+end
+gaussianFac = bgInh1o'*bgInh2o;
+if FLAG_CALCULATE_INH_BROAD
+ yspcInh = zeros(size(yspc));
+ for ii = 1:1
+ for jj = 1:1
+ bCenter_ = 1/2*(xx1o(ii) + xx2o(jj));
+ bDiff_ = 1/2*(xx1o(ii) - xx2o(jj));
+
+ Omega_ = hypot(bDiff_, (J + d0/2))'; % mT
+
+ bres1_ = bCenter_ - (J - d0) - Omega_;
+ bres2_ = bCenter_ + (J - d0) - Omega_;
+ bres3_ = bCenter_ - (J - d0) + Omega_;
+ bres4_ = bCenter_ + (J - d0) + Omega_;
+
+ yspc_ = gaussian(xx, bres1_, smallLw)/max(gaussian(xx, bres1_, smallLw));
+ yspc_ = yspc_ - gaussian(xx, bres2_, smallLw)/max(gaussian(xx, bres2_, smallLw));
+ yspc_ = yspc_ + gaussian(xx, bres3_, smallLw)/max(gaussian(xx, bres3_, smallLw));
+ yspc_ = yspc_ - gaussian(xx, bres4_, smallLw)/max(gaussian(xx, bres4_, smallLw));
+ yspcInh = yspc_*gaussianFac(ii, jj);
+
+ end
+ end
+end
+
+fig = figure(7);
+clf
+plot([xx(ires1), xx(ires2)], [0, 0], 'xk', 'HandleVisibility', 'off')
+box on
+hold on
+plot(xx, yspcInh/max(yspcInh)*scale, 'Color', plotColors(2, :))
+xline([bres1_, bres2_, bres3_, bres4_], '--', 'HandleVisibility', 'off')
+
+xlim(setaxlim(xx))
+ylim([-1, 1])
+xticks(343:1:350)
+yticks(-1:0.5:1)
+labelaxesfig(gca, 'Magnetic field / mT', 'trEPR intensity / a.u.')
+
+if gInhLw == gInhLwSmall && FLAG_SAVE_FIG
+ % savePath = fullfile(saveFolder, "trEPRsimul_newVsOld_07.png");
+ % exportgraphics(fig, savePath, 'Resolution', 600)
+end
+
+%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% FIGURE 8 AND 10
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+fig = figure(8);
+
+clf
+gaussianFac = bgInh1o'*bgInh2o;
+if FLAG_CALCULATE_INH_BROAD
+ yspcInh = zeros(size(yspc));
+ for ii = 1:numel(xx1o)
+ for jj = 1:numel(xx2o)
+ bCenter_ = 1/2*(xx1o(ii) + xx2o(jj));
+ bDiff_ = 1/2*(xx1o(ii) - xx2o(jj));
+
+ Omega_ = hypot(bDiff_, (J + d0/2))'; % mT
+
+ bres1_ = bCenter_ - (J - d0) - Omega_;
+ bres2_ = bCenter_ + (J - d0) - Omega_;
+ bres3_ = bCenter_ - (J - d0) + Omega_;
+ bres4_ = bCenter_ + (J - d0) + Omega_;
+
+ yspc_ = gaussian(xx, bres1_, smallLw)/max(gaussian(xx, bres1_, smallLw));
+ yspc_ = yspc_ - gaussian(xx, bres2_, smallLw)/max(gaussian(xx, bres2_, smallLw));
+ yspc_ = yspc_ + gaussian(xx, bres3_, smallLw)/max(gaussian(xx, bres3_, smallLw));
+ yspc_ = yspc_ - gaussian(xx, bres4_, smallLw)/max(gaussian(xx, bres4_, smallLw));
+ yspcInh = yspcInh+yspc_*gaussianFac(ii, jj);
+ end
+ end
+end
+
+plot([xx(ires1), xx(ires2)], [0, 0], 'xk', 'HandleVisibility', 'off')
+box on
+hold on
+plot(xx, yspcInh/max(yspcInh)*scale, 'k', 'DisplayName', 'Method 1')
+plot(xx, yspc/max(yspc)*scale, 'Color', plotColors(2, :), 'DisplayName', 'Method 2')
+
+xlim(setaxlim(xx))
+ylim([-1, 1])
+xticks(343:1:350)
+yticks(-1:0.5:1)
+labelaxesfig(gca, 'Magnetic field / mT', 'trEPR intensity / a.u.')
+legend()
+
+if FLAG_SAVE_FIG
+ if gInhLw == gInhLwSmall
+ % savePath = fullfile(saveFolder, "trEPRsimul_newVsOld_08.png");
+ % exportgraphics(fig, savePath, 'Resolution', 600)
+ % elseif gInhLw == gInhLwLarge
+ % savePath = fullfile(saveFolder, "trEPRsimul_newVsOld_10.png");
+ % exportgraphics(fig, savePath, 'Resolution', 600)
+ end
+end
+%}
+
diff --git a/reports/20241018_castleMonthlyUpdates/castle_05_01_oneSpc.png b/reports/20241018_castleMonthlyUpdates/castle_05_01_oneSpc.png
new file mode 100644
index 0000000000000000000000000000000000000000..19cba43deefced483f82a04afb1f77aaa3123a0f
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/castle_05_01_oneSpc.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/castle_05_02_twoSpc.png b/reports/20241018_castleMonthlyUpdates/castle_05_02_twoSpc.png
new file mode 100644
index 0000000000000000000000000000000000000000..642a16f79d643d01c55ea6638ed12f627468a01f
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/castle_05_02_twoSpc.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/castle_05_03_threeSpc.png b/reports/20241018_castleMonthlyUpdates/castle_05_03_threeSpc.png
new file mode 100644
index 0000000000000000000000000000000000000000..8ab3e15b0fa2bb201d56ac6391130bbc3d22bcd0
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/castle_05_03_threeSpc.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/castle_05_04_fourSpc.png b/reports/20241018_castleMonthlyUpdates/castle_05_04_fourSpc.png
new file mode 100644
index 0000000000000000000000000000000000000000..645d3153ffb85fd24fafd8f3a9d73ee62e998fc6
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/castle_05_04_fourSpc.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/castle_05_specShapeTimeAndTemperature.m b/reports/20241018_castleMonthlyUpdates/castle_05_specShapeTimeAndTemperature.m
new file mode 100644
index 0000000000000000000000000000000000000000..8b9f17b3193eaef01ad7ad294741d466eee4cda5
--- /dev/null
+++ b/reports/20241018_castleMonthlyUpdates/castle_05_specShapeTimeAndTemperature.m
@@ -0,0 +1,126 @@
+% FIGURES in slide 5: Transient EPR vs temperature vs time
+clearvars
+
+% Import
+loadFolder = "/net/storage/gianlum33/projects/zech_psi/data/processed/";
+loadFileName = loadFolder + "ZePSI-E-011";
+
+load(loadFileName)
+xx = x;
+yOld = y;
+yLoad = yOld(1, 1:8);
+
+loadFileName = loadFolder + "ZePSI-E-011-part2";
+load(loadFileName)
+
+yLoad(9:12) = y(4:-1:1); % 230 240 250 260
+
+nMeas = numel(yLoad);
+
+%%
+
+[nt, nB] = size(yLoad{1});
+
+y = yLoad;
+% for ii = 1:nMeas
+% for iB = 1:nB
+% y{ii}(:, iB) = datasmooth(yLoad{ii}(:, iB), 20, 'savgol');
+% end
+% end
+
+% Scrollable traces
+clf
+tiledlayout('flow', 'TileSpacing', 'compact', 'Padding', 'compact')
+for ii = 1:nMeas
+ nexttile
+ h = ScrollableAxes();
+ % ii = 8;
+ if ii < 9
+ plot(h, xx{1}, xx{2}, yLoad{ii});
+ else
+ plot(h, x{1}, x{2}, yLoad{ii});
+ end
+ hold on
+ title(string(60 + ii*20))
+ % plot(h, x{1}, x{2}, y{ii}');
+end
+
+%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% FIGURES castle_05_0x_{oneSpec|twoSpec|threeSpec|fourSpec}.png
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+temperOld = [80:20:260, 265];
+iTemps = [4, 8, 9, 10];
+iTimes = [500, 1500, 2000, 3000];
+cmap = inferno(numel(iTemps) + 1);
+fig = figure(3);
+fig.Position = [680 578 1.75*560 3/2*420];
+clf
+tL = tiledlayout('flow', 'TileSpacing', 'none', 'Padding', 'compact');
+for ii = 1:numel(iTemps)
+ iTemp_ = iTemps(ii);
+ for it = 1:numel(iTimes)
+ ax(it) = nexttile(it);
+ it_ = iTimes(it);
+
+ if ii ~= 4 || it == 1
+ yplot = rescaledata(yOld{iTemp_}(it_, :), 'maxabs');
+ xplot = xx{2};
+ % yplot = y{ii}(it, :);
+ else
+ yplot = datasmooth(yLoad{end}(floor(it_/2), :), 4, 'savgol');
+ yplot = rescaledata(yplot, 'maxabs');
+ xplot = x{2};
+ end
+ plot(xplot, yplot, ...
+ 'DisplayName', string(temperOld(iTemp_)) + " K", ...
+ 'Color', cmap(ii, :))
+
+ xlim(setaxlim(xx{2}, 1))
+ ylim(setaxlim(...
+ [rescaledata(yOld{iTemps(1)}(it_, :), 'maxabs'), ...
+ rescaledata(yOld{iTemps(2)}(it_, :), 'maxabs'), ...
+ rescaledata(yOld{iTemps(3)}(it_, :), 'maxabs'), ...
+ rescaledata(yOld{iTemps(4)}(it_, :), 'maxabs')] , 1.05))
+
+ xticks(339:342)
+ yticks(-1:1)
+ ylim(setaxlim([-1, 1], 1.05))
+ switch it
+ case 1
+ xticklabels('')
+ case 2
+ xticklabels('')
+ yticklabels('')
+ case 3
+ %
+ case 4
+ yticklabels('')
+ end
+
+ hold on
+ if ii == 1
+ text(341.5, 0.8, "t = " + string(round(xx{1}(it_)/1e3)) + " us", ...
+ 'FontSize', 14)
+ end
+ end
+
+
+ % rescaledata(yLoad{end}(floor(it_/2), :), 'maxabs')
+ % Get space above the plots
+ lgd1 = legend(ax(2), 'Location', 'northoutside', 'NumColumns', ii);
+ % Actual legend
+ % lgd1.Position(1) = 0.1;
+ % lgd = legend(ax(2), 'Location', 'northoutside', ...
+ % 'Position', lgd1.Position, 'NumColumns', ii);
+
+ labelaxesfig(tL, 'Magnetic field / mT', 'trEPR signal norm. / a.u.')
+
+ saveFolder = "/net/storage/gianlum33/projects/oop_ciss_calculations/reports/20241018_castleMonthlyUpdates";
+ saveStrEnd = {"oneSpc", "twoSpc", "threeSpc", "fourSpc"};
+ savePath = fullfile(saveFolder, ...
+ "castle_05_0" + string(ii) + "_" + saveStrEnd{ii} + ".png");
+ exportgraphics(gcf, savePath, 'Resolution', 600)
+end
+%%
+fig.Position(4) = fig.Position(4) + 100;
diff --git a/reports/20241018_castleMonthlyUpdates/castle_06_01_zoomSigAmpTemperaturePSI.png b/reports/20241018_castleMonthlyUpdates/castle_06_01_zoomSigAmpTemperaturePSI.png
new file mode 100644
index 0000000000000000000000000000000000000000..443c64b81f52662a3c87d51660c35cb80bd5d1bb
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/castle_06_01_zoomSigAmpTemperaturePSI.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/castle_extra_01_naaman.jpg b/reports/20241018_castleMonthlyUpdates/castle_extra_01_naaman.jpg
new file mode 100755
index 0000000000000000000000000000000000000000..c937e6e37eb7b5e342f4e66bafc0be75d1c3c621
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/castle_extra_01_naaman.jpg differ
diff --git a/reports/20241018_castleMonthlyUpdates/castle_extra_02_naaman.jpg b/reports/20241018_castleMonthlyUpdates/castle_extra_02_naaman.jpg
new file mode 100755
index 0000000000000000000000000000000000000000..2b822e38046278244657719b9272cab9cf278d6a
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/castle_extra_02_naaman.jpg differ
diff --git a/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_example_01.png b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_example_01.png
new file mode 100644
index 0000000000000000000000000000000000000000..f201451bf33eabf03c232144366bf56abc2562db
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_example_01.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_example_02.png b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_example_02.png
new file mode 100644
index 0000000000000000000000000000000000000000..af9e01a22eef49bf6f369747f45a9addbda03ebf
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_example_02.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_example_03.png b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_example_03.png
new file mode 100644
index 0000000000000000000000000000000000000000..c571e3ad157e675d65b0042f5c3ad372a611ab16
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_example_03.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_example_04.png b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_example_04.png
new file mode 100644
index 0000000000000000000000000000000000000000..dc47041bb10e887db7f6b8863cf0183d60c9325f
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_example_04.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_01.png b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_01.png
new file mode 100644
index 0000000000000000000000000000000000000000..64c09a91c3d1b9e5e91e571710357c0c8efcf2c6
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_01.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_02.png b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_02.png
new file mode 100644
index 0000000000000000000000000000000000000000..4e31ae04837e2e1e830a9bbbe87541aae171bd74
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_02.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_03.png b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_03.png
new file mode 100644
index 0000000000000000000000000000000000000000..f12fcfab4467dc66f96fc4a8ff0a9f146c862c4b
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_03.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_04.png b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_04.png
new file mode 100644
index 0000000000000000000000000000000000000000..b417e0a1f668682ab779bb46ff4fb44c00923b0d
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_04.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_05.png b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_05.png
new file mode 100644
index 0000000000000000000000000000000000000000..b6065132b3ee820e16ea61ca4e5cf0252554c29b
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_05.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_06.png b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_06.png
new file mode 100644
index 0000000000000000000000000000000000000000..3ec20cd41745e4a08270fa45702b2dcfe456bdec
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_06.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_07.png b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_07.png
new file mode 100644
index 0000000000000000000000000000000000000000..6ecee10b92b539f1932d0ebff521244ae3c16ccc
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_07.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_08.png b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_08.png
new file mode 100644
index 0000000000000000000000000000000000000000..ce5fad131c2104910f8c70085c0d9c4ce9acef91
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_08.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_09.png b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_09.png
new file mode 100644
index 0000000000000000000000000000000000000000..bd0da3efd0356fc620fbb21191f9f8ee8f11179d
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_09.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_10.png b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_10.png
new file mode 100644
index 0000000000000000000000000000000000000000..b6cb4c0b60edc4b8e6644841531fb6b9b7480566
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/not_used/trEPRsimul_newVsOld_10.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/trEPRsimul_comparisonOnZechPSI.m b/reports/20241018_castleMonthlyUpdates/trEPRsimul_comparisonOnZechPSI.m
new file mode 100644
index 0000000000000000000000000000000000000000..06ca75af31fc92167c784648efd055def9dd8573
--- /dev/null
+++ b/reports/20241018_castleMonthlyUpdates/trEPRsimul_comparisonOnZechPSI.m
@@ -0,0 +1,430 @@
+%
+% clearvars
+
+%% Easyspin simulation
+
+Sys.S = [1/2 1/2];
+Sys.g = [2.00308 2.00264 2.00226; ... % P700+
+ 2.00622 2.00507 2.00218]; % 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-
+
+Sys.J = (-1/2)*unitconvert(1e-3,'mT->MHz'); % MHz
+Sys.dip = unitconvert(0.170,'mT->MHz')/1.5; % MHz
+Sys.initState = 'singlet';
+
+fieldAxis = linspace(343.635, 347.635, 301); % GHz
+mwFreq = 9.7; % GHz
+
+Opt.GridSize = [91 0];
+Opt.GridSymmetry = '';
+
+% Field sweep
+Exp.CenterSweep = [mean(fieldAxis), max(fieldAxis) - min(fieldAxis)];
+Exp.mwFreq = mwFreq;
+Exp.Harmonic = 0;
+Exp.nPoints = numel(fieldAxis);
+
+Sys.lwpp = 0.35; % mT
+
+ySim = pepper(Sys, Exp, Opt);
+ySim = ySim/max(abs(ySim));
+
+clf
+plot(fieldAxis, ySim)
+
+%%
+
+% pathFolder = "/net/storage/gianlum33/";
+pathFolder = "S:\";
+pathFolder = pathFolder + "projects/zech_psi/data/processed/";
+fileName = "ZePSI-E-011.mat";
+
+iTime = 1000;
+importedData = load(fullfile(pathFolder, fileName));
+xdata = importedData.x{2};
+ydata = importedData.y{4}(iTime, :);
+ydata = ydata/max(abs(ydata));
+% figure()
+% clf
+hold on
+plot(xdatareal, ydata)
+hold on
+legend()
+
+%% My powder average (only for equivalent nuclei)
+
+mySys.calculategInhBroadening = 1;
+mySys.g = Sys.g;
+mySys.gFrame = Sys.gFrame;
+mySys.dip = mt2mhz(-0.170); % MHz
+mySys.J = mt2mhz(1e-3); % MHz
+mySys.eeFrame = Sys.eeFrame;
+mySys.nNuc = 3;
+mySys.A = [0, 0, 0, 9, 9.4, 12.8]; % MHz
+mySys.AFrame = [0 0 0 60 -90 0]*pi/180;
+
+mySys.rho = zeros(4, 4);
+% Up-down and down-up
+% mySys.rho(2, 2) = 1;
+% mySys.rho(3, 3) = 1;
+% Singlet
+mySys.rho(2, 2) = 1;
+mySys.rho(2, 3) = -1;
+mySys.rho(3, 2) = -1;
+mySys.rho(3, 3) = 1;
+mySys.rho = 1/2*mySys.rho;
+
+mySys.mwFreq = 9.6;
+mySys.x = xdata;
+mySys.nTheta = 40;
+mySys.nPhi = 20;
+mySys.nSamplegInh1 = 20;
+mySys.nSamplegInh2 = 20;
+mySys.coeffResLw = 0.2;
+
+mySys.lw1 = 0.67; % mT
+% mySys.lw1 = 0.1; % mT
+% mySys.lw2 = 0.25; % mT
+mySys.lw2 = 0.39; % mT
+% mySys.trlwpp1 = 0.75;
+% mySys.trlwpp2 = 0.32;
+mySys.xshift = 1.89;
+mySys.trlwpp1 = 0.7;
+mySys.trlwpp2 = 0.3;
+
+% Vary = struct( ...
+% 'lw1', mySys.lw1 - 0.005, ...
+% 'lw2', mySys.lw2 - 0.005, ...
+% 'xshift', 3);
+Vary = struct( ...
+ 'trlwpp1', mySys.trlwpp1 - 0.005, ...
+ 'trlwpp2', mySys.trlwpp2 - 0.005, ...
+ 'xshift', 3);
+
+% FitOpt.maxTime = 0.1;
+FitOpt.TolFun = 1e-4;
+FitOpt.TolEdgeLength = 1e-4;
+FitOpt.CalculateUncertainties = 1;
+FitOpt.Method = 'levmar';
+FitOpt.Gradient = 1e-4;
+FitOpt.TolStep = 1e-6;
+FitOpt.AutoScale = 'lsq';
+
+% tic
+% ySimSing = treprstickspectrum(mySys);
+% FitSing = esfit(ydata, @treprstickspectrum, {mySys}, {Vary}, FitOpt);
+% toc
+
+% mySys.rho = zeros(4, 4);
+% % (sigmaUpDown + sigmaDownUp)/2
+% mySys.rho(2, 2) = 1;
+% mySys.rho(3, 3) = 1;
+% mySys.rho = mySys.rho/2;
+
+tic
+% [ySimCiss, signalCiss] = treprstickspectrum(mySys);
+Fit = esfit(ydata, @treprstickspectrum, {mySys}, {Vary}, FitOpt);
+toc
+
+xdatareal = xdata + FitSing.pfit(3);
+% figure()
+clf
+% tiledlayout('flow', 'TileSpacing', 'compact', 'Padding', 'compact')
+% nexttile
+plot(xdatareal, ydata)
+hold on
+plot(xdatareal, FitSing.fit)
+plot(xdatareal, rescaledata(ySimCiss, ydata, 'lsq'))
+% for ii = 1:20
+% plot(xdatareal, rescaledata(signalCiss(:, ii), ydata, 'lsq'))
+% end
+% plot(xdatareal, rescaledata(ySim1 + ySim2, ydata, 'lsq'))
+% plot(xdatareal, Fit.fit)
+xlim(setaxlim(xdatareal))
+
+%%
+[thetas__, phis__] = createthetaphigrid(mySys.nTheta, mySys.nPhi);
+nVers__ = ang2vec(ones(mySys.nTheta, 1)*phis__, thetas__*ones(1, mySys.nPhi));
+
+clf
+plot(nVers__(3, 1:mySys.nTheta))
+
+%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% FIGURE castle_03_02_expDataFit160K.png
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+load("plotColors.mat")
+xplot = xdatareal;
+yplot = rescaledata(importedData.y{5}(1000, :), 'maxabs');
+
+% figure()
+clf
+plot(xplot, yplot, 'k')
+hold on
+box on
+plot(xplot, FitSing.fit, 'Color', plotColors(2, :))
+
+xlim(setaxlim(xplot, 1))
+ylim(setaxlim([yplot, FitSing.fit'], 1.05))
+xticks(340:344)
+yticks(-1:1)
+
+labelaxesfig(gca, 'Magnetic Field / mT', 'trEPR signal norm. / a.u.')
+legend("Exp. data", "Best-fit simul.")
+
+saveFolder = "S:\projects\oop_ciss_calculations/reports/20241018_castleMonthlyUpdates";
+savePath = fullfile(saveFolder, "castle_03_02_expDataFit160K.png");
+% exportgraphics(gcf, savePath, 'Resolution', 600)
+
+%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% FIGURE castle_04_05_cissPowder.png
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+load("plotColors.mat")
+xplot = xdatareal;
+fitfitplot = rescaledata(FitSing.fit, 'maxabs');
+% figure()
+clf
+plot(xdatareal, fitfitplot, 'Color', plotColors(2, :))
+hold on
+box on
+plot(xplot, rescaledata(ySimCiss, 'maxabs'), ...
+ 'Color', plotColors(6, :))
+
+xlim(setaxlim(xplot, 1))
+ylim(setaxlim([fitfitplot', rescaledata(ySimCiss, ydata, 'lsq')'], 1.05))
+xticks(340:344)
+yticks(-1:1)
+
+labelaxesfig(gca, 'Magnetic Field / mT', 'trEPR signal norm. / a.u.')
+legend("Simul. no CISS", "Simul. 100% CISS")
+
+saveFolder = "S:\projects\oop_ciss_calculations/reports/20241018_castleMonthlyUpdates";
+savePath = fullfile(saveFolder, "castle_04_05_cissPowder.png");
+% exportgraphics(gcf, savePath, 'Resolution', 600)
+
+%%
+
+function data = mytreprmf(Sys)
+
+ nSpc = numel(Sys.mwFreq_MF);
+ fields = fieldnames(Sys);
+ for ispc = 1:nspc
+ for ifield = 1:numel(fields)
+ f = fields{ifield};
+ if endsWith(f, '_MF')
+ f_ = f(1:end - 3);
+ Sys.(f_) = Sys.(f){ispc};
+ end
+ end
+ data{ispc} = mytrepr(Sys);
+ end
+end
+
+% function [ySim, gax, gInh1, gInh2, gax1, gInh11, gax2, gInh22] = mytrepr(Sys)
+function ySim = mytrepr(Sys)
+ fieldAxis = Sys.x + Sys.xshift;
+ Bmean = mean(fieldAxis);
+ dip = Sys.dip;
+ JJ = Sys.J;
+ mwFreq = Sys.mwFreq;
+ nTheta = Sys.nTheta;
+ nPhi = Sys.nPhi;
+ if isfield(Sys, 'nNuc')
+ nNuc = Sys.nNuc;
+ else
+ nNuc = 0;
+ end
+ nHfiLine = nNuc + 1;
+
+ % Parameters for g broadening
+ lw1 = mt2mhz(Sys.lw1);
+ lw2 = mt2mhz(Sys.lw2);
+ % The g-factor linewidth
+ glw1 = freq2gvalue(Bmean, lw1);
+ glw2 = freq2gvalue(Bmean, lw2);
+ % Spacing between g-values
+ dgax1 = glw1 / Sys.nSamplegInh1;
+ dgax2 = glw2 / Sys.nSamplegInh2;
+ resLw = min(lw1, lw2)*Sys.coeffResLw*1e-3; % GHz
+
+ % Grid
+ [thetas, phis] = createthetaphigrid(nTheta, nPhi);
+ nSolidAngle = nTheta*nPhi;
+
+ % Direction of B0
+ nVers = ang2vec(ones(nTheta, 1)*phis, thetas*ones(1, nPhi));
+ % Effective g-values
+ g1Tensor = rotatematrixeuangles(diag(Sys.g(1, :)), Sys.gFrame(1, :));
+ g1 = sqrt( sum( (g1Tensor*nVers).^2, 1));
+ g2Tensor = rotatematrixeuangles(diag(Sys.g(2, :)), Sys.gFrame(2, :));
+ g2 = sqrt( sum( (g2Tensor*nVers).^2, 1));
+ % Dipolar interaction
+ zD = erot(Sys.eeFrame)*[0, 0, 1]';
+ dd = dipinteraction(dip, nVers, zD);
+
+ if nNuc > 0
+ [APlus, AMinus] = calculateaplusaminusnohfi(...
+ Sys.A, Sys.AFrame, nVers);
+ else
+ APlus = zeros(nSolidAngle, 1);
+ AMinus = zeros(nSolidAngle, 1);
+ end
+
+ ySim = zeros(numel(fieldAxis), 1);
+ for ii = 1 % :nSolidAngle
+ nVers = ang2vec(ones(nTheta, 1)*phis, thetas*ones(1, nPhi));
+
+ [gax1, gInh1] = creategaxisinhbroadening(g1(ii), glw1, dgax1);
+ [gax2, gInh2] = creategaxisinhbroadening(g2(ii), glw2, dgax2);
+ ngax1 = numel(gax1);
+ ngax2 = numel(gax2);
+
+ gaussianFactor = gInh1'*gInh2;
+ gaussianFactor = gaussianFactor/sum(gaussianFactor, 'all');
+ gaussianFactor = gaussianFactor(:)';
+
+ gPlus = 1/2*(repmat(gax1', [1, ngax2]) + repmat(gax2, [ngax1, 1]));
+ gPlus = gPlus(:)';
+ gMinus = 1/2*(repmat(gax1', [1, ngax2]) - repmat(gax2, [ngax1, 1]));
+ gMinus = gMinus(:)';
+ %}
+
+ w0_ = gvalue2freq(fieldAxis', gPlus);
+ deltaw_ = gvalue2freq(fieldAxis', gMinus);
+ sintheta_ = sqrt( nVers(1, ii)^2 + nVers(2, ii)^2 );
+
+ quantNumNuc = -nNuc/2:1:nNuc/2;
+ signJmd = [-1, 1, -1, 1];
+ signOmega = [-1, -1, 1, 1];
+ signIntens = [1, -1, 1, -1];
+ if nNuc > 0
+ pascalMatrix = pascal(nHfiLine);
+ % Antidiag
+ pascalFactor = pascalMatrix(nHfiLine:nHfiLine - 1:end - 1);
+ else
+ pascalFactor = 1;
+ end
+
+ for itrans = 1:4
+ for ihfi = 1:nHfiLine
+ w0__ = w0_ + quantNumNuc(ihfi)*APlus(ii);
+ deltaw__ = deltaw_ + quantNumNuc(ihfi)*AMinus(ii);
+ Omega__ = hypot(deltaw__, JJ + dd(ii)/2);
+ wReson = w0__ + signJmd(itrans)*(JJ - dd(ii)) + ...
+ signOmega(itrans)*Omega__;
+
+ intensityReson = signIntens(itrans)*1/8*...
+ (deltaw__.^2)./(Omega__.^2);
+ trSignal__ = gaussianresonancebsweep(...
+ wReson*1e-3, mwFreq, resLw, "lwpp");
+ trSignal__ = intensityReson.*trSignal__*...
+ pascalFactor(ihfi).*gaussianFactor;
+ ySim = ySim + sum(trSignal__, 2)*sintheta_;
+ end
+ end
+ %}
+ end
+ %}
+% ySim = ySim/max(abs(ySim));
+
+end
+
+% function [ySim, signal] = treprstickspectrum(Sys)
+function ySim = treprstickspectrum(Sys)
+
+ fieldAxis = Sys.x + Sys.xshift;
+ nField = numel(fieldAxis);
+ dip = Sys.dip;
+ JJ = Sys.J;
+ trlwpp1 = Sys.trlwpp1;
+ trlwpp2 = Sys.trlwpp2;
+ mwFreq = Sys.mwFreq;
+ nTheta = Sys.nTheta;
+ nPhi = Sys.nPhi;
+ rhoInit = Sys.rho;
+ if isfield(Sys, 'nNuc')
+ nNuc = Sys.nNuc;
+ else
+ nNuc = 0;
+ end
+ nHfiLine = nNuc + 1;
+
+ % Grid
+ [thetas, phis] = createthetaphigrid(nTheta, nPhi);
+ longPhis = ones(nTheta, 1)*phis;
+ longThetas = thetas*ones(1, nPhi);
+ nSolidAngle = nTheta*nPhi;
+
+
+ % Direction of B0
+ nVers = ang2vec(longPhis, longThetas);
+ % Effective g-values
+ g1Tensor = rotatematrixeuangles(diag(Sys.g(1, :)), Sys.gFrame(1, :));
+ g1 = sqrt( sum( (g1Tensor*nVers).^2, 1));
+ g2Tensor = rotatematrixeuangles(diag(Sys.g(2, :)), Sys.gFrame(2, :));
+ g2 = sqrt( sum( (g2Tensor*nVers).^2, 1));
+ % Dipolar interaction
+ zD = erot(Sys.eeFrame)*[0, 0, 1]';
+ dd = dipinteraction(dip, nVers, zD);
+
+ % Hyperfine
+ if nNuc > 0
+ [APlus, AMinus] = calculateaplusaminusnohfi(...
+ Sys.A, Sys.AFrame, nVers);
+ else
+ APlus = zeros(nSolidAngle, 1);
+ AMinus = zeros(nSolidAngle, 1);
+ end
+
+ signal = zeros(nField, nSolidAngle);
+ parfor ii = 1:nSolidAngle
+ % disp(floor(ii/nSolidAngle*100))
+ gPlus = 1/2*(g1(ii) + g2(ii));
+ gMinus = 1/2*(g1(ii) - g2(ii));
+
+ w0_ = gvalue2freq(fieldAxis, gPlus);
+ deltaw_ = gvalue2freq(fieldAxis, gMinus);
+
+ quantNumNuc = -nNuc/2:1:nNuc/2;
+ if nNuc > 0
+ pascalMatrix = pascal(nHfiLine);
+ % Antidiag
+ pascalFactor = pascalMatrix(nHfiLine:nHfiLine - 1:end - 1);
+ pascalFactor = pascalFactor'/sum(pascalFactor);
+ else
+ pascalFactor = 1;
+ end
+
+
+ for itrans = 1:4
+
+ for ihfi = 1:nHfiLine
+ w0__ = w0_ + quantNumNuc(ihfi)*APlus(ihfi);
+ deltaw__ = deltaw_ + quantNumNuc(ihfi)*AMinus(ihfi);
+ Omega__ = hypot(deltaw__, JJ + dd(ii)/2);
+ pChar = mean(0.5*(1 + deltaw__./Omega__));
+ if itrans == 3 || itrans == 4
+ trlwpp = pChar*trlwpp1 + (1. - pChar)*trlwpp2;
+ else
+ trlwpp = pChar*trlwpp2 + (1. - pChar)*trlwpp1;
+ end
+ wReson = myeigenenergies(w0__, deltaw__, JJ, dd(ii), itrans);
+
+ intensityReson = intensityresonance(rhoInit, nVers(:, ii), deltaw__, JJ, dd(ii), itrans);
+ signal__ = gaussianresonancebsweep( ...
+ wReson*1e-3, mwFreq, mt2mhz(trlwpp)*1e-3, "lwpp");
+ signal(:, ii) = signal(:, ii) + intensityReson'.*signal__'*pascalFactor(ihfi);
+
+ end
+ end
+ end
+ ySim = averageoversolidangle(signal, nTheta, nPhi, 2);
+ ySim = real(ySim);
+end
+
+
diff --git a/reports/20241018_castleMonthlyUpdates/trEPRsimul_newVsOld.m b/reports/20241018_castleMonthlyUpdates/trEPRsimul_newVsOld.m
new file mode 100644
index 0000000000000000000000000000000000000000..b4eddd23ec37440538be17e9857906a628389a2c
--- /dev/null
+++ b/reports/20241018_castleMonthlyUpdates/trEPRsimul_newVsOld.m
@@ -0,0 +1,362 @@
+%
+saveFolder = "/net/storage/gianlum33/projects/oop_ciss_calculations/reports/20241018_castleMonthlyUpdates";
+load('plotColors.mat')
+
+FLAG_CALCULATE_INH_BROAD = 1;
+FLAG_SAVE_FIG = 1;
+
+g1 = 2.006;
+g2 = 2.002;
+bfieldg1 = mhz2mt(9.7, g1)*1e3;
+bfieldg2 = mhz2mt(9.7, g2)*1e3;
+
+d0 = mhz2mt(-10); % mT
+J = 0;
+
+nPoint = 1500;
+bWidth = 4.5; % mT
+
+gInhLwSmall = 0.1;
+gInhLwLarge = 1;
+trlwppSmall = gInhLwSmall - 0.01;
+trlwppLarge = gInhLwLarge - 0.2;
+
+gInhLw = gInhLwSmall;
+% gInhLw = gInhLwLarge;
+trlwpp = trlwppSmall;
+% trlwpp = trlwppLarge;
+
+bCenter = 1/2*(bfieldg2 + bfieldg1);
+bDiff = 1/2*(bfieldg2 - bfieldg1);
+
+Omega = hypot(bDiff, (J + d0/2))'; % mT
+bMin = bCenter - bWidth/2;
+bMax = bCenter + bWidth/2;
+xx = linspace(bMin, bMax, nPoint);
+
+% Find bres1 and bres2 in xx
+[~, ires1] = min(abs(xx - bfieldg1));
+[~, ires2] = min(abs(xx - bfieldg2));
+
+
+fig = figure(1);
+clf
+plot([xx(ires1), xx(ires2)], [0, 0], 'xk', 'HandleVisibility', 'off')
+box on
+yline(0, 'HandleVisibility', 'off')
+xline([bfieldg1, bfieldg2], ':', 'HandleVisibility', 'off')
+
+ann1 = "g = " + string(g1);
+xann1 = 344.6;
+xann2 = 346.25;
+yann = -0.1;
+text(xann1, yann, ann1, 'FontSize', 14)
+ann2 = "g = " + string(g2);
+text(xann2, yann, ann2, 'FontSize', 14)
+
+xlim(setaxlim(xx))
+xticks(343:1:350)
+yticks(-1:0.5:1)
+labelaxesfig(gca, 'Magnetic field / mT', 'trEPR intensity / a.u.')
+
+if gInhLw == gInhLwSmall && FLAG_SAVE_FIG
+ savePath = fullfile(saveFolder, "trEPRsimul_newVsOld_01.png");
+ exportgraphics(fig, savePath, 'Resolution', 600)
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% FIGURE 2
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+bres1 = bCenter - (J - d0) - Omega;
+bres2 = bCenter + (J - d0) - Omega;
+bres3 = bCenter - (J - d0) + Omega;
+bres4 = bCenter + (J - d0) + Omega;
+
+fig = figure(2);
+clf
+plot([xx(ires1), xx(ires2)], [0, 0], 'xk', 'HandleVisibility', 'off')
+box on
+yline(0, 'HandleVisibility', 'off')
+xline([bres1, bres2, bres3, bres4], '--', 'HandleVisibility', 'off')
+
+labelaxesfig(gca, 'Magnetic field / mT', 'trEPR intensity / a.u.')
+ax = gca;
+axCenter = ax.Position(1) + ax.Position(3)/2;
+xOmega = [-1, 1]*Omega/bWidth*ax.Position(3) + axCenter;
+xDiff1 = [-1, 1]*d0/bWidth*ax.Position(3) + xOmega(1);
+xDiff2 = [-1, 1]*d0/bWidth*ax.Position(3) + xOmega(2);
+annotation('doublearrow', xOmega, [.45, .45])
+annotation('doublearrow', xDiff1, [.7, .7])
+annotation('doublearrow', xDiff2, [.7, .7])
+
+xlim(setaxlim(xx))
+xticks(343:1:350)
+yticks(-1:0.5:1)
+
+if gInhLw == gInhLwSmall && FLAG_SAVE_FIG
+ savePath = fullfile(saveFolder, "trEPRsimul_newVsOld_02.png");
+ exportgraphics(fig, savePath, 'Resolution', 600)
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% FIGURE 3
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+scale = 0.75;
+yspc1 = scale*gaussian(xx, bres1, trlwpp)/max(gaussian(xx, bres1, trlwpp));
+yspc2 = scale*gaussian(xx, bres2, trlwpp)/max(gaussian(xx, bres2, trlwpp));
+yspc3 = scale*gaussian(xx, bres3, trlwpp)/max(gaussian(xx, bres3, trlwpp));
+yspc4 = scale*gaussian(xx, bres4, trlwpp)/max(gaussian(xx, bres4, trlwpp));
+yspc = (yspc1 - yspc2 + yspc3 - yspc4);
+fig = figure(3);
+clf
+plot([xx(ires1), xx(ires2)], [0, 0], 'xk', 'HandleVisibility', 'off')
+box on
+hold on
+plot(0, 0, 'HandleVisibility', 'off') % Chaotic way to skip one color
+plot(xx, yspc1, 'DisplayName', 'Reson. 1')
+plot(xx, -yspc2, 'DisplayName', 'Reson. 2')
+plot(xx, yspc3, 'DisplayName', 'Reson. 3')
+plot(xx, -yspc4, 'DisplayName', 'Reson. 4')
+plot(xx, yspc, 'k', 'DisplayName', 'Sum')
+% yline(0, 'HandleVisibility', 'off')
+xline([bres1, bres2, bres3, bres4], '--', 'HandleVisibility', 'off')
+
+labelaxesfig(gca, 'Magnetic field / mT', 'trEPR intensity / a.u.')
+legend()
+
+xlim(setaxlim(xx))
+ylim([-1, 1])
+xticks(343:1:350)
+yticks(-1:0.5:1)
+
+if gInhLw == gInhLwSmall && FLAG_SAVE_FIG
+ savePath = fullfile(saveFolder, "trEPRsimul_newVsOld_03.png");
+ exportgraphics(fig, savePath, 'Resolution', 600)
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% FIGURE 4
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ scale = 0.75;
+
+fig = figure(4);
+clf
+plot([xx(ires1), xx(ires2)], [0, 0], 'xk', 'HandleVisibility', 'off')
+box on
+hold on
+plot(xx, yspc, 'k', 'DisplayName', 'Method 1')
+% yline(0, 'HandleVisibility', 'off')
+
+labelaxesfig(gca, 'Magnetic field / mT', 'trEPR intensity / a.u.')
+legend()
+
+xlim(setaxlim(xx))
+ylim([-1, 1])
+xticks(343:1:350)
+yticks(-1:0.5:1)
+
+if gInhLw == gInhLwSmall && FLAG_SAVE_FIG
+ savePath = fullfile(saveFolder, "trEPRsimul_newVsOld_04.png");
+ exportgraphics(fig, savePath, 'Resolution', 600)
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% FIGURE 5 AND FIGURE 9
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+bgInh1 = gaussian(xx, bfieldg1, gInhLw);
+bgInh1 = bgInh1/max(bgInh1)*0.5;
+bgInh2 = gaussian(xx, bfieldg2, gInhLw);
+bgInh2 = bgInh2/max(bgInh2)*0.5;
+
+fig = figure(5);
+clf
+plot([xx(ires1), xx(ires2)], [0, 0], 'xk', 'HandleVisibility', 'off')
+box on
+hold on
+plot(xx, bgInh1, '--', 'Color', plotColors(2, :), 'DisplayName', 'g-value distr.')
+plot(xx, bgInh2, '--', 'Color', plotColors(2, :), 'HandleVisibility', 'off')
+yline(0, 'HandleVisibility', 'off')
+xline([bfieldg1, bfieldg2], ':', 'HandleVisibility', 'off')
+
+text(xann1, yann, ann1, 'FontSize', 14)
+text(xann2, yann, ann2, 'FontSize', 14)
+
+xlim(setaxlim(xx))
+ylim([-1, 1])
+xticks(343:1:350)
+yticks(-1:0.5:1)
+labelaxesfig(gca, 'Magnetic field / mT', 'trEPR intensity / a.u.')
+legend()
+
+if FLAG_SAVE_FIG
+ if gInhLw == gInhLwSmall
+ savePath = fullfile(saveFolder, "trEPRsimul_newVsOld_05.png");
+ exportgraphics(fig, savePath, 'Resolution', 600)
+ elseif gInhLw == gInhLwLarge
+ savePath = fullfile(saveFolder, "trEPRsimul_newVsOld_09.png");
+ exportgraphics(fig, savePath, 'Resolution', 600)
+ end
+end
+
+%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% FIGURE 6
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+iCut = bgInh1 > 0.045;
+bgInh1o = bgInh1(iCut);
+xx1o = xx(iCut);
+iCut = 1:2:numel(xx1o);
+bgInh1o = bgInh1o(iCut);
+xx1o = xx1o(iCut);
+iCut = 1:5:numel(bgInh1o);
+bgInh1oPlot = bgInh1o(iCut);
+xx1oPlot = xx1o(iCut);
+
+iCut = bgInh2 > 0.045;
+bgInh2o = bgInh2(iCut);
+xx2o = xx(iCut);
+iCut = 1:2:numel(xx2o);
+bgInh2o = bgInh2o(iCut);
+xx2o = xx2o(iCut);
+iCut = 1:5:numel(bgInh2o);
+bgInh2oPlot = bgInh2o(iCut);
+xx2oPlot = xx2o(iCut);
+
+fig = figure(6);
+clf
+plot([xx(ires1), xx(ires2)], [0, 0], 'xk', 'HandleVisibility', 'off')
+box on
+hold on
+plot(xx, bgInh1, '--', 'Color', plotColors(2, :), 'DisplayName', 'g-value distr.')
+plot(xx, bgInh2, '--', 'Color', plotColors(2, :), 'HandleVisibility', 'off')
+plot(xx1oPlot, bgInh1oPlot, 'o', 'Color', plotColors(2, :), 'DisplayName', 'g-value sampling')
+plot(xx2oPlot, bgInh2oPlot, 'o', 'Color', plotColors(2, :), 'HandleVisibility', 'off')
+yline(0, 'HandleVisibility', 'off')
+xline([bfieldg1, bfieldg2], ':', 'HandleVisibility', 'off')
+
+text(xann1, yann, ann1, 'FontSize', 14)
+text(xann2, yann, ann2, 'FontSize', 14)
+
+xlim(setaxlim(xx))
+ylim([-1, 1])
+xticks(343:1:350)
+yticks(-1:0.5:1)
+labelaxesfig(gca, 'Magnetic field / mT', 'trEPR intensity / a.u.')
+legend()
+
+if gInhLw == gInhLwSmall && FLAG_SAVE_FIG
+ savePath = fullfile(saveFolder, "trEPRsimul_newVsOld_06.png");
+ exportgraphics(fig, savePath, 'Resolution', 600)
+end
+
+%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% FIGURE 7
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if gInhLw == gInhLwSmall
+ smallLw = gInhLw/5;
+else
+ smallLw = gInhLw/10;
+end
+gaussianFac = bgInh1o'*bgInh2o;
+if FLAG_CALCULATE_INH_BROAD
+ yspcInh = zeros(size(yspc));
+ for ii = 1:1
+ for jj = 1:1
+ bCenter_ = 1/2*(xx1o(ii) + xx2o(jj));
+ bDiff_ = 1/2*(xx1o(ii) - xx2o(jj));
+
+ Omega_ = hypot(bDiff_, (J + d0/2))'; % mT
+
+ bres1_ = bCenter_ - (J - d0) - Omega_;
+ bres2_ = bCenter_ + (J - d0) - Omega_;
+ bres3_ = bCenter_ - (J - d0) + Omega_;
+ bres4_ = bCenter_ + (J - d0) + Omega_;
+
+ yspc_ = gaussian(xx, bres1_, smallLw)/max(gaussian(xx, bres1_, smallLw));
+ yspc_ = yspc_ - gaussian(xx, bres2_, smallLw)/max(gaussian(xx, bres2_, smallLw));
+ yspc_ = yspc_ + gaussian(xx, bres3_, smallLw)/max(gaussian(xx, bres3_, smallLw));
+ yspc_ = yspc_ - gaussian(xx, bres4_, smallLw)/max(gaussian(xx, bres4_, smallLw));
+ yspcInh = yspc_*gaussianFac(ii, jj);
+
+ end
+ end
+end
+
+fig = figure(7);
+clf
+plot([xx(ires1), xx(ires2)], [0, 0], 'xk', 'HandleVisibility', 'off')
+box on
+hold on
+plot(xx, yspcInh/max(yspcInh)*scale, 'Color', plotColors(2, :))
+xline([bres1_, bres2_, bres3_, bres4_], '--', 'HandleVisibility', 'off')
+
+xlim(setaxlim(xx))
+ylim([-1, 1])
+xticks(343:1:350)
+yticks(-1:0.5:1)
+labelaxesfig(gca, 'Magnetic field / mT', 'trEPR intensity / a.u.')
+
+if gInhLw == gInhLwSmall && FLAG_SAVE_FIG
+ savePath = fullfile(saveFolder, "trEPRsimul_newVsOld_07.png");
+ exportgraphics(fig, savePath, 'Resolution', 600)
+end
+
+%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% FIGURE 8 AND 10
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+fig = figure(8);
+
+clf
+gaussianFac = bgInh1o'*bgInh2o;
+if FLAG_CALCULATE_INH_BROAD
+ yspcInh = zeros(size(yspc));
+ for ii = 1:numel(xx1o)
+ for jj = 1:numel(xx2o)
+ bCenter_ = 1/2*(xx1o(ii) + xx2o(jj));
+ bDiff_ = 1/2*(xx1o(ii) - xx2o(jj));
+
+ Omega_ = hypot(bDiff_, (J + d0/2))'; % mT
+
+ bres1_ = bCenter_ - (J - d0) - Omega_;
+ bres2_ = bCenter_ + (J - d0) - Omega_;
+ bres3_ = bCenter_ - (J - d0) + Omega_;
+ bres4_ = bCenter_ + (J - d0) + Omega_;
+
+ yspc_ = gaussian(xx, bres1_, smallLw)/max(gaussian(xx, bres1_, smallLw));
+ yspc_ = yspc_ - gaussian(xx, bres2_, smallLw)/max(gaussian(xx, bres2_, smallLw));
+ yspc_ = yspc_ + gaussian(xx, bres3_, smallLw)/max(gaussian(xx, bres3_, smallLw));
+ yspc_ = yspc_ - gaussian(xx, bres4_, smallLw)/max(gaussian(xx, bres4_, smallLw));
+ yspcInh = yspcInh+yspc_*gaussianFac(ii, jj);
+ end
+ end
+end
+
+plot([xx(ires1), xx(ires2)], [0, 0], 'xk', 'HandleVisibility', 'off')
+box on
+hold on
+plot(xx, yspcInh/max(yspcInh)*scale, 'k', 'DisplayName', 'Method 1')
+plot(xx, yspc/max(yspc)*scale, 'Color', plotColors(2, :), 'DisplayName', 'Method 2')
+
+xlim(setaxlim(xx))
+ylim([-1, 1])
+xticks(343:1:350)
+yticks(-1:0.5:1)
+labelaxesfig(gca, 'Magnetic field / mT', 'trEPR intensity / a.u.')
+legend()
+
+if FLAG_SAVE_FIG
+ if gInhLw == gInhLwSmall
+ savePath = fullfile(saveFolder, "trEPRsimul_newVsOld_08.png");
+ exportgraphics(fig, savePath, 'Resolution', 600)
+ elseif gInhLw == gInhLwLarge
+ savePath = fullfile(saveFolder, "trEPRsimul_newVsOld_10.png");
+ exportgraphics(fig, savePath, 'Resolution', 600)
+ end
+end
\ No newline at end of file
diff --git a/reports/20241018_castleMonthlyUpdates/zechPSI_trEPRsigAmp_vsTemperature_01.png b/reports/20241018_castleMonthlyUpdates/zechPSI_trEPRsigAmp_vsTemperature_01.png
new file mode 100644
index 0000000000000000000000000000000000000000..f88901ac6558eded4e5da6a087279a1cc0989d5e
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/zechPSI_trEPRsigAmp_vsTemperature_01.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/zechPSI_trEPRsigAmp_vsTemperature_02.png b/reports/20241018_castleMonthlyUpdates/zechPSI_trEPRsigAmp_vsTemperature_02.png
new file mode 100644
index 0000000000000000000000000000000000000000..67ecc79e980b20d1626110d87805932d129ecc3e
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/zechPSI_trEPRsigAmp_vsTemperature_02.png differ
diff --git a/reports/20241018_castleMonthlyUpdates/zechPSI_trEPRsigAmp_vsTemperature_03.png b/reports/20241018_castleMonthlyUpdates/zechPSI_trEPRsigAmp_vsTemperature_03.png
new file mode 100755
index 0000000000000000000000000000000000000000..86477049a7a0f5263602d6fd3ebd74d1ba3f2f91
Binary files /dev/null and b/reports/20241018_castleMonthlyUpdates/zechPSI_trEPRsigAmp_vsTemperature_03.png differ
diff --git a/util/averageoversolidangle.m b/util/averageoversolidangle.m
index 8562386d0a8df6ec390bb9fc77c18de13f5beda7..f54ea0c78ace5fd4315b144fa53d576499c10cd7 100644
--- a/util/averageoversolidangle.m
+++ b/util/averageoversolidangle.m
@@ -1,4 +1,10 @@
function outV = averageoversolidangle(inV, nTheta, nPhi, dim)
+arguments
+ inV
+ nTheta
+ nPhi
+ dim
+end
[thetas, ~] = createthetaphigrid(nTheta, nPhi);
diff --git a/util/averageoversolidanglemask.m b/util/averageoversolidanglemask.m
new file mode 100755
index 0000000000000000000000000000000000000000..312d4b11a4c62dce7da6c50ba8e1adcb058c673b
--- /dev/null
+++ b/util/averageoversolidanglemask.m
@@ -0,0 +1,73 @@
+function outV = averageoversolidanglemask(inV, nTheta, nPhi, dim, mask)
+arguments
+ inV
+ nTheta
+ nPhi
+ dim
+ mask = []
+end
+
+mask = checkparams(inV, nTheta, nPhi, dim, mask);
+sizeV = size(inV);
+[thetas, ~] = createthetaphigrid(nTheta, nPhi);
+
+if numel(dim) == 1
+ % Put dim as first dimension
+ idxNonAvg = [1:(dim - 1), (dim + 1):numel(sizeV)];
+ inV = permute(inV, [dim, idxNonAvg]); % Size is (nTheta*nPhi, ...)
+
+ % Create logical mask
+ logicMask = zeros(1, nTheta*nPhi);
+ logicMask(mask) = 1;
+ logicMask = boolean(logicMask);
+
+ % Create weight
+ thetas = repmat(thetas, [nPhi, 1]);
+ weight = sin(thetas); % Size is (nTheta*nPhi, 1)
+ weight = weight(logicMask)/sum(weight(logicMask));
+
+ % Average
+ outV = sum(weight.*inV(logicMask, :), 1);
+
+ % Reshape
+ if numel(size(outV)) > 2
+ outV = reshape(outV, sizeV(idxNonAvg));
+ end
+else
+% weight = sin(thetas)/sum(sin(thetas)*ones(1, nPhi), 'all');
+%
+% weight = reshape(weight, [ones(1, min(dim(1) - 1, 1)), nTheta]);
+% outV = sum(weight(mask).*inV(mask), dim);
+% outV = squeeze(outV);
+ error('Still needs to be implemented')
+end
+end
+
+function mask = checkparams(inV, nTheta, nPhi, dim, mask)
+
+% Check numel(dim)
+switch numel(dim)
+ case 1
+ if size(inV, dim) ~= nTheta*nPhi
+ error("size(inV, dim) must equal nTheta*nPhi.")
+ end
+ case 2
+ if size(inV, dim(1)) ~= nTheta
+ error("size(inV, dim(1)) must equal nTheta.")
+ elseif size(inV, dim(2)) ~= nPhi
+ error("size(inV, dim(2)) must equal nPhi.")
+ end
+ otherwise
+ error("numel(dim) must be equal to 1 or 2.")
+end
+
+% Check mask
+if isempty(mask)
+ if numel(dim) == 1
+ mask = boolean(ones(1, nTheta*nPhi));
+ else
+ mask = boolean(ones(nTheta, nPhi));
+ end
+end
+
+end
\ No newline at end of file
diff --git a/util/gaussianresonancebsweep.m b/util/gaussianresonancebsweep.m
index 43a73daaa856270cb231d704a70d480baaaea630..e242575be1d7dfee1d09695792ddb03fcd54c184 100755
--- a/util/gaussianresonancebsweep.m
+++ b/util/gaussianresonancebsweep.m
@@ -2,7 +2,7 @@ function y1 = gaussianresonancebsweep(wReson, mwFreq, c, mode)
arguments
wReson
mwFreq
- c (1, 1) double
+ c
mode string = "lwpp"
end
diff --git a/util/intensityresonance.m b/util/intensityresonance.m
new file mode 100644
index 0000000000000000000000000000000000000000..7f3bb43be53ff9ba18e181b1138e2c5364134bb9
--- /dev/null
+++ b/util/intensityresonance.m
@@ -0,0 +1,98 @@
+function intensityResonance = intensityresonance(rhoInit, nVers, deltaw, J, d, itrans)
+
+ % Rotate rhoInit depending on theta, phi
+ rho0 = rotatezfirstysecond(rhoInit, nVers);
+
+ % rho should be in the product basis
+ rho1 = rotateproduct2coupled(rho0);
+ % rho0 = rhoInit;
+
+ alpha = 1/2*atan(deltaw./(J + d/2));
+
+ transProb = [sin(alpha).^2/2; ... % 1 to 2
+ cos(alpha).^2/2; ... ss % 3 to 4
+ cos(alpha).^2/2; ... % 1 to 3
+ sin(alpha).^2/2]; % 2 to 4
+
+ UMatrix = diagonalizingmatrix(alpha);
+ rho = pagemtimes(UMatrix, pagemtimes(rho1, pagetranspose(UMatrix)));
+ populs = mydiag3d(rho);
+ populDiff = [populs(1, :) - populs(2, :); ...
+ populs(3, :) - populs(4, :); ...
+ populs(1, :) - populs(3, :); ...
+ populs(2, :) - populs(4, :)];
+
+ if abs(nVers(3)) < 0.1 && itrans == 1
+% nVers
+% rho0
+% rho1
+% trace(rho(:, :, 1))
+% populDiff(:, 1)
+ end
+ intensityResonance = -populDiff(itrans, :).*transProb(itrans, :);
+end
+
+function rho = rotatezfirstysecond(rhoInit, nVers)
+
+ theta = atan(hypot(nVers(1), nVers(2))/nVers(3));
+ phi = atan(nVers(2)/nVers(1));
+
+ szapszb = 1/2*diag([1, 0, 0, -1]);
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%
+ % syapsyb =
+ % 0 -i/2 -i/2 0
+ % i/2 0 0 -i/2
+ % i/2 0 0 -i/2
+ % 0 i/2 i/2 0
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ syapsyb = zeros(4, 4);
+ syapsyb(1, 2:3) = [-1, -1];
+ syapsyb(2, [1, 4]) = [1, -1];
+ syapsyb(3, [1, 4]) = [1, -1];
+ syapsyb(4, 2:3) = [1, 1];
+ syapsyb = syapsyb*1i/2;
+
+ rotz = expm(1i*szapszb*(-phi)); % Rot 1: around z of minus phi (phi!)
+ rotzT = expm(-1i*szapszb*(-phi));
+ roty = expm(1i*syapsyb*(-theta)); % Rot 2: around y of minus theta
+ rotyT = expm(-1i*syapsyb*(-theta));
+
+ rho = rotzT*rhoInit*rotz;
+ rho = rotyT*rho*roty;
+
+end
+
+function rho = rotateproduct2coupled(rhoInit)
+ % Transform rho to the T+, S, T0, T- basis
+ rotProd2Coupled = diag([1, 1/sqrt(2), 1/sqrt(2), 1]);
+ rotProd2Coupled(2, 3) = 1/sqrt(2);
+ rotProd2Coupled(3, 2) = -1/sqrt(2);
+
+ rho = rotProd2Coupled'*rhoInit*rotProd2Coupled;
+end
+
+function U = diagonalizingmatrix(alpha)
+
+ % U = [1 0 0 0;
+ % 0 cos(alpha) sin(alpha) 0;
+ % 0 -sin(alpha) cos(alpha) 0;
+ % 0 0 0 1];
+
+ U = zeros(4, 4, numel(alpha));
+ U(1, 1, :) = 1;
+ U(2, 2, :) = cos(alpha);
+ U(2, 3, :) = sin(alpha);
+ U(3, 2, :) = -sin(alpha);
+ U(3, 3, :) = cos(alpha);
+ U(4, 4, :) = 1;
+
+end
+
+function diag = mydiag3d(matr)
+
+ diag = squeeze(matr(1, :, :));
+ for ii = 1:4 % 4 = first and second dimension of matr
+ diag(ii, :) = matr(ii, ii, :);
+ end
+end
\ No newline at end of file
diff --git a/util/myeigenenergies.m b/util/myeigenenergies.m
new file mode 100644
index 0000000000000000000000000000000000000000..e85383e13106fef6c76ec24dc09b7345d0d3fc24
--- /dev/null
+++ b/util/myeigenenergies.m
@@ -0,0 +1,9 @@
+function energ = myeigenenergies(w0, deltaw, J, d, itrans)
+
+ signJmd = [-1, 1, -1, 1];
+ signOmega = [-1, -1, 1, 1];
+
+ Omega = hypot(deltaw, J + d/2);
+ energ = w0 + signJmd(itrans)*(J - d) + signOmega(itrans)*Omega;
+end
+
diff --git a/zech_dPSItrEPR_gInhBroadeningVsUsual.m b/zech_dPSItrEPR_gInhBroadeningVsUsual.m
new file mode 100644
index 0000000000000000000000000000000000000000..6bd33ac363afc1d8c44e01265bf84d61627bc309
--- /dev/null
+++ b/zech_dPSItrEPR_gInhBroadeningVsUsual.m
@@ -0,0 +1,278 @@
+%
+clearvars
+
+%% Easyspin simulation
+
+% Taken from Easyspin example
+Sys.S = [1/2 1/2];
+Sys.g = [2.00308 2.00264 2.00226; ... % P700+
+ 2.00622 2.00507 2.00218]; % 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-
+
+Sys.J = (-1/2)*unitconvert(1e-3,'mT->MHz'); % MHz
+Sys.dip = unitconvert(0.170,'mT->MHz')/1.5; % MHz
+Sys.initState = 'singlet';
+
+fieldAxis = linspace(343.635, 347.635, 301); % GHz
+mwFreq = 9.7; % GHz
+
+Opt.GridSize = [91 0];
+Opt.GridSymmetry = '';
+
+% Field sweep
+Exp.CenterSweep = [mean(fieldAxis), max(fieldAxis) - min(fieldAxis)];
+Exp.mwFreq = mwFreq;
+Exp.Harmonic = 0;
+Exp.nPoints = numel(fieldAxis);
+
+Sys.lwpp = 0.35; % mT
+
+ySim = pepper(Sys, Exp, Opt);
+ySim = ySim/max(abs(ySim));
+
+clf
+plot(fieldAxis, ySim)
+
+%%
+
+% pathFolder = "/net/storage/gianlum33/";
+pathFolder = "S:\";
+pathFolder = pathFolder + "projects/oop_ciss_calculations/data/digitized/";
+pathX = pathFolder + "expData_zech_aStructuralModelFor_xBandDeut.mat";
+pathQ = pathFolder + "expData_zech_aStructuralModelFor_qBandDeut.mat";
+pathW = pathFolder + "expData_zech_aStructuralModelFor_wBandDeut.mat";
+
+importedData = load(pathX);
+xdataXdeut = importedData.x;
+ydataXdeut = importedData.y;
+ydataXdeut = ydataXdeut - (ydataXdeut(1) + ydataXdeut(end))/2;
+importedData = load(pathQ);
+xdataQdeut = importedData.x;
+ydataQdeut = importedData.y;
+ydataQdeut = ydataQdeut - (ydataQdeut(1) + ydataQdeut(end))/2;
+importedData = load(pathW);
+xdataWdeut = importedData.x;
+ydataWdeut = importedData.y;
+
+importedData = load('S:\projects\zech_psi\data\processed\ZePSI-E-007015.mat');
+% clf
+% h = ScrollableAxes();
+% plot(h, x{2}, x{1}, y');
+xdataX = importedData.x{2}'/10; % mT
+ydataX = importedData.y(815, :);
+ydataX = ydataX/max(abs(ydataX));
+pathQ2 = pathFolder + "expData_zech_aStructuralModelFor_qBand.mat";
+importedData = load(pathQ2);
+xdataQ = importedData.x;
+ydataQ = importedData.y;
+
+%% My powder average (only for equivalent nuclei)
+
+mySys.g = Sys.g;
+mySys.gFrame = Sys.gFrame;
+mySys.dip = mt2mhz(-0.170); % MHz
+mySys.J = mt2mhz(1e-3); % MHz
+mySys.eeFrame = Sys.eeFrame;
+mySys.nNuc = 0;
+mySys.A = [0, 0, 0, 9, 9.4, 12.8]; % MHz
+mySys.AFrame = [0 0 0 60 -90 0]*pi/180;
+mySys.lw1 = 0.3; % mT
+% mySys.lw1 = 0.7499; % mT
+% mySys.lw2 = 0.25; % mT
+mySys.lw2 = 0.3; % mT
+
+% FitOpt.maxTime = 0.1;
+FitOpt.TolFun = 1e-4;
+FitOpt.TolEdgeLength = 1e-4;
+FitOpt.CalculateUncertainties = 0;
+FitOpt.Method = 'levmar';
+FitOpt.Gradient = 1e-3;
+FitOpt.AutoScale = 'lsq';
+
+Vary.lw1 = mySys.lw1 - 0.005;
+Vary.lw2 = mySys.lw2 - 0.005;
+Vary.xshift = 0.1;
+
+ydata = {ydataXdeut, ydataQdeut, ydataWdeut};
+mySys.mwFreq_MF = {9.7, 34, 94.11};
+ % xdataX + 5.5, xdataQ, xdataW + 3.5
+mySys.xshift = -0.0245;
+mySys.x_MF = {xdataXdeut - 1, xdataQdeut, xdataWdeut + 0.};
+mySys.nTheta_MF = {40, 20, 50}; % {40, 40, 20}; %, 50};
+mySys.nPhi_MF = {20, 10, 30}; % {20, 20, 20}; %, 20};
+mySys.nSamplegInh1_MF = {40, 20, 20};
+mySys.nSamplegInh2_MF = {40, 20, 20};
+mySys.coeffResLw_MF = {0.05, 0.1, 0.1};
+
+nSim = 15;
+xshifts = linspace(-0.04, 0.04, nSim);
+wb = waitbar(0, 'Simulting');
+ydata = {ydataXdeut};
+for ii = 1
+ waitbar(ii/nSim, wb, ...
+ "Simulation " + string(ii) + " out of " + string(nSim) + "...")
+ mySys.xshift = -0.0245;
+ tic
+ ySim = mytreprmf(mySys);
+ Fit = esfit(ydata, @mytreprmf, {mySys}, {Vary}, FitOpt);
+ toc
+
+ figure()
+ clf
+ tiledlayout('flow', 'TileSpacing', 'compact', 'Padding', 'compact')
+ for ifig = 1
+ x_ = mySys.x_MF{ifig};
+ nexttile
+ plot(x_, ydata{ifig})
+ hold on
+ plot(x_ + mySys.xshift, Fit.fit)
+ plot(x_ + mySys.xshift, rescaledata(ySim{ifig}, ydata{ifig}, 'lsq'))
+% xlim(setaxlim(mySys.x_MF{ii}))
+% xlim(setaxlim(1:numel(ydata{ii})))
+% ylim(setaxlim(ydata{ii}, 1.05))
+ title(string(xshifts(ii)) + " : " + ...
+ string(Fit.pfit(1)) + " " + string(Fit.pfit(2)))
+
+ end
+end
+close(wb); clear('wbar')
+% xline(mhz2mt(mySys.mwFreq_MF{3}, mySys.g(1))*1e3, 'r')
+% xline(mhz2mt(mySys.mwFreq_MF{3}, mySys.g(3))*1e3, 'r')
+% xline(mhz2mt(mySys.mwFreq_MF{3}, mySys.g(5))*1e3, 'r')
+% xline(mhz2mt(mySys.mwFreq_MF{3}, mySys.g(2))*1e3, 'b')
+% xline(mhz2mt(mySys.mwFreq_MF{3}, mySys.g(4))*1e3, 'b')
+% xline(mhz2mt(mySys.mwFreq_MF{3}, mySys.g(6))*1e3, 'b')
+% tic
+% ydata = {yX, yQ, yW};
+% Fit = esfit(ydata, @mytreprmf, {mySys}, {Vary}, FitOpt);
+
+%%
+ySim0 = mytreprmf(mySys);
+%%
+
+function data = mytreprmf(Sys)
+
+ nSpc = numel(Sys.mwFreq_MF);
+ fields = fieldnames(Sys);
+ for ispc = 1
+ for ifield = 1:numel(fields)
+ f = fields{ifield};
+ if endsWith(f, '_MF')
+ f_ = f(1:end - 3);
+ Sys.(f_) = Sys.(f){ispc};
+ end
+ end
+ data{ispc} = mytrepr(Sys);
+ end
+end
+
+% function [ySim, gax, gInh1, gInh2, gax1, gInh11, gax2, gInh22] = mytrepr(Sys)
+function ySim = mytrepr(Sys)
+ fieldAxis = Sys.x + Sys.xshift;
+ Bmean = mean(fieldAxis);
+ dip = Sys.dip;
+ JJ = Sys.J;
+ mwFreq = Sys.mwFreq;
+ nTheta = Sys.nTheta;
+ nPhi = Sys.nPhi;
+ if isfield(Sys, 'nNuc')
+ nNuc = Sys.nNuc;
+ else
+ nNuc = 0;
+ end
+ nHfiLine = nNuc + 1;
+
+ % Parameters for g broadening
+ lw1 = mt2mhz(Sys.lw1);
+ lw2 = mt2mhz(Sys.lw2);
+ % The g-factor linewidth
+ glw1 = freq2gvalue(Bmean, lw1);
+ glw2 = freq2gvalue(Bmean, lw2);
+ % Spacing between g-values
+ dgax1 = glw1 / Sys.nSamplegInh1;
+ dgax2 = glw2 / Sys.nSamplegInh2;
+ resLw = min(lw1, lw2)*Sys.coeffResLw*1e-3; % GHz
+
+ % Grid
+ [thetas, phis] = createthetaphigrid(nTheta, nPhi);
+ nSolidAngle = nTheta*nPhi;
+
+ % Direction of B0
+ nVers = ang2vec(ones(nTheta, 1)*phis, thetas*ones(1, nPhi));
+ % Effective g-values
+ g1Tensor = rotatematrixeuangles(diag(Sys.g(1, :)), Sys.gFrame(1, :));
+ g1 = sqrt( sum( (g1Tensor*nVers).^2, 1));
+ g2Tensor = rotatematrixeuangles(diag(Sys.g(2, :)), Sys.gFrame(2, :));
+ g2 = sqrt( sum( (g2Tensor*nVers).^2, 1));
+ % Dipolar interaction
+ zD = erot(Sys.eeFrame)*[0, 0, 1]';
+ dd = dipinteraction(dip, nVers, zD);
+
+ if nNuc > 0
+ [APlus, AMinus] = calculateaplusaminusnohfi(...
+ Sys.A, Sys.AFrame, nVers);
+ else
+ APlus = zeros(nSolidAngle, 1);
+ AMinus = zeros(nSolidAngle, 1);
+ end
+
+ ySim = zeros(numel(fieldAxis), 1);
+ parfor ii = 1:nSolidAngle
+ nVers = ang2vec(ones(nTheta, 1)*phis, thetas*ones(1, nPhi));
+
+ [gax1, gInh1] = creategaxisinhbroadening(g1(ii), glw1, dgax1);
+ [gax2, gInh2] = creategaxisinhbroadening(g2(ii), glw2, dgax2);
+ ngax1 = numel(gax1);
+ ngax2 = numel(gax2);
+
+ gaussianFactor = gInh1'*gInh2;
+ gaussianFactor = gaussianFactor/sum(gaussianFactor, 'all');
+ gaussianFactor = gaussianFactor(:)';
+
+ gPlus = 1/2*(repmat(gax1', [1, ngax2]) + repmat(gax2, [ngax1, 1]));
+ gPlus = gPlus(:)';
+ gMinus = 1/2*(repmat(gax1', [1, ngax2]) - repmat(gax2, [ngax1, 1]));
+ gMinus = gMinus(:)';
+ %}
+
+ w0_ = gvalue2freq(fieldAxis', gPlus);
+ deltaw_ = gvalue2freq(fieldAxis', gMinus);
+ sintheta_ = sqrt( nVers(1, ii)^2 + nVers(2, ii)^2 );
+
+ quantNumNuc = -nNuc/2:1:nNuc/2;
+ signJmd = [-1, 1, -1, 1];
+ signOmega = [-1, -1, 1, 1];
+ signIntens = [1, -1, 1, -1];
+ if nNuc > 0
+ pascalMatrix = pascal(nHfiLine);
+ % Antidiag
+ pascalFactor = pascalMatrix(nHfiLine:nHfiLine - 1:end - 1);
+ else
+ pascalFactor = 1;
+ end
+
+ for itrans = 1:4
+ for ihfi = 1:nHfiLine
+ w0__ = w0_ + quantNumNuc(ihfi)*APlus(ii);
+ deltaw__ = deltaw_ + quantNumNuc(ihfi)*AMinus(ii);
+ Omega__ = hypot(deltaw__, JJ+dd(ii)/2);
+ wReson = w0__ + signJmd(itrans)*(JJ - dd(ii)) + ...
+ signOmega(itrans)*Omega__;
+
+ intensityReson = signIntens(itrans)*1/8*...
+ (deltaw__.^2)./(Omega__.^2);
+ trSignal__ = gaussianresonancebsweep(...
+ wReson*1e-3, mwFreq, resLw, "lwpp");
+ trSignal__ = intensityReson.*trSignal__*...
+ pascalFactor(ihfi).*gaussianFactor;
+ ySim = ySim + sum(trSignal__, 2)*sintheta_;
+ end
+ end
+ %}
+ end
+ %}
+% ySim = ySim/max(abs(ySim));
+
+end
diff --git a/zech_gBroadening_leastSquareFit.m b/zech_gBroadening_leastSquareFit.m
index 92b4f229a52d6cc621793a42d58984770ae0bd7b..9a33fbbb032fa1a58b69a9e1ae99f44a81e3fe58 100755
--- a/zech_gBroadening_leastSquareFit.m
+++ b/zech_gBroadening_leastSquareFit.m
@@ -1,7 +1,6 @@
%
clearvars
-addpath(genpath('S:\soft\matlab\'))
-addpath(genpath('/net/storage/gianlum33/soft/matlab/'))
+
addpath('S:\projects\oop_ciss_calculations\util\')
%% Easyspin simulation
diff --git a/zech_psiOopEseem_powderAverage_inhBroadening.m b/zech_psiOopEseem_powderAverage_inhBroadening.m
index 91557d4c151dc1f4da36bec86ca8d48d321cf058..f3d69ae9043dc2e076fa4e46e9c718ec1bc56bec 100755
--- a/zech_psiOopEseem_powderAverage_inhBroadening.m
+++ b/zech_psiOopEseem_powderAverage_inhBroadening.m
@@ -1,7 +1,5 @@
%
clearvars
-addpath(genpath('S:\soft\matlab\'))
-addpath(genpath('/net/storage/gianlum33/soft/matlab/'))
addpath(genpath('/net/storage/gianlum33/projects/oop_ciss_calculations/util/'))
addpath(genpath('S:\projects/oop_ciss_calculations/util/'))
@@ -15,9 +13,12 @@ fitName = 'data/digitized/zech_p46_oopEseem_fit.csv';
importedData = readtable(expName);
x = importedData{:, 1}*pi/180;
y = importedData{:, 2};
+y = rescaledata(y, 'maxabs');
+
importedData = readtable(fitName);
xx = importedData{:, 1}*pi/180;
yy = importedData{:, 2};
+yy = rescaledata(yy, 'maxabs');
crystalEseemZech = @(p, beta) -( ...
0.5*sin(2*beta) .* sin(2*p).^4 + ...
@@ -32,9 +33,9 @@ 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-')
+plot(x*180/pi, y, 'o-')
hold on
-plot(xx*180/pi, rescaledata(yy, 'maxabs'))
+plot(xx*180/pi, yy)
plot(xx*180/pi, initFit)
plot(xx*180/pi, initFit2)
plot(xx*180/pi, expectFit)
@@ -49,14 +50,6 @@ Sys0.gFrame = [-10 -128 -83; ...
0 0 0]*pi/180;
Sys0.eeFrame = [0 90 0]*pi/180; % zD direction is along -x of A1-
-% tau = 0.; % us
-
-% Assume finite bandwidth of the excitation pulse
-% pulseBw = 0.5e-3; % T
-% gPulseBw = planck*mwFreq/bmagn*(pulseBw/B0/B0);
-% gPulsePosition = planck*mwFreq/bmagn/B0;
-% excitRange = [-1/2, +1/2]*gPulseBw + gPulsePosition;
-
Sys.g = Sys0.g;
Sys.gFrame = Sys0.gFrame;
Sys.dip = mt2mhz(-0.170); % MHz
@@ -68,7 +61,7 @@ Sys.AFrame = [0 0 0 60 -90 0]*pi/180;
Sys.mwFreq = 9.7;
Sys.x = linspace(344.5, 347, 301); % mT
-Sys.turningAngleAxis = xx;
+Sys.firstPulseTurningAngles = xx;
Sys.nTheta = 90; %, 50};
Sys.nPhi = 45; %, 20};
Sys.lw1 = 0.5; % mT
@@ -78,260 +71,230 @@ Sys.nSamplegInh2 = 100;
%%
-[a0, b0] = oopEseemVsBeta(Sys, 0);
-tic
-[a1, b1] = oopEseemVsBeta(Sys, 1);
-toc
+[a0, b0] = oopeseemvsbeta(Sys, 0);
tic
-[a2, b2] = oopEseemVsBeta(Sys, 2);
+[a2, b2] = oopeseemvsbeta(Sys, 1);
toc
%%
aa0 = averageoversolidangle(a0, Sys.nTheta, Sys.nPhi, 2);
aa1 = averageoversolidangle(a1, Sys.nTheta, Sys.nPhi, 2);
-aa2 = averageoversolidangle(a2, Sys.nTheta, Sys.nPhi, 2);
+% aa2 = averageoversolidangle(a2, Sys.nTheta, Sys.nPhi, 2);
bb0 = averageoversolidangle(abs(b0), Sys.nTheta, Sys.nPhi, 2);
bb1 = averageoversolidangle(abs(b1), Sys.nTheta, Sys.nPhi, 2);
-bb2 = averageoversolidangle(abs(b2), Sys.nTheta, Sys.nPhi, 2);
+% bb2 = averageoversolidangle(abs(b2), Sys.nTheta, Sys.nPhi, 2);
% sum(abs(a1(:) - a2(:)))
% sum(abs(b1 - b2))
clf
plot(xx*180/pi, rescaledata(aa0, 'maxabs'))
hold on
plot(xx*180/pi, rescaledata(aa1, 'maxabs'), 'o')
-plot(xx*180/pi, rescaledata(aa2, 'maxabs'), 'x')
+% plot(xx*180/pi, rescaledata(aa2, 'maxabs'), 'x')
disp("xiMean0 = " + string(bb0*180/pi) + newline + ...
- "xiMean1 = " + string(bb1*180/pi) + newline + ...
- "xiMean2 = " + string(bb2*180/pi))
+ "xiMean1 = " + string(bb1*180/pi) + newline);% + ...
+% "xiMean2 = " + string(bb2*180/pi))
%%
-xplot = linspace(0, 13, 1000);
-% tiledlayout('flow', 'TileSpacing', 'compact', 'Padding', 'compact')
-% nexxtile
-
-clf
-plot(xplot, abs(mixinganglexi(xplot', Sys.J, ddd(1))), 'o')
-hold on
-plot(xplot, abs(mixinganglexi(xplot', Sys.J, ddd(2))), 'o')
-% plot(xplot, abs(mixinganglexi(xplot', Sys.J, ddd(3))), 'o')
-% plot(xplot, abs(mixinganglexi(xplot', Sys.J, ddd(4))), 'o')
-clf
-disp("nTheta = " + string(Sys.nTheta))
-for ith = 3 %1:Sys.nTheta
- plot(xplot, abs(mixinganglexi(xplot', Sys.J, ddd(ith)))*180/pi, 'o')
- hold on
- switch ith
- case 1
- color = 'k';
- case 2
- color = 'b';
- end
- for ihfi = 1:(Sys.nNuc + 1)
- switch ihfi
- case 1
- lstyle = '--';
- case 2
- lstyle = '-';
- case 3
- lstyle = '-.';
- end
- xline(abs(domeg(ith, ihfi)), ...
- 'Color', color, 'LineStyle', lstyle)
- yline(abs(mixinganglexi(domeg(ith, ihfi), Sys.J, ddd(ith)))*180/pi, ...
- 'Color', color, 'LineStyle', lstyle)
- end
- disp("xi central = " + string(xiii(ith, 2)*180/pi) + ...
- ", xi mean = " + string(mean(abs(xiii(ith, :)))*180/pi))
+wb = waitbar(0, 'Hi', 'Name', string(datetime));
+nSim = 1001;
+Sys.nTheta = 40;
+Sys.nPhi = 20;
+samplings = linspace(10, 1000, nSim);
+samplings = round(samplings);
+tic
+clear oop xi oopPowder xiPlot
+xiMean = zeros(1, nSim);
+for isim = 1:nSim
+ waitbar(isim/nSim, wb, sprintf('%d out of %d', isim, nSim))
+ Sys.nSamplegInh1 = samplings(isim);
+ Sys.nSamplegInh2 = samplings(isim);
+ [oop{isim}, xi{isim}] = oopeseemvsbeta(Sys);
+ xiMean(isim) = averageoversolidangle(...
+ abs(xi{isim}), Sys.nTheta, Sys.nPhi, 2);
+ oopPowder{isim} = averageoversolidangle(...
+ oop{isim}, Sys.nTheta, Sys.nPhi, 2);
+ xiPlot{isim} = reshape(xi{isim}, [Sys.nTheta, Sys.nPhi])*180/pi;
+% if isim == 1 % Pre-allocation
+% oop = repmat(oop, 1, nSim);
+% xi = repmat(xi, 1, nSim);
+% oopPowder = repmat(oopPowder, 1, nSim);
+% xiPlot = repmat(xiPlot, 1, nSim);
+% end
end
+toc
+delete(wb); clear('wb');
-% Sys.nNuc = 0;
-% [oopNoLw, xiNoLw] = oopEseemVsBeta(Sys);
-% xiMeanNoLw = averageoversolidangle(abs(xiNoLw), Sys.nTheta, Sys.nPhi, 2);
-% oopPowderNoLw = averageoversolidangle(oopNoLw, Sys.nTheta, Sys.nPhi, 2);
-% xiNoLwPlot = reshape(xiNoLw, [Sys.nTheta, Sys.nPhi])*180/pi;
-
%%
-Sys.lw1 = 0.01; % mT
-Sys.lw2 = 0.01; % mT
-Sys.nSamplegInh1 = 10;
-Sys.nSamplegInh2 = 10;
-Sys.nTheta = 3;
-Sys.nPhi = 1;
-Sys.nNuc = 0;
-[sigsig, xixi] = oopEseemVsBeta(Sys, true);
+clf
+plot(x*180/pi, y, 'o-')
+% plot(nThetas, xiMeana*180/pi, 'o')
+hold on
+plot(samplings, xiMean*180/pi, 'o')
+clf
+plot(xx*180/pi, oopPowder{1} - oopPowder{end})
+hold on
%%
-Sys.nTheta = 90;
-Sys.nPhi = 45;
-wb = waitbar(0, 'Hi');
-nnNuc = 2000;
-nNucs = 1:nnNuc;
-for inNuc = 1:nnNuc
- waitbar(inNuc/nnNuc, wb, sprintf('%d out of %d', inNuc, nnNuc))
- Sys.nNuc = nNucs(inNuc);
- [oopNoLwa{inNuc}, xiNoLwa{inNuc}] = oopEseemVsBeta(Sys);
- xiMeanNoLwa(inNuc) = averageoversolidangle(abs(xiNoLwa{inNuc}), Sys.nTheta, Sys.nPhi, 2);
- oopPowderNoLwa{inNuc} = averageoversolidangle(oopNoLwa{inNuc}, Sys.nTheta, Sys.nPhi, 2);
- xiNoLwPlota{inNuc} = reshape(xiNoLwa{inNuc}, [Sys.nTheta, Sys.nPhi])*180/pi;
-end
+% tau = 0.; % us
-%%
+% Assume finite bandwidth of the excitation pulse
+pulseBw = 0.50; % mT
+pulseFieldPosition = 346.1; % mean(Sys.x);
+Sys.excitRange = [-1/2, +1/2]*pulseBw + pulseFieldPosition;
+
+Sys.nTheta = 40;
+Sys.nPhi = 30;
+[oop, xixi, g1, g2, idxAngle, idxAngle1, idxAngle2] = oopeseemvsbeta(Sys, 0, 1);
-xplot = [0, nNucs(1:222)];
-yplot = [xiMeanNoLw, xiMeanNoLwa];
+figure(1)
+g1 = reshape(g1, [Sys.nTheta, Sys.nPhi]);
clf
-plot(xplot, yplot)
+imagesc(g1)
+colorbar
+colormap(viridis())
+hold on
+for ii = 1:numel(idxAngle1)
+ [xSquare, ySquare] = drawsquaresexcitrange(idxAngle1(ii), Sys.nTheta);
+ plot(xSquare, ySquare, 'k', 'LineWidth', 2)
+end
+for ii = 1:numel(idxAngle)
+ [xSquare, ySquare] = drawsquaresexcitrange(idxAngle(ii), Sys.nTheta);
+ plot(xSquare, ySquare, 'r', 'LineWidth', 2, 'LineStyle', '--')
+end
figure(2)
-tiledlayout('flow', 'TileSpacing', 'compact', 'Padding', 'compact')
-for ii = 1:9
- nexttile
- if ii == 1
- imagesc(abs(xiNoLwPlot))
- colorbar
- colormap(viridis())
- else
- imagesc(abs(xiNoLwPlota{ii}))
- colorbar
- colormap(viridis())
- end
- title(ii - 1)
+g2 = reshape(g2, [Sys.nTheta, Sys.nPhi]);
+clf
+imagesc(g2)
+colorbar
+colormap(viridis())
+hold on
+for ii = 1:numel(idxAngle2)
+ [xSquare, ySquare] = drawsquaresexcitrange(idxAngle2(ii), Sys.nTheta);
+ plot(xSquare, ySquare, 'k', 'LineWidth', 2)
+end
+for ii = 1:numel(idxAngle)
+ [xSquare, ySquare] = drawsquaresexcitrange(idxAngle(ii), Sys.nTheta);
+ plot(xSquare, ySquare, 'r', 'LineWidth', 2, 'LineStyle', '--')
end
-%%
-
-cmappa = viridis(10);
figure(3)
clf
-plot(xx*180/pi, rescaledata(oopPowderNoLw, 'maxabs'), 'o', 'DisplayName', '0')
+xixiplot = reshape(abs(xixi), [Sys.nTheta, Sys.nPhi]);
+imagesc(xixiplot)
+colorbar
+colormap(viridis())
hold on
-imap = 0;
-for ii = 101:110
- imap = imap + 1;
- if mod(ii, 2)
- plot(xx*180/pi, rescaledata(oopPowderNoLwa{ii}, 'maxabs'), ...
- 'DisplayName', string(ii), 'Color', cmappa(imap, :))
- else
- plot(xx*180/pi, rescaledata(oopPowderNoLwa{ii}, 'maxabs'), '--', ...
- 'DisplayName', string(ii), 'Color', cmappa(imap, :))
- end
+for ii = 1:numel(idxAngle)
+ [xSquare, ySquare] = drawsquaresexcitrange(idxAngle(ii), Sys.nTheta);
+ plot(xSquare, ySquare, 'r', 'LineWidth', 2, 'LineStyle', '--')
end
-legend()
%%
-% clf
-% tiledlayout('flow', 'TileSpacing', 'compact', 'Padding', 'compact')
-% nexttile
-% imagesc(abs(xiNoLwPlot))
-% colorbar
-% colormap(viridis())
-% nexttile
-% imagesc(abs(xiNoLwPlot2))
-% colorbar
-% colormap(viridis())
+xiMeanInfBw = averageoversolidanglemask(abs(xixi), Sys.nTheta, Sys.nPhi, 2, []);
+xiMeanBw = averageoversolidanglemask(abs(xixi), Sys.nTheta, Sys.nPhi, 2, idxAngle);
-% g inh broadening
-Sys.nNuc = 0;
-Sys.lw1 = 0.35; % mT
-Sys.lw2 = 0.35; % mT
-Sys.nSamplegInh1 = 10;
-Sys.nSamplegInh2 = 10;
+oopInfBw = averageoversolidanglemask(oop, Sys.nTheta, Sys.nPhi, 2, []);
+oopBw = averageoversolidanglemask(oop, Sys.nTheta, Sys.nPhi, 2, idxAngle);
+figure(4)
clf
-plot(x*180/pi, rescaledata(y, 'maxabs'), 'o-')
+plot(x*180/pi, y, '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(oopPowderNoLw, 'maxabs'), 'o')
-plot(xx*180/pi, rescaledata(oopPowderNoLw1, 'maxabs'), 'x')
-plot(xx*180/pi, rescaledata(oopPowderNoLw2, 'maxabs'), 'Marker', '+', 'LineStyle', 'none')
-plot(xx*180/pi, rescaledata(oopPowderNoLw3, 'maxabs'), 'Marker','square', 'LineStyle', 'none')
+plot(xx*180/pi, yy, '--')
+plot(xx*180/pi, rescaledata(oopInfBw, 'maxabs'))
+plot(xx*180/pi, rescaledata(oopBw, 'maxabs'))
+yline(0, 'Color', 'yellow')
+xline(90, 'Color', 'yellow')
+
+yyInterp = interp1(xx, yy, x);
+oopInterp = interp1(xx, rescaledata(oopBw, 'maxabs'), x);
+plot(x*180/pi, yyInterp, 'x')
+plot(x*180/pi, oopInterp, 'x')
+rmsdZech = sum( (y - yyInterp).^2)/numel(y);
+rmsdMe = sum( (y - oopInterp).^2)/numel(y);
%%
-xplot = nNucs(20:548);
-yplot = xiMeanNoLwa(20:548);
-model = @(p, x) p(1)*exp(-(x/p(2)).^p(3)) + p(4);
-model2 = @(p, x) p(1)*(1./x).^p(2) + p(3);
-p0 = [0.1, 20, 0.5, 0.02];
-p02 = [0.1, 0.95, 0.02];
+pulseBw = 10; % mT
+pulseFieldPosition = 346.1; % mean(Sys.x);
+Sys.excitRange = [-1/2, +1/2]*pulseBw + pulseFieldPosition;
+gExcitRange = planck*Sys.mwFreq/bmagn./Sys.excitRange*1e12;
+gExcitRange = flip(gExcitRange);
-pfit = nlinfit(xplot, yplot, model, p0);
-pfit2 = nlinfit(xplot, yplot, model2, p02);
-clf
-plot(xplot, yplot)
-hold on
-plot(xplot, model(pfit, xplot))
-yyaxis right
-plot(xplot, yplot - model(pfit, xplot))
-yline(0)
-
-%%
-
-Sys.nTheta = 2;
-Sys.nPhi = 1;
-Sys.nSamplegInh1 = 100;
-Sys.nSamplegInh2 = 100;
-Sys.nNuc = 0;
-[g1, y1, g2, y2, gmi, savexi, gfac, xiMultiplicity] = oopEseemVsBeta(Sys, 1);
+Sys.nTheta = 40;
+Sys.nPhi = 30;
+tic
+[oopNoLw, xiNoLw, ~] = oopeseemvsbeta(Sys, 0, 0);
+toc
+tic
+[oopLw, xiLw, ~] = oopeseemvsbeta(Sys, 1, 0);
+toc
+tic
+[oopBw, xiBw, idxAngle0] = oopeseemvsbeta(Sys, 0, 1);
+toc
+tic
+[oopLwBw, xiLwBw, ~] = oopeseemvsbeta(Sys, 1, 1);
+toc
+xiMeanNoLw = averageoversolidanglemask(abs(xiNoLw), Sys.nTheta, Sys.nPhi, 2, []);
+xiMeanLw = averageoversolidanglemask(abs(xiLw), Sys.nTheta, Sys.nPhi, 2, []);
+xiMeanBw = averageoversolidanglemask(abs(xiBw), Sys.nTheta, Sys.nPhi, 2, idxAngle0);
+xiMeanLwBw = averageoversolidanglemask(abs(xiLwBw), Sys.nTheta, Sys.nPhi, 2, []);
+oopPowderNoLw = averageoversolidanglemask(oopNoLw, Sys.nTheta, Sys.nPhi, 2, []);
+oopPowderLw = averageoversolidanglemask(oopBw, Sys.nTheta, Sys.nPhi, 2, idxAngle);
+oopPowderBw = averageoversolidanglemask(oopLw, Sys.nTheta, Sys.nPhi, 2, []);
+oopPowderLwBw = averageoversolidanglemask(oopLwBw, Sys.nTheta, Sys.nPhi, 2, idxAngle);
+
+xiNoLw = reshape(xiNoLw, [Sys.nTheta, Sys.nPhi]);
+xiBw = reshape(xiBw, [Sys.nTheta, Sys.nPhi]);
+xiLw = reshape(xiLw, [Sys.nTheta, Sys.nPhi]);
+xiLwBw = reshape(xiLwBw, [Sys.nTheta, Sys.nPhi]);
figure(1)
-clf
-tiledlayout('flow', 'TileSpacing', 'compact', 'Padding', 'compact')
+tiledlayout(2, 2, 'TileSpacing', 'compact', 'Padding', 'compact')
nexttile
-plot(g1, y1, 'o', g2, y2, 'o')
-nexttile
-imagesc(y1'*y2/sum(y1'*y2, 'all'))
-colorbar
-colormap(gca, viridis())
-nexttile
-imagesc(gmi)
-colorbar
-colormap(gca, crameri('vik', 'pivot', 0))
-nexttile
-xiplot = reshape(savexi, [1 + 2*Sys.nSamplegInh1, 1 + 2*Sys.nSamplegInh2]);
-imagesc(xiplot)
-colorbar
-colormap(gca, crameri('vik', 'pivot', 0))
-
-sum(abs(savexi.*gfac), 'all')
-
-[g1, y1, g2, y2, gmi, savexi, gfac, xiMultiplicity] = oopEseemVsBeta(Sys, 2);
-
-figure(2)
-clf
-tiledlayout('flow', 'TileSpacing', 'compact', 'Padding', 'compact')
-nexttile
-plot(g1, y1, 'o', g2, y2, 'o')
+imagesc(abs(xiNoLw)); colorbar; colormap(viridis())
+title('xiNoLw')
nexttile
-imagesc(y1'*y2/sum(y1'*y2, 'all'))
-colorbar
-colormap(gca, viridis())
+imagesc(abs(xiBw)); colorbar; colormap(viridis())
+title('xiBw')
nexttile
-plot(gmi)
+imagesc(abs(xiLw)); colorbar; colormap(viridis())
+title('xiLw')
nexttile
-plot(gmi, savexi, 'o')
-
-sum(abs(savexi.*xiMultiplicity))
+imagesc(abs(xiLwBw)); colorbar; colormap(viridis())
+title('xiBwLw')
%%
+
% function [gax1, gInh1, gax2, gInh2, gMinusAx, savexi, gaussianFactor, xiMultiplicity] = oopEseemVsBeta(Sys, isInhBroadening)
-function [signal, xi] = oopEseemVsBeta(Sys, isInhBroadening)
+% function [signal, xi, gax1, gax2, spinPairInRange, gaussianFactor, xiMultiplicity, xiMultiplicity2, gMinusAx] = oopeseemvsbeta(Sys, isInhBroadening, isFiniteExcitBw)
+function [signal, xi, idxAngleInExcit] = oopeseemvsbeta(Sys, isInhBroadening, isFiniteExcitBw)
arguments
Sys
isInhBroadening = 0
+ isFiniteExcitBw = 0
end
+ Sys.isInhBroadening = isInhBroadening;
+ Sys.isFiniteExcitBw = isFiniteExcitBw;
+ idxAngleInExcit = [];
+ idxAngleInExcit1 = [];
+ idxAngleInExcit2 = [];
+
+ errorMsg = checkfieldssys(Sys);
+ if ~isempty(errorMsg)
+ error(errorMsg)
+ end
+
fieldAxis = Sys.x;
Bmean = mean(fieldAxis);
- turningAngleAxis = Sys.turningAngleAxis;
+ firstPulseTurningAngles = Sys.firstPulseTurningAngles;
dip = Sys.dip;
JJ = Sys.J;
mwFreq = Sys.mwFreq;
@@ -343,6 +306,10 @@ function [signal, xi] = oopEseemVsBeta(Sys, isInhBroadening)
nNuc = 0;
end
nHfiLine = nNuc + 1;
+ if isFiniteExcitBw
+ gExcitRange = planck*mwFreq/bmagn./Sys.excitRange*1e12;
+ gExcitRange = flip(gExcitRange); % Increasing value
+ end
% Parameters for g broadening
lw1 = mt2mhz(Sys.lw1);
@@ -378,11 +345,27 @@ function [signal, xi] = oopEseemVsBeta(Sys, isInhBroadening)
AMinus = zeros(nSolidAngle, 1);
end
- signal = zeros(numel(turningAngleAxis), nSolidAngle);
+ signal = zeros(numel(firstPulseTurningAngles), nSolidAngle);
xi = zeros(1, nSolidAngle);
if ~isInhBroadening
for ii = 1:nSolidAngle
+ if isFiniteExcitBw
+ g1flag = isginexcitrange(g1(ii), gExcitRange);
+ g2flag = isginexcitrange(g2(ii), gExcitRange);
+ if g1flag
+ idxAngleInExcit1(end + 1) = ii;
+ end
+ if g2flag
+ idxAngleInExcit2(end + 1) = ii;
+ end
+% if ~g1flag || ~g2flag
+% continue
+% end
+ if g1flag && g2flag
+ idxAngleInExcit(end + 1) = ii;
+ end
+ end
gPlus = 1/2*(g1(ii) + g2(ii));
gMinus = 1/2*(g1(ii) - g2(ii));
@@ -399,7 +382,7 @@ function [signal, xi] = oopEseemVsBeta(Sys, isInhBroadening)
pascalFactor = 1;
end
- signal__ = zeros(numel(turningAngleAxis), nHfiLine);
+ signal__ = zeros(numel(firstPulseTurningAngles), nHfiLine);
xi__ = zeros(1, nHfiLine);
for ihfi = 1:nHfiLine
% w0__ = w0_ + quantNumNuc(ihfi)*APlus(ii);
@@ -409,70 +392,15 @@ function [signal, xi] = oopEseemVsBeta(Sys, isInhBroadening)
xi__(ihfi) = mixinganglexi(deltaw__, JJ, dd(ii));
% savexi(ii, ihfi) = xi__(ihfi);
signal__(:, ihfi) = ...
- crystalEseem(xi__(ihfi), turningAngleAxis);
+ crystalEseem(xi__(ihfi), firstPulseTurningAngles);
end
% Average over hfi probability
xi(ii) = xi__*pascalFactor;
signal(:, ii) = signal__*pascalFactor;
end
- elseif isInhBroadening == 1
- parfor ii = 1:nSolidAngle
- nVers = ang2vec(ones(nTheta, 1)*phis, thetas*ones(1, nPhi));
-
- [gax1, gInh1] = creategaxisinhbroadening(g1(ii), glw1, dgax1);
- [gax2, gInh2] = creategaxisinhbroadening(g2(ii), glw2, dgax2);
- ngax1 = numel(gax1);
- ngax2 = numel(gax2);
-
- gaussianFactor = gInh1'*gInh2;
- gaussianFactor = gaussianFactor/sum(gaussianFactor, 'all');
- gaussianFactor = gaussianFactor(:)';
-
- gPlus = ...
- 1/2*(repmat(gax1', [1, ngax2]) + repmat(gax2, [ngax1, 1]));
- gPlus = gPlus(:)';
- gMinus = ...
- 1/2*(repmat(gax1', [1, ngax2]) - repmat(gax2, [ngax1, 1]));
- gMinusAx = gMinus;
- gMinus = gMinus(:)';
- %}
-
- w0_ = gvalue2freq(Bmean, gPlus);
- deltaw_ = gvalue2freq(Bmean, gMinus);
- % sintheta_ = sqrt( nVers(1, ii)^2 + nVers(2, ii)^2 );
-
- quantNumNuc = -nNuc/2:1:nNuc/2;
- if nNuc > 0
- pascalMatrix = pascal(nHfiLine);
- % Antidiag
- pascalFactor = pascalMatrix(nHfiLine:nHfiLine - 1:end - 1);
- else
- pascalFactor = 1;
- end
-
- signal__ = zeros(numel(turningAngleAxis), nHfiLine, numel(gaussianFactor));
- xi__ = zeros(1, nHfiLine, numel(gaussianFactor));
- for ihfi = 1:nHfiLine
- % w0__ = w0_ + quantNumNuc(ihfi)*APlus(ii);
- deltaw__ = deltaw_ + quantNumNuc(ihfi)*AMinus(ii);
-
- xi__(1, ihfi, :) = mixinganglexi(deltaw__, JJ, dd(ii));
- savexi = xi__(ihfi, :);
- signal__(:, ihfi, :) = ...
- crystalEseem(xi__(1, ihfi, :), turningAngleAxis);
- end
- xi(ii) = squeeze(sum(...
- xi__.*reshape(gaussianFactor, [1, size(gaussianFactor)]), ...
- 3 ...
- ))*pascalFactor;
- signal(:, ii) = squeeze(sum(...
- signal__.*reshape(gaussianFactor, [1, size(gaussianFactor)]), ...
- 3 ...
- ))*pascalFactor;
- end
- elseif isInhBroadening == 2
- parfor ii = 1:nSolidAngle
+ else
+ for ii = 1:nSolidAngle
nVers = ang2vec(ones(nTheta, 1)*phis, thetas*ones(1, nPhi));
[gax1, gInh1] = creategaxisinhbroadening(g1(ii), glw1, dgax1);
@@ -492,10 +420,19 @@ function [signal, xi] = oopEseemVsBeta(Sys, isInhBroadening)
gaussianFactor = gInh1'*gInh2;
gaussianFactor = gaussianFactor/sum(gaussianFactor, 'all');
+
+ if isFiniteExcitBw
+ flag1 = isginexcitrange(gax1, gExcitRange);
+ flag2 = isginexcitrange(gax2, gExcitRange);
+ spinPairInRange = flag1'*flag2;
+ else
+ spinPairInRange = ones(ngax1, ngax2);
+ end
diagsToExtract = -(ngax1 - 1):(ngax2 - 1); % All the diagonals
xiMultiplicity = sum(spdiags(...
- gaussianFactor, diagsToExtract), 1);
-
+ gaussianFactor.*spinPairInRange, diagsToExtract), 1);
+ xiMultiplici
+
quantNumNuc = -nNuc/2:1:nNuc/2;
if nNuc > 0
pascalMatrix = pascal(nHfiLine);
@@ -505,7 +442,7 @@ function [signal, xi] = oopEseemVsBeta(Sys, isInhBroadening)
pascalFactor = 1;
end
- signal__ = zeros(numel(turningAngleAxis), nHfiLine, ngMinusAx);
+ signal__ = zeros(numel(firstPulseTurningAngles), nHfiLine, ngMinusAx);
xi__ = zeros(1, nHfiLine, ngMinusAx);
for ihfi = 1:nHfiLine
% w0__ = w0_ + quantNumNuc(ihfi)*APlus(ii);
@@ -514,7 +451,7 @@ function [signal, xi] = oopEseemVsBeta(Sys, isInhBroadening)
xi__(1, ihfi, :) = mixinganglexi(deltaw__, JJ, dd(ii));
% savexi = xi__(ihfi, :);
signal__(:, ihfi, :) = ...
- crystalEseem(xi__(1, ihfi, :), turningAngleAxis);
+ crystalEseem(xi__(1, ihfi, :), firstPulseTurningAngles);
end
xi(ii) = squeeze(sum(...
xi__.*reshape(xiMultiplicity, [1, 1, ngMinusAx]), ...
@@ -528,976 +465,55 @@ function [signal, xi] = oopEseemVsBeta(Sys, isInhBroadening)
end
end
-
function signal = crystalEseem(p, beta)
% p = xi = 2alpha
% beta = xx is the turning angle axis
signal = - (1/2 * sin(2*beta) .* cos(p).^4 + ...
1/4 * sin(beta) .* sin(2*p).^2);
end
-%{
-
-% Temporal prefactor sin(2*(dd - JJ)*tau)
-if tau ~= 0
- timePrefac = sin(2*(dd - JJ)*tau);
- timePrefac = insertDimensionPos1(timePrefac, 1);
-else
- timePrefac = ones(1, nTheta, nPhi);
-end
-signalNoLw = signalNoLw.*timePrefac;
-
-toc
-
-%%
-
-%% Optimized for loop
-
-lwArray = linspace(0.1, 1, 10);
-% for ilw = 5:9
-for ilw = 1:numel(lwArray)
-% ilw = 2;
-lw1 = mt2mhz(lwArray(ilw)); % MHz
-lw2 = mt2mhz(lwArray(ilw)); % MHz
-% The g-factor linewidth
-glw1 = w2gFunc(B0, lw1);
-glw2 = w2gFunc(B0, lw2);
-
-% Spacing between g-values
-dgax = 1e-6;
-% dgaxInterp = 5e-6;
-
-signal = zeros(numel(xx), nTheta*nPhi);
-asignal = zeros(numel(xx), nTheta*nPhi);
-xi = zeros(1, nTheta*nPhi);
-alpha = zeros(1, nTheta*nPhi);
-
-computTimes = zeros(nTheta, 1);
-
-% clear('xiPhi_', 'signalPhi_', 'alphaPhi_', 'asignalPhi_')
-tic
-parfor ii = 1:nTheta*nPhi
-% tic
-% % 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);
-
-% signalPhi_ = zeros(numel(xx), nPhi);
-% asignalPhi_ = zeros(numel(xx), nPhi);
-% xiPhi_ = zeros(1, nPhi);
-% alphaPhi_ = zeros(1, nPhi);
-% ddPhi_ = dd(ith, :);
-% g1Phi_ = g1(ith, :);
-% g2Phi_ = g2(ith, :);
-% phi_ = phis(iph);
-
- % 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
- [gax, gInh1, gInh2] = ...
- creategaxis(g1(ii), g2(ii), glw1, glw2, dgax);
- deltags = gax - gax(1);
- % "Probability matrix" of a certain spin interacting with another
- gaussianFactors = gInh1'*gInh2;
-
- % 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)
- diagsToExtract = -(numel(gInh1) - 1):(numel(gInh1) - 1);
-% xiMultiplicity = spdiags(gaussianFactors, diagsToExtract);
- xiMultiplicity = sum(spdiags(gaussianFactors, diagsToExtract), 1);
- xiMultiplicity = xiMultiplicity(end/2 + 0.5:end) + ...
- xiMultiplicity(end/2 + 0.5:-1:1);
-% xiMultiplicity(1) = xiMultiplicity(1)/2;
-
- % Interpolate
-% deltagsInterp = deltags;
-
- % [gaxInterp, ~, ~] = ...
-% creategaxis(g1(ii), g2(ii), glw1, glw2, dgaxInterp);
-% deltagsInterp = gaxInterp - gaxInterp(1);
-% xiMultiplicity(1) = xiMultiplicity(1)*2;
-% xiMultiplicityInterp = ...
-% interp1(deltags, xiMultiplicity, deltagsInterp, ...
-% 'linear', 'extrap');
-% xiMultiplicityInterp = ...
-% xiMultiplicityInterp/sum(xiMultiplicityInterp);
- % xiMultiplicity2(1) = 1/2*xiMultiplicity2(1);
-
- % Mixing angles
- xi_ = atan( (dd(ii) + 2*JJ)*1e6 * ...
- planck/(bmagn*B0*1e-3) ./ deltags);
-% planck/(bmagn*B0*1e-3) ./ deltagsInterp );
-% alpha_ = 1/2*atan( (bmagn*B0*1e-3)/planck.*(deltagsInterp)/2 ./ ...
-% ((JJ + dd(ii)/2)*1e6));
- % Average over inh. broadening
- xiWeight = xiMultiplicity/sum(xiMultiplicity); % Interp)/
- xi(ii) = sum(abs(xi_).*xiWeight);
-% alpha(ii) = sum(abs(alpha_).*xiMultiplicityInterp)/...
-% sum(xiMultiplicityInterp);
- % Signal
- signalInh = crystalEseem(xi_, xx);
-% asignalInh = crystalEseemZech(xi_, xx);
- % Average over inh. broadening
- signal(:, ii) = sum(signalInh.*xiWeight, 2);
-% asignal(:, ii) = sum(asignalInh.*xiMultiplicityInterp, 2)/ ...
-% sum(xiMultiplicityInterp);
-
-
-% ngaxLarge = numel(gaxLarge);
-% ngax = numel(gaussianFactors(1, :));
-
-% gaxIdxs = round(ngaxLarge/2 - ngax/2):round(ngaxLarge/2 + ngax/2);
-% gax = gaxLarge(gaxIdxs);
-% gInh1 = gInh1Large(gaxIdxs);
-% gInh2 = gInh2Large(gaxIdxs);
-% xi_ = atan( (ddPhi_(iph) + 2*JJ)*1e6 * ...
-% planck/(2*pi*(bmagn*B0)) ./ deltags );
-
-
-
-end
- % Store
-% signal(:, ith, :) = signalPhi_;
-% asignal(:, ith, :) = asignalPhi_;
-% xi(ith, :) = xiPhi_;
-% alpha(ith, :) = alphaPhi_;
-% computTimes(ith) = toc;
-% computTimes(ith)
-% end
-
-% Average over solid angle
-solidAngleWeight = sin(thetas)/sum(sin(thetas)*ones(1, nPhi), 'all');
-xi = reshape(xi, [nTheta, nPhi]);
-alpha = reshape(alpha, [nTheta, nPhi]);
-signal = reshape(signal, [numel(xx), nTheta, nPhi]);
-asignal = reshape(asignal, [numel(xx), nTheta, nPhi]);
-xiMean = sum(solidAngleWeight.*abs(xi), 'all');
-signalPowder = sum(squeeze(sum(solidAngleWeight'.*signal, 3)), 2);
-alphaMean = sum(solidAngleWeight.*abs(alpha), 'all');
-asignalPowder = sum(squeeze(sum(solidAngleWeight'.*asignal, 3)), 2);
-timeDuration = toc;
-disp(append(newline, ...
- 'Lwpp = ', string(lwArray(ilw)), ' mT', newline, ...
- 'xiMean = ', string(xiMean*180/pi), ' deg', newline, ...
- 'Time duration = ', string(timeDuration/60), ' min'))
-% disp(xiMean*180/pi)
-
-end
-%%
-clf
-imagesc(abs(xi)*180/pi)
-colorbar
-colormap(viridis())
-
-%%
-
-for ii = 1:12
-% disp((abs(xiNoLw(ii)) - xi(ii))*180/pi)
-end
-mean(((xi(1:12) - abs(xiNoLw(1:12))))*180/pi)
-%%
-ii = 1;
-% Consider inhomogenous broadening
- % Create g-axis in order to 'sample' the inh. broadened g-values
- [gax, gInh1, gInh2] = ...
- creategaxis(g1(ii), g2(ii), glw1, glw2, dgax);
- deltags = gax - gax(1);
- % "Probability matrix" of a certain spin interacting with another
- gaussianFactors = sparse(gInh1'*gInh2);
-
- % 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)
- diagsToExtract = -(numel(gInh1) - 1):(numel(gInh1) - 1);
- % xiMultiplicity = spdiags(gaussianFactors, diagsToExtract);
- xiMultiplicity = sum(spdiags(gaussianFactors, diagsToExtract), 1);
- xiMultiplicity = xiMultiplicity(end/2 + 0.5:end) + ...
- xiMultiplicity(end/2 + 0.5:-1:1);
- % xiMultiplicity(1) = xiMultiplicity(1)/2;
-
- % Interpolate
- [gaxInterp, ~, ~] = ...
- creategaxis(g1(ii), g2(ii), glw1, glw2, dgaxInterp);
- deltagsInterp = gaxInterp - gaxInterp(1);
- % xiMultiplicity(1) = xiMultiplicity(1)*2;
- xiMultiplicityInterp = ...
- interp1(deltags, xiMultiplicity, deltagsInterp, ...
- 'linear', 'extrap');
- xiMultiplicityInterp = ...
- xiMultiplicityInterp/sum(xiMultiplicityInterp);
- % xiMultiplicity2(1) = 1/2*xiMultiplicity2(1);
-
- % Mixing angles
- xi_ = atan( (dd(ii) + 2*JJ)*1e6 * ...
- planck/(bmagn*B0*1e-3) ./ deltagsInterp );
- alpha_ = 1/2*atan( (bmagn*B0*1e-3)/planck.*(deltagsInterp)/2 ./ ...
- ((JJ + dd(ii)/2)*1e6));
- % Average over inh. broadening
- xi(ii) = sum(abs(xi_).*xiMultiplicityInterp);
- alpha(ii) = sum(abs(alpha_).*xiMultiplicityInterp);
- % Signal
- signalInh = crystalEseem(xi_, xx);
- asignalInh = crystalEseemZech(xi_, xx);
- % Average over inh. broadening
- signal(:, ii) = sum(signalInh.*xiMultiplicityInterp, 2)/ ...
- sum(xiMultiplicityInterp);
- asignal(:, ii) = sum(asignalInh.*xiMultiplicityInterp, 2)/ ...
- sum(xiMultiplicityInterp);
-
-clf
-plot(deltags, xiMultiplicity, 'o')
-yyaxis right
-plot(deltags, xi_)
-
-% figure(1)
-% clf
-% imagesc(gaussianFactors)
-% colorbar
-% [~, idx1] = max(gaussianFactors(:, 100));
-% [~, idx2] = max(gaussianFactors(100, :));
-% xline(idx2, 'r')
-% yline(idx1, 'r')
-% figure(2)
-% clf
-% plot(gax, gaussianFactors(idx1, :)/max(gaussianFactors(idx1, :)))
-% hold on
-% plot(gax, gInh2/max(gInh2), 'o')
-%%
-% savePath = append("powderAverage_inhBroadening_glw_0", string(ilw), ".mat");
-% save(savePath, 'xi', 'signal', 'xiMean', 'signalPowder')
-
-% delete(wbarTh)
-% clear('wbarTh')
-% end
-
-% clf
-% plot(xx, signal_{ith, iph})
-% hold on
-% plot(xx, signalNoLw_{ith, iph})
-
-% max(signal_{ith, iph})
-
-%{
-%% Optimized for loop
-
-% lwArray = logspace(-1, 2, 10);
-% for ilw = 5:9
- lw1 = mt2mhz(0.35); % MHz
- lw2 = mt2mhz(0.35); % MHz
- % The g-factor linewidth
- glw1 = w2gFunc(B0, lw1);
- glw2 = w2gFunc(B0, lw2);
-
- % Spacing between g-values
- dgax = 1e-5;
- dgaxInterp = 1e-7;
-
- signal = zeros(numel(xx), nTheta, nPhi);
- asignal = zeros(numel(xx), nTheta, nPhi);
- xi = zeros(nTheta, nPhi);
- alpha = zeros(nTheta, nPhi);
-
- computTimes = zeros(nTheta, 1);
-
- clear('xiPhi_', 'signalPhi_', 'alphaPhi_', 'asignalPhi_')
- % tic
- for ith = 1:nTheta
- tic
- % 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);
-
- signalPhi_ = zeros(numel(xx), nPhi);
- asignalPhi_ = zeros(numel(xx), nPhi);
- xiPhi_ = zeros(1, nPhi);
- alphaPhi_ = zeros(1, nPhi);
- ddPhi_ = dd(ith, :);
- g1Phi_ = g1(ith, :);
- g2Phi_ = g2(ith, :);
- parfor iph = 1:nPhi
- phi_ = phis(iph);
-
- % 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
- [gax, gInh1, gInh2] = ...
- creategaxis(g1Phi_(iph), g2Phi_(iph), glw1, glw2, dgax);
- deltags = gax - gax(1);
- % "Probability matrix" of a certain spin interacting with another
- gaussianFactors = gInh1'*gInh2;
-
- % 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)
- diagsToExtract = -(numel(gInh1) - 1):(numel(gInh1) - 1);
- % xiMultiplicity = spdiags(gaussianFactors, diagsToExtract);
- xiMultiplicity = sum(spdiags(gaussianFactors, diagsToExtract), 1);
- xiMultiplicity = xiMultiplicity(end/2 + 0.5:end) + ...
- xiMultiplicity(end/2 + 0.5:-1:1);
- % xiMultiplicity(1) = xiMultiplicity(1)/2;
-
- % Interpolate
- [gaxInterp, ~, ~] = ...
- creategaxis(g1Phi_(iph), g2Phi_(iph), glw1, glw2, dgaxInterp);
- deltagsInterp = gaxInterp - gaxInterp(1);
- % xiMultiplicity(1) = xiMultiplicity(1)*2;
- xiMultiplicityInterp = ...
- interp1(deltags, xiMultiplicity, deltagsInterp, ...
- 'linear', 'extrap');
- xiMultiplicityInterp = ...
- xiMultiplicityInterp/sum(xiMultiplicityInterp);
- % xiMultiplicity2(1) = 1/2*xiMultiplicity2(1);
-
- % Mixing angles
- xi_ = atan( (ddPhi_(iph) + 2*JJ)*1e6 * ...
- planck/(bmagn*B0) ./ deltagsInterp );
- alpha_ = 1/2*atan( (bmagn*B0)/planck.*(deltagsInterp)/2 ./ ...
- ((JJ + ddPhi_(iph)/2)*1e6));
- % Average over inh. broadening
- xiPhi_(iph) = sum(abs(xi_).*xiMultiplicityInterp)/...
- sum(xiMultiplicityInterp);
- alphaPhi_(iph) = sum(abs(alpha_).*xiMultiplicityInterp)/...
- sum(xiMultiplicityInterp);
- % Signal
- signalInh = crystalEseem(xi_, xx);
- asignalInh = crystalEseemZech(xi_, xx);
- % Average over inh. broadening
- signalPhi_(:, iph) = sum(signalInh.*xiMultiplicityInterp, 2)/ ...
- sum(xiMultiplicityInterp);
- asignalPhi_(:, iph) = sum(asignalInh.*xiMultiplicityInterp, 2)/ ...
- sum(xiMultiplicityInterp);
-
-
- % ngaxLarge = numel(gaxLarge);
- % ngax = numel(gaussianFactors(1, :));
-
- % gaxIdxs = round(ngaxLarge/2 - ngax/2):round(ngaxLarge/2 + ngax/2);
- % gax = gaxLarge(gaxIdxs);
- % gInh1 = gInh1Large(gaxIdxs);
- % gInh2 = gInh2Large(gaxIdxs);
- % xi_ = atan( (ddPhi_(iph) + 2*JJ)*1e6 * ...
- % planck/(2*pi*(bmagn*B0)) ./ deltags );
-
-
+function errorMsg = checkfieldssys(Sys)
+
+% excitRange
+errorMsgContent = "excitRange should be a vector of " + ...
+ "two elements of increasing value.";
+errorMsg = '';
+if Sys.isFiniteExcitBw
+ if isempty(Sys.excitRange)
+ errorMsg = errorMsgContent;
+ elseif ~isa(Sys.excitRange, 'double')
+ errorMsg = errorMsgContent;
+ else
+ if size(Sys.excitRange) ~= 2
+ errorMsg = errorMsgContent;
+ elseif Sys.excitRange(1) >= Sys.excitRange(2)
+ errorMsg = errorMsgContent;
end
- % Store
- signal(:, ith, :) = signalPhi_;
- asignal(:, ith, :) = asignalPhi_;
- xi(ith, :) = xiPhi_;
- alpha(ith, :) = alphaPhi_;
- computTimes(ith) = toc;
-% computTimes(ith)
end
- % toc
-
- % Average over solid angle
- xiMean = sum(sin(thetas').*abs(xi), 'all')/ ...
- sum(sin(thetas')*ones(1, nPhi), 'all');
- signalPowder = sum(squeeze(sum(sin(thetas).*signal, 3)), 2)/ ...
- sum(sin(thetas')*ones(1, nPhi), 'all');
- alphaMean = sum(sin(thetas').*abs(alpha), 'all')/ ...
- sum(sin(thetas')*ones(1, nPhi), 'all');
- asignalPowder = sum(squeeze(sum(sin(thetas).*asignal, 3)), 2)/ ...
- sum(sin(thetas')*ones(1, nPhi), 'all');
-
-% savePath = append("powderAverage_inhBroadening_glw_0", string(ilw), ".mat");
-% save(savePath, 'xi', 'signal', 'xiMean', 'signalPowder')
-
- delete(wbarTh)
- clear('wbarTh')
-% end
-
-% clf
-% plot(xx, signal_{ith, iph})
-% hold on
-% plot(xx, signalNoLw_{ith, iph})
-
-% max(signal_{ith, iph})
-%}
-%%
-
-%{
-%% Check if it is the same to use xi or alpha
-
-dip = unitconvert(-0.177,'mT->MHz'); % MHz, -4.9 MHz
-JJ = unitconvert(1e-3,'mT->MHz'); % MHz
-
-% Expected input: B0 in mT
-% Output: frequency nu in GHz
-g2wFunc = @(B0, g) bmagn*B0/planck*g*1e-12;
-% Expected input: deltaw in GHz, J in MHz, d in MHz
-% Output: Omega in GHz
-OmegaFunc = @(deltaw, J, d) sqrt(deltaw.^2 + (J*1e-3 + d/2*1e-3).^2);
-% Dipolar interaction
-% Makes use of the fact that cos(thetaD) = dotProduct(B0, zD)
-% Expected input:
-% dip: 1 x 1
-% nVers: 3 x nTheta x nPhi
-% zD: 3 x 1
-dipFunc = @(dip, nVers, zD) dip*((sum(nVers.*zD)).^2 - 1/3);
-% Input vec of dimensions a1 x a2 x a3 x ... and return vec of dimensions
-% n x a1 x a2 x a3 x ... with vec(i, :, :, ...) = vec(j, :, :, ...) for
-% every i and j
-insertDimensionPos1 = @(vec, n) ...
- repmat(reshape(vec, [1, size(vec)]), [n, ones(1, numel(size(vec)))]);
-
-% Mixing angle, p = [dd, JJ, deltaw] in [MHz, MHz, GHz]
-xiFunc = @(dd, JJ, deltaw) atan( (dd/2 + JJ).*1e-3 ./ deltaw);
-alphaFunc = @(dd, JJ, deltaw) 1/2*atan( deltaw ./ ((JJ + dd/2)*1e-3));
-
-tau = 0.; % us
-
-% Assume finite bandwidth of the excitation pulse
-pulseBw = 0.5e-3; % T
-gPulseBw = planck*mwFreq/bmagn*(pulseBw/B0/B0);
-gPulsePosition = planck*mwFreq/bmagn/B0;
-excitRange = [-1/2, +1/2]*gPulseBw + gPulsePosition;
-
-% Theta and phi grid
-thetas = (0:3:180)*pi/180;
-nTheta = numel(thetas);
-phis = (0:3:360)*pi/180;
-nPhi = numel(phis);
-
-tic
-
-clear('dd', 'nVers', 'nVersInFrame1', 'g1', 'g2', 'xiNoLw', 'xiNoLw3D', ...
- 'alphaNoLw', 'alphaNoLw3D', ....
- 'signalNoLw', 'signalPowderNoLw', 'asignalNoLw', 'asignalPowderNoLw')
-
-% Direction of B0
-clear('nVers')
-nVers(1, :, :) = sin(thetas')*cos(phis);
-nVers(2, :, :) = sin(thetas')*sin(phis);
-nVers(3, :, :) = cos(thetas')*ones(1, nPhi);
-% Effective g-values
-g2 = squeeze(sqrt( sum( (Sys.g(2, :)'.*nVers).^2)));
-euMatrixg1 = erot(Sys.gFrame(1, :));
-g1TensorInFrame2 = euMatrixg1'*diag(Sys.g(1, :))*euMatrixg1;
-g1 = squeeze(sqrt( sum( (pagemtimes(g1TensorInFrame2, nVers)).^2)));
-% Dipolar interaction
-zD = erot(Sys.eeFrame)*[0, 0, 1]';
-dd = squeeze(dipFunc(dip, nVers, zD));
-
-deltaw = g2wFunc(B0, (g1 - g2)/2); % GHz
-
-%
-% Mixing angles
-%
-alphaNoLw = alphaFunc(dd, JJ, deltaw);
-xiNoLw = xiFunc(dd, JJ, deltaw);
-% Average over all the solid angles with proper normalization
-solidAngleWeight = sin(thetas')/sum(sin(thetas')*ones(1, nPhi), 'all');
-alphaMeanNoLw = sum(solidAngleWeight.*abs(alphaNoLw), 'all');
-xiMeanNoLw = sum(solidAngleWeight.*abs(xiNoLw), 'all');
-
-%
-% Signal
-%
-asignalNoLw = crystalEseemZech(insertDimensionPos1(alphaNoLw, 1), xx);
-signalNoLw = crystalEseem(insertDimensionPos1(xiNoLw, 1), xx);
-
-% Temporal prefactor sin(2*(dd - JJ)*tau)
-if tau ~= 0
- timePrefac = sin(2*(dd - JJ)*tau);
- timePrefac = insertDimensionPos1(timePrefac, 1);
-else
- timePrefac = ones(1, nTheta, nPhi);
end
-asignalNoLw = asignalNoLw.*timePrefac;
-signalNoLw = signalNoLw.*timePrefac;
-
-% Sum every angular contribution to the signal
-asignalPowderNoLw = sum(squeeze(sum(solidAngleWeight'.*asignalNoLw, 3)), 2);
-signalPowderNoLw = sum(squeeze(sum(solidAngleWeight'.*signalNoLw, 3)), 2);
-
-toc
-%}
-
-%%
-
-% Fit to determine the best mixing angle if it were a crystal
-fitModel = @(p) p(1)*crystalEseem(p(2), xx);
-p0 = [2, 6.9*pi/180];
-vary = [2, 40*pi/180];
-FitOpt.x = xx;
-
-ydata = rescaledata(signalPowderNoLw, 'maxabs');
-Fit = esfit(ydata, fitModel, p0, vary, FitOpt);
-
-clf
-plot(xx, ydata, xx, Fit.fit);
-yline(0, '--')
-yyaxis right
-plot(xx, Fit.residuals)
-yline(0, '--')
-
-Fit.pfit(2)*180/pi % Best fit crystal mixing angle
-
-%{
-%%
-
-
-figure(6)
-mySignal = rescaledata(signalPowderNoLw, 'maxabs');
-testManualFit = rescaledata(crystalEseem(6.18*pi/180, xx), 'maxabs');
-crystalMean = rescaledata(crystalEseem(xiNoLwMean, xx), 'maxabs');
-
-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, 'r--', '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')
-
-% [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
-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()
-
-%}
-
-
-%%
-
-lwArrayTest = linspace(0.1, 50, 100);
-ith = 15;
-iph = 15;
-% ith = 342;
-% iph = 357;
-% for ilw = 1:100
-% lw1 = lwArrayTest(ilw)*1e6; % Hz
-% lw2 = lwArrayTest(ilw)*1e6; % Hz
- lw1 = 5e6;
- lw2 = 5e6;
- glw1 = planck/bmagn/B0*lw1;
- glw2 = planck/bmagn/B0*lw2;
- signalPhi_ = zeros(numel(xx), nPhi);
- dgax = 1e-5;
- dgaxInterp = 1e-7;
- xiPhi_ = zeros(1, nPhi);
- ddPhi_ = dd(ith, :);
- g1Phi_ = g1(ith, :);
- g2Phi_ = g2(ith, :);
- [gax, gInh1, gInh2] = ...
- creategaxis(g1Phi_(iph), g2Phi_(iph), glw1, glw2, dgax);
- deltags = gax - gax(1);
- % "Probability matrix" of a certain spin interacting with another
- gaussianFactors = gInh1'*gInh2;
-
- % 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)
- diagsToExtract = -(numel(gInh1) - 1):(numel(gInh1) - 1);
- % xiMultiplicity = spdiags(gaussianFactors, diagsToExtract);
- xiMultiplicity = sum(spdiags(gaussianFactors, diagsToExtract), 1);
- xiMultiplicity = xiMultiplicity(end/2 + 0.5:end) + ...
- xiMultiplicity(end/2 + 0.5:-1:1);
-% xiMultiplicity(1) = xiMultiplicity(1)/2;
-
- % Interpolate
- [gaxInterp, ~, ~] = ...
- creategaxis(g1Phi_(iph), g2Phi_(iph), glw1, glw2, dgaxInterp);
- deltagsInterp = gaxInterp - gaxInterp(1);
-% xiMultiplicity(1) = xiMultiplicity(1)*2;
- xiMultiplicityInterp = ...
- interp1(deltags, xiMultiplicity, deltagsInterp, ...
- 'linear', 'extrap');
- % xiMultiplicity2(1) = 1/2*xiMultiplicity2(1);
-
- % Mixing angles
- xi_ = atan( (ddPhi_(iph) + 2*JJ)*1e6 * ...
- planck/(bmagn*B0) ./ deltagsInterp );
- % Average over inh. broadening
- xiPhi_(iph) = sum(abs(xi_).*xiMultiplicityInterp)/...
- sum(xiMultiplicityInterp);
- mixingAngle(ilw) = xiPhi_(iph)*180/pi;
- [valMax, idxMax] = max(xiMultiplicityInterp);
- positionMax(ilw) = deltagsInterp(idxMax);
- [~, idxFwhm] = min(abs(valMax/2 - xiMultiplicityInterp));
- xiMultFwhm(ilw) = abs(deltagsInterp(idxFwhm) - positionMax(ilw));
- % Signal
- signalInh = crystalEseem(xi_, xx);
- % Average over inh. broadening
- signalPhi_(:, iph) = sum(signalInh.*xiMultiplicityInterp, 2)/ ...
- sum(xiMultiplicityInterp);
-
-% end
-clf
-plot(deltags, gInh1, 'o-', deltags, gInh2, 'o-')
-plot(deltags, xiMultiplicity, 'o', deltagsInterp, xiMultiplicityInterp)
-yyaxis right
-plot(deltagsInterp, xi_, 'o')
-
-%%
-%{
-%% Test
-
-% tempInhMultiplicity = sum(spdiags(gaussianFactors), 1);
-% InhMultiplicity = aa1(end/2 + 0.5:end) + aa1(end/2 + 0.5:-1:1);
-clf
-ngax = numel(gax);
-plot(1:2*ngax - 1, xiMultiplicity, 1:ngax, xiMultiplicity)
-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;
-
-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);
-% signal_{ith, iph} = sum(signalSingleSumg2, 2)*sin(theta_);
-% signalPowder = signalPowder + signal_{ith, iph};
- end
- theta_ = thetas(ith);
- signalPowder = signalPowder + signal_{ith, iph}*sin(theta_);
end
-delete(wbarTh)
-% delete(wbarPh)
-%}
-
-% signalPowder = zeros(numel(xx), 1);
-% for ith = 1:nTheta
-% theta_ = thetas(ith);
-% for iph = 1:nPhi
-% signalPowder = signalPowder + signal_{ith, iph}*sin(theta_);
-% end
+function flag = isginexcitrange(g, gExcitRange)
+flag = ones(size(g));
+idxOutRange = g < gExcitRange(1) | g > gExcitRange(2);
+flag(idxOutRange) = 0;
+% if g < gExcitRange(1) || g > gExcitRange(2)
+% flag = 0; % g outside the range
+% else
+% flag = 1;
% end
-
-
-% 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);
-% xitest(ith, :) = repmat(cos(theta_)^4, 1, nPhi);
-% xitest(ith, :) = repmat(1, 1, nPhi);
-% end
-for ith = 1:nTheta
- theta_ = thetas(ith);
- xiMean = xiMean + sum(xi(ith, :))*sin(theta_);
+function [xSquare, ySquare] = drawsquaresexcitrange(idxAngle, nTheta)
+xSquare0 = floor(idxAngle/nTheta) + 0.5;
+if rem(idxAngle, nTheta) == 0
+ xSquare0 = xSquare0 - 1;
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 [gax, gInh1, gInh2] = creategaxis(g1, g2, glw1, glw2, dgax)
- % Here glw1 and glw2 are the lwpp, therefore the variance will be
- % lwpp*2
- xxFactor = 1;
- xmin = min([g1 - xxFactor*glw1, g2 - xxFactor*glw2]);
- xmax = max([g1 + xxFactor*glw1, g2 + xxFactor*glw2]);
- gax = xmin:dgax:xmax;
-
- % Inhomogeneously broadened gaussian distributions
- gInh1 = gaussian(gax, g1, glw1);
- gInh2 = gaussian(gax, g2, glw2);
- % Normalization
-% gaxForNormaliz1 = ...
-% (gax(round(end/2)) - 5*glw1):dgax:(gax(round(end/2)) + 5*glw1);
-% gaxForNormaliz2 = ...
-% (gax(round(end/2)) - 5*glw2):dgax:(gax(round(end/2)) + 5*glw2);
-% normFactorgInh1 = sum(gaussian(gaxForNormaliz1, g1, glw1));
-% normFactorgInh2 = sum(gaussian(gaxForNormaliz2, g2, glw2));
-% gInh1 = gInh1/normFactorgInh1;
-% gInh2 = gInh2/normFactorgInh2;
+ySquare0 = mod(idxAngle, nTheta) - 0.5;
+if ySquare0 == -0.5
+ ySquare0 = ySquare0 + nTheta;
end
-
-% function centeredMatrix = createcenteredmatrix(gInh1, gInh2)
-% nPoints = numel(gInh1);
-% multipl = gInh1'*gInh2;
-% [~, idxMax1] = max(multipl(:, round(nPoints)));
-% [~, idxMax2] = max(multipl(round(nPoints), :));
-% halfWidthFin = min([idxMax1, idxMax2]) - 2;
-%
-% centeredMatrix = ...
-% multipl(idxMax1 - halfWidthFin:idxMax1 + halfWidthFin, ...
-% idxMax2 - halfWidthFin:idxMax2 + halfWidthFin);
-% end
-
-
-
-%}
-
-
-
-
-
-
-
-
+xSquare = [xSquare0, xSquare0, xSquare0 + 1, xSquare0 + 1, xSquare0];
+ySquare = [ySquare0, ySquare0 + 1, ySquare0 + 1, ySquare0, ySquare0];
+end
\ No newline at end of file
diff --git a/zech_psiOopEseem_trEPRstickSpectra.m b/zech_psiOopEseem_trEPRstickSpectra.m
index a768171692d4f9a473465a37e95e23874e17530e..0e70ae8556d44187cb787b01ad2dc1c67284efb7 100755
--- a/zech_psiOopEseem_trEPRstickSpectra.m
+++ b/zech_psiOopEseem_trEPRstickSpectra.m
@@ -69,9 +69,9 @@ FitOpt.maxTime = 0.5;
%
% % ySimNoHfi = pepper(Sys, Exp);
-ySim = pepper(Sys, Exp, Opt);
+% ySim = pepper(Sys, Exp, Opt);
% yySim = ySim;
-ySim = ySim'/max(abs(ySim));
+% ySim = ySim'/max(abs(ySim));
% Ham = ham(Sys, [0, 0, B0]);
% [Ham0, mux, muy, muz] = ham(Sys);