Calculations/ODT-Calculator/+Potentials/generateAstigmaticCrossedBeamPotential.m

39 lines
2.0 KiB
Matlab

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