91 lines
3.3 KiB
Matlab
91 lines
3.3 KiB
Matlab
function [AveragePCD] = extractAveragePeakColumnDensity(folder_path, run_index, radius, minPeakHeight, SuppressPlotFlag)
|
|
set(0,'defaulttextInterpreter','latex')
|
|
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
|
|
|
|
format long
|
|
|
|
% Load data
|
|
filePath = sprintf(horzcat(folder_path, '/Run_%03i/psi_gs.mat'), run_index);
|
|
|
|
try
|
|
Data = load(filePath, 'psi', 'Params', 'Transf', 'Observ');
|
|
catch ME
|
|
warning('Failed to load file: %s\n%s', filePath, ME.message);
|
|
AveragePCD = NaN;
|
|
return;
|
|
end
|
|
|
|
Params = Data.Params;
|
|
Transf = Data.Transf;
|
|
Observ = Data.Observ;
|
|
|
|
if isgpuarray(Data.psi)
|
|
psi = gather(Data.psi);
|
|
else
|
|
psi = Data.psi;
|
|
end
|
|
if isgpuarray(Data.Observ.residual)
|
|
Observ.residual = gather(Data.Observ.residual);
|
|
else
|
|
Observ.residual = Data.Observ.residual;
|
|
end
|
|
|
|
% Axes scaling and coordinates in micrometers
|
|
x = Transf.x * Params.l0 * 1e6;
|
|
y = Transf.y * Params.l0 * 1e6;
|
|
z = Transf.z * Params.l0 * 1e6;
|
|
|
|
dz = z(2)-z(1);
|
|
|
|
% Compute probability density |psi|^2
|
|
n = abs(psi).^2;
|
|
nxy = squeeze(trapz(n*dz,3));
|
|
|
|
state = nxy';
|
|
peaks = extractPeaks(state, radius, minPeakHeight);
|
|
peakHeights = state(peaks);
|
|
AveragePCD = mean(peakHeights);
|
|
[row, col] = find(peaks);
|
|
|
|
if ~SuppressPlotFlag
|
|
figure(1);
|
|
clf
|
|
set(gcf,'Position', [100, 100, 900, 900])
|
|
plotxy = pcolor(x,y,state);
|
|
set(plotxy, 'EdgeColor', 'none');
|
|
hold on;
|
|
plot(x(col), y(row), 'w+', 'MarkerSize', 8, 'LineWidth', 1.5);
|
|
cbar1 = colorbar;
|
|
cbar1.Label.Interpreter = 'latex';
|
|
% cbar1.Ticks = []; % Disable the ticks
|
|
colormap(gca, Helper.Colormaps.plasma())
|
|
xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16)
|
|
ylabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16)
|
|
title(['<CD> = ' num2str(AveragePCD)], ...
|
|
'Interpreter', 'tex', 'FontSize', 16)
|
|
set(gca, 'FontSize', 16); % For tick labels only
|
|
end
|
|
end
|
|
|
|
function [peaks_mask] = extractPeaks(image, radius, minPeakHeight)
|
|
peaks = imregionalmax(image);
|
|
peaks(image < minPeakHeight) = 0;
|
|
|
|
[row, col] = find(peaks); % row and col are the coordinates of the detected peaks
|
|
unique_peaks = true(size(row)); % Assume all peaks are unique initially
|
|
|
|
% Check distances between each peak and remove duplicates
|
|
for i = 1:length(row)
|
|
if unique_peaks(i) % If this peak hasn't been excluded
|
|
% Find the distance from this peak to all other peaks
|
|
dist = sqrt((row(i) - row).^2 + (col(i) - col).^2); % Euclidean distance
|
|
|
|
% Exclude any peak within the radius
|
|
unique_peaks(dist < radius & dist > 0) = false; % Exclude other peaks in radius
|
|
end
|
|
end
|
|
|
|
% Save the peaks
|
|
peaks_mask = false(size(image));
|
|
peaks_mask(sub2ind(size(image), row(unique_peaks), col(unique_peaks))) = true;
|
|
end |