function U = generateModulatedSingleGaussianBeamPotential(positions, waists, P, options)
    
    alpha      = options.Polarizability;
    wavelength = options.Wavelength;
    mod_amp    = options.ModulationAmplitude;

    % Define constants
    eps0 = 8.854187817e-12;    % Permittivity of free space (F/m)
    c    = 299792458;          % Speed of light in vacuum (m/s)
    a0   = 5.29177210903e-11;  % Bohr radius (m)

    % Modulation amplitude
    mod_amp = mod_amp * waists(1);

    % Number of points in positions
    n_points = length(positions(1, :));

    % Calculate modulation function
    [dx, x_mod] = Helpers.calculateModulationFunction(mod_amp, n_points, 'arccos');

    % Calculate amplitude A
    A = 2 * P ./ (pi * Calculator.calculateWaist(positions(2, :), waists(1), wavelength) .* Calculator.calculateWaist(positions(2, :), waists(2), wavelength));

    % U_tilde constant term (same as before)
    U_tilde = (1 / (2 * eps0 * c)) * alpha * (4 * pi * eps0 * a0^3);

    % Initialize dU matrix
    dU = zeros(2 * n_points, n_points);

    % Compute dU values for each modulated x position
    for i = 1:length(x_mod)
        dU = [dU; exp(-2 * ((x_mod(i) - positions(1, :)) ./ Calculator.calculateWaist(positions(2, :), waists(1), wavelength)).^2)];
    end

    % Calculate the potential U using the trapezoidal rule for numerical integration
    U = - U_tilde * A .* (1 / (2 * mod_amp)) .* trapz(dx, dU, 1); % Integration along the first dimension (axis=0)
end