function U = generateAstigmaticSingleGaussianBeamPotential(positions, waists, P, options)
    
    alpha      = options.Polarizability;
    wavelength = options.Wavelength;
    del_y      = options.foci_disp_crossed;

    % Calculate the beam waists at adjusted positions (due to astigmatism)
    w_x = w(positions(2, :) - (del_y / 2), waists(1), wavelength); % Waist in x direction
    w_z = w(positions(2, :) + (del_y / 2), waists(2), wavelength); % Waist in z direction

    % Calculate amplitude A
    A = 2 * P ./ (pi * w_x .* w_z);

    % U_tilde is a constant term (relating to polarizability, permittivity, and speed of light)
    U_tilde = (1 / (2 * eps0 * c)) * alpha * (4 * pi * eps0 * a0^3);

    % Gaussian potential calculation
    U = - U_tilde * A .* exp(-2 * ((positions(1, :) ./ w_x).^2 + (positions(3, :) ./ w_z).^2));
end