function [ModulationFlag] = determineDensityModulation(psi, Params, Transf)
    % Axes scaling and coordinates in micrometers
    x               = Transf.x * Params.l0 * 1e6; 
    dx              = x(2)-x(1); 
    % Compute probability density |psi|^2
    n               = abs(psi).^2;
    nyz             = squeeze(trapz(n*dx,1));
    
    densityProfile  = sum(nyz,2);
    
    % We need to first smoothen the density profile and subtract the smoothened profile from the
    % original density profile to - 
    % 1. Remove the dominant low-frequency content.
    % 2. Isolate the high-frequency components (like a sinusoidal modulation).
    % FFT of the residual then highlights the periodic part, making the
    % modulation easy to detect.

    % Step 1: Smooth the density profile (Gaussian smoothing)
    smoothedProfile = smooth(densityProfile, 10);
    
    % Step 2: Compute the residual (original - smoothed)
    residual        = densityProfile - smoothedProfile; % We do this 
    
    % Step 3: Compute the Fourier Transform of the residual
    N               = length(residual);
    Y               = fft(residual);
    P2              = abs(Y/N);       % Two-sided spectrum
    P1              = P2(1:N/2+1);    % Single-sided spectrum
    P1(2:end-1)     = 2*P1(2:end-1);  % Correct for the energy in the negative frequencies
    
    % Step 4: Check for significant peaks in the Fourier spectrum
    % We check if the peak frequency is above a certain threshold
    threshold       = 1E-3;           % This can be adjusted based on the expected modulation strength
    peakValue       = max(P1);
    
    if peakValue > threshold
        ModulationFlag = true;   % Indicates sinusoidal modulation
    else
        ModulationFlag = false;  % Indicates otherwise
    end
end