function U = generateAstigmaticCrossedBeamPotential(positions, waists, P, options)
    
    alpha      = options.Polarizability;
    wavelength = options.Wavelength;
    delta      = options.delta;
    del_y      = options.foci_disp_crossed;
    del_y_1 = del_y(1);
    del_y_2 = del_y(2);

    foci_shift = options.foci_shift;
    focus_shift_beam_1 = foci_shift(1);
    focus_shift_beam_2 = foci_shift(2);

    beam_disp = options.beam_disp;
    beam_1_disp = repmat(beam_disp{1}, 1, size(positions, 2));
    beam_2_disp = repmat(beam_disp{2}, 1, size(positions, 2));

    % Calculate beam 1 potential
    beam_1_positions = positions + beam_1_disp;
    A_1 = 2*P(1) / (pi * w(beam_1_positions(2,:) - (del_y_1/2) + focus_shift_beam_1, waists{1}(1), wavelength) ...
                     * w(beam_1_positions(2,:) + (del_y_1/2) + focus_shift_beam_1, waists{1}(2), wavelength));
    U_1_tilde = (1 / (2 * eps0 * c)) * alpha * (4 * pi * eps0 * a0^3);
    U_1 = - U_1_tilde * A_1 .* exp(-2 * ((beam_1_positions(1,:) ./ w(beam_1_positions(2,:) - (del_y_1/2) + focus_shift_beam_1, waists{1}(1), wavelength)).^2 ...
                    + (beam_1_positions(3,:) ./ w(beam_1_positions(2,:) + (del_y_1/2) + focus_shift_beam_1, waists{1}(2), wavelength)).^2));

    % Rotation matrix for beam 2
    R = rotation_matrix([0, 0, 1], deg2rad(delta));
    beam_2_positions = R * (positions + beam_2_disp);

    % Calculate beam 2 potential
    A_2 = 2*P(2) / (pi * w(beam_2_positions(2,:) - (del_y_2/2) + focus_shift_beam_2, waists{2}(1), wavelength) ...
                     * w(beam_2_positions(2,:) + (del_y_2/2) + focus_shift_beam_2, waists{2}(2), wavelength));
    U_2_tilde = (1 / (2 * eps0 * c)) * alpha * (4 * pi * eps0 * a0^3);
    U_2 = - U_2_tilde * A_2 .* exp(-2 * ((beam_2_positions(1,:) ./ w(beam_2_positions(2,:) - (del_y_2/2) + focus_shift_beam_2, waists{2}(1), wavelength)).^2 ...
                    + (beam_2_positions(3,:) ./ w(beam_2_positions(2,:) + (del_y_2/2) + focus_shift_beam_2, waists{2}(2), wavelength)).^2));

    % Total potential
    U = U_1 + U_2;
end