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 ------------------ %

    % Apply 2D Fourier Transform
    fftResult    = fft2(double(nxyScaled));
    fftMagnitude = abs(fftshift(fftResult));  % Shift zero frequency to center

    % Get the size of the image
    [Ny, Nx]     = size(nxyScaled);

    % Set the DC component to zero (ignore it)
    fftMagnitude(Ny/2 + 1, Nx/2 + 1) = 0;  % Assuming fftshift was used

    % Calculate the sampling intervals (real-space distance)
    dx = mean(diff(x));  % Sampling interval in the X direction
    dy = mean(diff(y));  % Sampling interval in the Y direction
    
    % Calculate wavenumber increment (frequency axes)
    dkx = 1 / (Nx * dx);  % Increment in the X direction (in micrometers^-1)
    dky = 1 / (Ny * dy);  % Increment in the Y direction (in micrometers^-1)
    
    % Create the wavenumber axes
    kx = (-Nx/2:Nx/2-1) * dkx;  % Wavenumber axis in X (micrometers^-1)
    ky = (-Ny/2:Ny/2-1) * dky;  % Wavenumber axis in Y (micrometers^-1)
    
    [G1, G2, reciprocalAngle, periodX, periodY] = extractLatticeProperties(fftMagnitude', kx, ky);
    G1      = G1 * Params.l0;
    G2      = G2 * Params.l0;
    periodX = periodX/Params.l0;
    periodY = periodY/Params.l0;
    
    %------------------ 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, 'parula')
    % 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
    nexttile;  % Equivalent to subplot('Position', [0.05, 0.55, 0.28, 0.4]);
    plotxy    = pcolor(kx*Params.l0,ky*Params.l0,fftMagnitude');
    set(plotxy, 'EdgeColor', 'none');
    cbar1 = colorbar;
    cbar1.Label.Interpreter = 'latex';
    colormap(gca, 'parula')
    

    % Overlay the reciprocal lattice vectors
    hold on;
    quiver(0, 0, G1(1), G1(2), 0, 'r', 'LineWidth', 2, 'MaxHeadSize', 1);  % G1 vector
    quiver(0, 0, G2(1), G2(2), 0, 'b', 'LineWidth', 2, 'MaxHeadSize', 1);  % G2 vector
    
    % Display the angle between the vectors
    angleDeg = rad2deg(reciprocalAngle);  % Convert the angle to degrees
    midX     = (G1(1) + G2(1)) / 2;
    midY     = (G1(2) + G2(2)) / 2;
    
    % Add the angle text annotation to the plot
    text(midX, midY, sprintf('Angle = %.2f deg', angleDeg), 'Color', 'white', 'FontSize', 16, 'FontWeight', 'bold');
    
    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);

end

function [G1, G2, reciprocalAngle, periodX, periodY] = extractLatticeProperties(fftMagnitude, kx, ky)
    
    % Define the index ranges for the positive side (first quadrant)
    posX = kx > 0;  % Positive kx
    posY = ky > 0;  % Positive ky
    
    % Restrict the FFT magnitude to only the positive quadrant (kx > 0, ky > 0)
    fftMagnitudePos = fftMagnitude(posY, posX);
    
    % Find peaks in the Fourier domain using findpeaks
    [peaks, peakIdx] = findpeaks(fftMagnitudePos(:));  % Find all peaks in the positive quadrant
    [~, sortIdx] = sort(peaks, 'descend'); % Sort peaks in descending order
    
    % Select the two largest distinct peaks
    peakIdx1 = peakIdx(sortIdx(1));  % Index of the largest peak
    peakIdx2 = peakIdx(sortIdx(2));  % Index of the second largest peak
    
    % Convert peak indices to subscripts in the original FFT size
    [peakY1, peakX1] = ind2sub(size(fftMagnitudePos), peakIdx1);
    [peakY2, peakX2] = ind2sub(size(fftMagnitudePos), peakIdx2);
    
    % Convert these subscripts back to the full FFT index space
    peakY1 = find(posY, 1, 'first') + peakY1 - 1;
    peakX1 = find(posX, 1, 'first') + peakX1 - 1;
    
    peakY2 = find(posY, 1, 'first') + peakY2 - 1;
    peakX2 = find(posX, 1, 'first') + peakX2 - 1;

    % Extract the wavenumbers corresponding to the peaks
    kx1 = kx(peakX1); ky1 = ky(peakY1);
    kx2 = kx(peakX2); ky2 = ky(peakY2);
    
    % Reciprocal lattice vectors
    G1 = [kx1, ky1];  % Reciprocal lattice vector 1
    G2 = [kx2, ky2];  % Reciprocal lattice vector 2

    % Calculate the angle between the reciprocal lattice vectors using dot product
    dotProduct = dot(G1, G2);  % Dot product of G1 and G2
    magnitudeG1 = norm(G1);  % Magnitude of G1
    magnitudeG2 = norm(G2);  % Magnitude of G2
    
    % Angle between the reciprocal lattice vectors (in radians)
    reciprocalAngle = acos(dotProduct / (magnitudeG1 * magnitudeG2));
    
    % Correct periodicity calculations (swap kx and ky for periodicity)
    periodX = 1 / abs(kx1);  % Periodicity along x-axis (use kx1)
    periodY = 1 / abs(ky1);  % Periodicity along y-axis (use ky1)
    
end