95 lines
3.9 KiB
Matlab
95 lines
3.9 KiB
Matlab
function [contrast, periodX, periodY] = analyzeGSWavefunction(folder_path, run_index)
|
|
|
|
set(0,'defaulttextInterpreter','latex')
|
|
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
|
|
|
|
% Load data
|
|
Data = load(sprintf(horzcat(folder_path, '/Run_%03i/psi_gs.mat'),run_index),'psi','Transf','Observ','Params','VParams');
|
|
|
|
Transf = Data.Transf;
|
|
Observ = Data.Observ;
|
|
Params = Data.Params;
|
|
VParams = Data.VParams;
|
|
|
|
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
|
|
|
|
% Set long format for output
|
|
format long
|
|
|
|
% Coordinates in micrometers
|
|
x = Transf.x * Params.l0;
|
|
y = Transf.y * Params.l0;
|
|
|
|
% Compute probability density |psi|^2
|
|
nxy = abs(psi).^2;
|
|
nxyScaled = nxy*(Params.add*10^6)^2;
|
|
|
|
%------------------ Calculate contrast ------------------ %
|
|
|
|
% Find unique maximum intensity values
|
|
uniqueMaxValues = unique(nxyScaled(nxyScaled == max(nxyScaled(:))));
|
|
if length(uniqueMaxValues) > 1
|
|
maxIntensity = median(uniqueMaxValues); % Choose the median of max values
|
|
else
|
|
maxIntensity = uniqueMaxValues; % If only one, take the unique value
|
|
end
|
|
|
|
% Find unique minimum intensity values
|
|
uniqueMinValues = unique(nxyScaled(nxyScaled == min(nxyScaled(:))));
|
|
if length(uniqueMinValues) > 1
|
|
minIntensity = median(uniqueMinValues); % Choose the median of min values
|
|
else
|
|
minIntensity = uniqueMinValues; % If only one, take the unique value
|
|
end
|
|
|
|
contrast = (maxIntensity - minIntensity) / (maxIntensity + minIntensity);
|
|
|
|
%------------------ Lattice Properties ------------------ %
|
|
|
|
[kx, ky, fftMagnitude, lattice_type, periodX, periodY, freq_x, freq_y, rotation_angle] = Scripts.extractLatticeProperties(nxyScaled, x, y);
|
|
|
|
%------------------ Plotting ------------------ %
|
|
|
|
figure('Position', [100, 100, 1600, 800]);
|
|
clf
|
|
t = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x3 grid
|
|
|
|
% Plot |psi(x,y)|^2 (density in x and y plane)
|
|
nexttile; % Equivalent to subplot('Position', [0.05, 0.55, 0.28, 0.4]);
|
|
plotxy = pcolor(x/Params.l0,y/Params.l0,nxyScaled');
|
|
set(plotxy, 'EdgeColor', 'none');
|
|
cbar1 = colorbar;
|
|
cbar1.Label.Interpreter = 'latex';
|
|
colormap(gca, Helper.Colormaps.plasma())
|
|
% clim(ax1,[0.00,0.3]);
|
|
ylabel(cbar1,'$na_{dd}^2$','FontSize',16,'Rotation',270)
|
|
xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
|
|
ylabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
|
|
title(['$|\Psi(x,y)|^2$ - Contrast: ', num2str(contrast, '%.3f'), '; Period X = ', num2str(periodX, '%.2f'), '$ \mu$m', '; Period Y = ', num2str(periodY, '%.2f'), '$ \mu$m'], 'Interpreter', 'latex', 'FontSize', 14)
|
|
|
|
% Plot 2-D FFT with detected peaks
|
|
nexttile; % Equivalent to subplot('Position', [0.05, 0.55, 0.28, 0.4]);
|
|
imagesc(kx*Params.l0,ky*Params.l0, log(1 + fftMagnitude));
|
|
cbar1 = colorbar;
|
|
cbar1.Label.Interpreter = 'latex';
|
|
colormap(gca, 'jet')
|
|
hold on;
|
|
plot(freq_x*Params.l0, freq_y*Params.l0, 'ro', 'MarkerSize', 10, 'LineWidth', 2);
|
|
hold off;
|
|
xlabel('$k_x l_o$', 'Interpreter', 'latex', 'FontSize', 14)
|
|
ylabel('$k_y l_o$', 'Interpreter', 'latex', 'FontSize', 14)
|
|
title('$\mathcal{F}\{|\Psi(x,y)|^2\}$', 'Interpreter', 'latex', 'FontSize', 16);
|
|
|
|
sgtitle(['$\omega_z = 2 \pi \times$', num2str(round(Params.wz/(2*pi)), '%.2f'), ' Hz; $\theta = $', num2str(rad2deg(Params.theta), '%.2f'), '$^\circ$; ', sprintf('Detected Lattice Type: %s\n', lattice_type)], 'Interpreter', 'latex', 'FontSize', 18, 'FontWeight', 'bold')
|
|
|
|
end |