New script to propagate Gaussian beams using ABCD matrices and the q-parameter, adjusted parameters for variational solver, minor additions to script to make estimations for static lattice.
This commit is contained in:
parent
8352c13059
commit
897743531e
@ -23,7 +23,7 @@ OptionsStruct.ResidualTolerance = 1E-05;
|
||||
OptionsStruct.NoiseScaleFactor = 0.05;
|
||||
|
||||
OptionsStruct.MaxIterations = 10;
|
||||
OptionsStruct.VariationalWidth = 1.3;
|
||||
OptionsStruct.VariationalWidth = 1.5;
|
||||
OptionsStruct.WidthLowerBound = 0.01;
|
||||
OptionsStruct.WidthUpperBound = 12;
|
||||
OptionsStruct.WidthCutoff = 5e-3;
|
||||
@ -50,7 +50,7 @@ OptionsStruct = struct;
|
||||
OptionsStruct.NumberOfAtoms = 101250;
|
||||
OptionsStruct.DipolarPolarAngle = deg2rad(15);
|
||||
OptionsStruct.DipolarAzimuthAngle = 0;
|
||||
OptionsStruct.ScatteringLength = 65.00;
|
||||
OptionsStruct.ScatteringLength = 71.00;
|
||||
|
||||
OptionsStruct.TrapFrequencies = [50, 50, 500];
|
||||
OptionsStruct.TrapPotentialType = 'Harmonic';
|
||||
@ -65,7 +65,7 @@ OptionsStruct.ResidualTolerance = 1E-05;
|
||||
OptionsStruct.NoiseScaleFactor = 0.05;
|
||||
|
||||
OptionsStruct.MaxIterations = 10;
|
||||
OptionsStruct.VariationalWidth = 1.3;
|
||||
OptionsStruct.VariationalWidth = 1.5;
|
||||
OptionsStruct.WidthLowerBound = 0.01;
|
||||
OptionsStruct.WidthUpperBound = 12;
|
||||
OptionsStruct.WidthCutoff = 5e-3;
|
||||
|
||||
@ -33,7 +33,7 @@ a_s = 100 * BohrRadius; % Scattering length in meters
|
||||
N = 1e5; % Number of atoms
|
||||
v_x = 10;
|
||||
v_y = 10;
|
||||
v_z = 1E3;
|
||||
v_z = 1500;
|
||||
|
||||
[R_x, R_y, R_z] = calculateThomasFermiRadiusFromTrapParameters(v_x, v_y, v_z, N, a_s, Dy164Mass);
|
||||
fprintf('Thomas-Fermi radii: R_x = %.2f µm, R_y = %.2f µm, R_z = %.2f µm\n', R_x*1e6, R_y*1e6, R_z*1e6);
|
||||
@ -41,24 +41,24 @@ fprintf('Thomas-Fermi radii: R_x = %.2f µm, R_y = %.2f µm, R_z = %.2f µm\n',
|
||||
%%
|
||||
|
||||
% Define the parameters
|
||||
Power_values = linspace(1, 10, 50); % Laser powers (W), from 1 W to 10 W
|
||||
waist_values = linspace(10e-6, 100e-6, 50); % Beam waists (m), from 10 µm to 100 µm
|
||||
Power_values = linspace(0.1, 1.5, 50); % Laser powers (W), from 1 W to 10 W
|
||||
waist_values = linspace(50e-6, 90e-6, 50); % Beam waists (m), from 10 µm to 100 µm
|
||||
|
||||
% Constants
|
||||
a_s = 100 * BohrRadius; % Scattering length in meters
|
||||
N = 1e5; % Number of atoms
|
||||
AngleBetweenBeams = 9.44;
|
||||
Wavelength = 532e-9; % Laser wavelength in meters
|
||||
LatticeConstant = 13e-6;
|
||||
alpha = 180 * (AtomicUnitOfPolarizability / (2 * SpeedOfLight * VacuumPermittivity));
|
||||
R_z_matrix = zeros(length(waist_values), length(Power_values)); % Preallocate a matrix to store the R_z values
|
||||
|
||||
R_z_matrix = zeros(length(waist_values), length(Power_values)); % Preallocate a matrix to store the R_z values
|
||||
omega_eff_matrix = zeros(length(waist_values), length(Power_values)); % Preallocate a matrix to store the omega_z values
|
||||
% Loop over all combinations of P and w0
|
||||
for i = 1:length(waist_values)
|
||||
for j = 1:length(Power_values)
|
||||
% Calculate the Thomas-Fermi radii for each P and w0
|
||||
[R_z, LatticeConstant, AxialTrapFrequency] = calculateThomasFermiRadiusFromLatticeParameters(Power_values(j), waist_values(i), AngleBetweenBeams, Wavelength, N, a_s, Dy164Mass, alpha);
|
||||
[R_z, omega_eff] = calculateThomasFermiRadiusFromLatticeParameters(Power_values(j), waist_values(i), LatticeConstant, N, a_s, Dy164Mass, alpha);
|
||||
% Store the result in the matrix
|
||||
R_z_matrix(i, j) = R_z * 1e6; % Convert R_z to µm
|
||||
omega_eff_matrix(i, j) = (omega_eff/(2*pi)) * 1e-3; % Convert omega_z to kHz
|
||||
end
|
||||
end
|
||||
|
||||
@ -75,8 +75,22 @@ cbar = colorbar; % Add a colorbar
|
||||
xlabel('Laser Power (W)','FontSize',16);
|
||||
ylabel('Beam Waist (\mum)', 'Interpreter', 'tex','FontSize',16);
|
||||
ylabel(cbar,'R_z (\mum)', 'Interpreter', 'tex','FontSize',16,'Rotation',270)
|
||||
title(['\bf d = ' num2str(LatticeConstant * 1E6) ' \mum ; \bf \nu_z = ' num2str(round((AxialTrapFrequency/ (2*pi)) * 1E-3),3) ' kHz'],'fontsize',16, 'Interpreter', 'tex');
|
||||
title(['\bf Thomas-Fermi Radius @ Lattice Spacing = ' num2str(LatticeConstant * 1E6) ' \mum'],'fontsize',16, 'Interpreter', 'tex');
|
||||
|
||||
% Plot the resulting matrix
|
||||
figure(2)
|
||||
clf
|
||||
set(gcf,'Position',[50 50 950 750])
|
||||
imagesc(Power_values, waist_values * 1e6, omega_eff_matrix); % Specify x and y data for axes
|
||||
hold on
|
||||
[contour_matrix, contour_handle] = contour(Power_values, waist_values * 1e6, omega_eff_matrix, 'LineColor', 'white', 'LineWidth', 1.5);
|
||||
clabel(contour_matrix, contour_handle, 'FontSize', 12, 'Color', 'black', 'LabelSpacing', 600); % Adjust LabelSpacing for better placement
|
||||
set(gca, 'YDir', 'normal'); % Correct the y-axis direction
|
||||
cbar = colorbar; % Add a colorbar
|
||||
xlabel('Laser Power (W)','FontSize',16);
|
||||
ylabel('Beam Waist (\mum)', 'Interpreter', 'tex','FontSize',16);
|
||||
ylabel(cbar,'\nu_z (kHz)', 'Interpreter', 'tex','FontSize',16,'Rotation',270)
|
||||
title(['\bf Vertical Trap Frequency @ Lattice Spacing = ' num2str(LatticeConstant * 1E6) ' \mum'],'fontsize',16, 'Interpreter', 'tex');
|
||||
|
||||
%%
|
||||
function [R_x, R_y, R_z] = calculateThomasFermiRadiusFromTrapParameters(v_x, v_y, v_z, N, a, m)
|
||||
@ -113,7 +127,7 @@ function [R_x, R_y, R_z] = calculateThomasFermiRadiusFromTrapParameters(v_x, v_y
|
||||
R_z = sqrt(2 * mu / (m * omega_z^2));
|
||||
end
|
||||
|
||||
function [R_TF, d_lat, omega_eff] = calculateThomasFermiRadiusFromLatticeParameters(P, w_z, theta, lambda, N, a, m, alpha)
|
||||
function [R_TF, omega_eff] = calculateThomasFermiRadiusFromLatticeParameters(P, w_z, d_lat, N, a, m, alpha)
|
||||
% Calculate Thomas-Fermi radius of a BEC in a Gaussian potential
|
||||
%
|
||||
% Inputs:
|
||||
@ -136,7 +150,7 @@ function [R_TF, d_lat, omega_eff] = calculateThomasFermiRadiusFromLatticeParamet
|
||||
V_o = 4 * U_single;
|
||||
|
||||
% Calculate lattice constant
|
||||
d_lat = lambda / (2 * cos(deg2rad(theta)/2));
|
||||
% d_lat = lambda / (2 * cos(deg2rad(theta)/2));
|
||||
|
||||
% Calculate effective trapping frequency
|
||||
omega_eff = (pi / d_lat) * sqrt(2 * V_o / m);
|
||||
|
||||
300
Estimations/GaussianABCDCalculator.m
Normal file
300
Estimations/GaussianABCDCalculator.m
Normal file
@ -0,0 +1,300 @@
|
||||
%% Define a beam and visualize it propagating
|
||||
% Define "optical system" using matrices. For this first example, we will just use a free space propagation matrix, and let the beam propagate a distance
|
||||
|
||||
w0 = 1E-3; % 1mm beam waist
|
||||
lam = 421E-9; % wavelength
|
||||
zR = Zr(w0, lam); % Rayleigh range in m
|
||||
z0 = 0; % location of waist in m
|
||||
|
||||
syms d; % Define 'd' as a symbolic variable
|
||||
M = propagator(d);
|
||||
|
||||
[R, w] = q1_inv_func(0, w0, lam, M); % Do all the ABCD and q-parameter math, and return the waist and radius of curvature functions
|
||||
|
||||
plotResult(w, d, linspace(0,10))
|
||||
|
||||
%% Add a lens to the system
|
||||
|
||||
w0 = 1E-3; % 1mm beam waist
|
||||
lam = 421E-9; % wavelength
|
||||
zR = Zr(w0, lam); % Rayleigh range in m
|
||||
z0 = 0; % location of waist in m
|
||||
|
||||
syms d; % Define 'd' as a symbolic variable
|
||||
M = matmultiplier(propagator(d), lens(.5), propagator(1));
|
||||
|
||||
[R, w] = q1_inv_func(0, w0, lam, M);
|
||||
|
||||
plotResult(w, d, linspace(0,1,50))
|
||||
|
||||
%% Add another to make a beam expander, that enlarges the beam waist by a factor 5 and then collimates it
|
||||
|
||||
w0 = 1E-3; % 1mm beam waist
|
||||
lam = 421E-9; % wavelength
|
||||
zR = Zr(w0, lam); % Rayleigh range in m
|
||||
z0 = 0; % location of waist in m
|
||||
expansion_factor = 5;
|
||||
|
||||
syms d1 d2 d3 f1 f2
|
||||
|
||||
M = matmultiplier(propagator(d3), lens(f2), propagator(d2), lens(f1), propagator(d1));
|
||||
|
||||
[R, w] = q1_inv_func(0, w0, lam, M);
|
||||
|
||||
% Substitute values into w
|
||||
w = subs(w, {d1, d3}, {1, 0});
|
||||
|
||||
% Solve for f1
|
||||
f1_eq = solve(w - expansion_factor*w0 == 0, f1);
|
||||
f1_eq_subs = subs(f1_eq, d2, 1);
|
||||
% Convert solutions to numeric values
|
||||
f1_val_numeric = double(f1_eq_subs);
|
||||
fprintf('f1 = %.2f m or %.2f m, for a lens separation of 1 meter\n', f1_val_numeric(1), f1_val_numeric(2))
|
||||
|
||||
for i = 1:length(f1_val_numeric)
|
||||
|
||||
f1_val = f1_val_numeric(i);
|
||||
R = subs(R, {d1, d2, d3, f1}, {1, 1, 0, f1_val});
|
||||
f2_val = solve(1/R == 0, f2);
|
||||
fprintf('f2 = %.2f, for a collimated beam, 5x the original waist, after propagating 1m to the first lens of f1 = %.2f, and propagating another 1m to the second lens\n', double(f2_val), double(f1_val))
|
||||
end
|
||||
|
||||
chosen_f1_val = f1_val_numeric(1);
|
||||
|
||||
fprintf('Chosen f1 = %.2f m \n', chosen_f1_val)
|
||||
|
||||
M = matmultiplier(propagator(d3), lens(f2_val), propagator(1), lens(chosen_f1_val), propagator(1));
|
||||
|
||||
[R, w] = q1_inv_func(0, w0, lam, M);
|
||||
|
||||
plotResult(w,d3)
|
||||
|
||||
estimated_expansion_factor = subs(w, d3, 0)/ w0;
|
||||
fprintf('beam is w = %.2f x w0\n', double(estimated_expansion_factor))
|
||||
|
||||
beam_size_change = ((subs(w,d3,10) - subs(w,d3,0)) / subs(w,d3,0)) * 100;
|
||||
fprintf('Over 10 m after second lens, beam changes by %.5f percent\n', double(beam_size_change))
|
||||
|
||||
%% Function definitions
|
||||
|
||||
% Input Ray parameter, i.e. height and angle
|
||||
function ret = ray(y, theta)
|
||||
% Parameters
|
||||
% ----------
|
||||
% y : float or integer in meters
|
||||
% The vertical height of a ray.
|
||||
% theta : float or integer in radians
|
||||
% The angle of divergence of the ray.
|
||||
%
|
||||
% Returns
|
||||
% -------
|
||||
% mat : 2x1 matrix
|
||||
% [
|
||||
% [y],
|
||||
% [theta]
|
||||
% ]
|
||||
|
||||
ret = [y; theta];
|
||||
end
|
||||
|
||||
% Ray Transfer Matrix for ideal lens with focal length f
|
||||
function ret = lens(f)
|
||||
% Parameters
|
||||
% ----------
|
||||
% f : float or integer in meters
|
||||
% Thin lens focal length in meters
|
||||
%
|
||||
% Returns
|
||||
% -------
|
||||
% mat : 2x2 matrix
|
||||
% [
|
||||
% [ 1, 0],
|
||||
% [-1/f, 1]
|
||||
% ]
|
||||
|
||||
ret = [1, 0; -1/f, 1];
|
||||
end
|
||||
|
||||
% Ray Transfer Matrix for propagation over a distance d
|
||||
function ret = propagator(d)
|
||||
% Parameters
|
||||
% ----------
|
||||
% d : float or integer
|
||||
% Distance light is propagating along the z-axis.
|
||||
%
|
||||
% Returns
|
||||
% -------
|
||||
% mat: 2x2 matrix
|
||||
% [
|
||||
% [1, d],
|
||||
% [0, 1]
|
||||
% ]
|
||||
|
||||
ret = [1, d; 0, 1];
|
||||
end
|
||||
|
||||
% Multiplying the matrices together. mat1 is the last matrix the light interacts with.
|
||||
function ret = matmultiplier(mat, varargin)
|
||||
% Parameters
|
||||
% ----------
|
||||
% mat1 : 2x2 ABCD matrix
|
||||
% Last matrix light interacts with.
|
||||
% varargin : 2x2 ABCD matrices
|
||||
% From left to right, the matrices should be entered such that the leftmost matrix interacts
|
||||
% with light temporally after the rightmost matrix.
|
||||
%
|
||||
% Returns
|
||||
% -------
|
||||
% Mat : 2x2 matrix
|
||||
% The ABCD matrix describing the whole optical system.
|
||||
|
||||
ret = mat;
|
||||
for i = 1:length(varargin)
|
||||
ret = ret * varargin{i};
|
||||
end
|
||||
end
|
||||
|
||||
% Adding Gaussian beam parameters
|
||||
function zr = Zr(wo, lam)
|
||||
% Parameters
|
||||
% ----------
|
||||
% wo : float or integer
|
||||
% Beam waist radius in meters.
|
||||
% lam : float or integer
|
||||
% Wavelength of light in meters.
|
||||
%
|
||||
% Returns
|
||||
% -------
|
||||
% zr : float or integer
|
||||
% Rayleigh range for given beam waist and wavelength.
|
||||
|
||||
zr = pi * wo^2 / lam;
|
||||
end
|
||||
|
||||
% Calculate beam waist radius
|
||||
function w0 = W0(zr, lam)
|
||||
% Parameters
|
||||
% ----------
|
||||
% zr : float or integer
|
||||
% Rayleigh range in meters.
|
||||
% lam : float or integer
|
||||
% Wavelength of light in meters.
|
||||
%
|
||||
% Returns
|
||||
% -------
|
||||
% w0 : float or integer
|
||||
% Beam waist radius in meters.
|
||||
|
||||
w0 = sqrt(lam * zr / pi);
|
||||
end
|
||||
|
||||
function [z_out, zr_out] = q1_func(z, w0, lam, mat)
|
||||
% Parameters
|
||||
% ----------
|
||||
% z : float or integer
|
||||
% Position of the beam waist in meters.
|
||||
% w0 : float or integer
|
||||
% Radial waist size in meters (of the embedded Gaussian, i.e. W0/M).
|
||||
% lam : float or integer
|
||||
% Wavelength of light in meters.
|
||||
% mat : 2x2 matrix
|
||||
% The ABCD matrix describing the optical system.
|
||||
%
|
||||
% Returns
|
||||
% -------
|
||||
% z_out : float or integer
|
||||
% Position of the beam waist after the optical system.
|
||||
% zr_out : float or integer
|
||||
% Rayleigh range of the beam after the optical system.
|
||||
|
||||
A = mat(1, 1);
|
||||
B = mat(1, 2);
|
||||
C = mat(2, 1);
|
||||
D = mat(2, 2);
|
||||
|
||||
% Calculate Rayleigh range for the given beam waist and wavelength
|
||||
zr = Zr(w0, lam);
|
||||
|
||||
% Calculate real and imaginary parts
|
||||
real_part = (A * C * (z^2 + zr^2) + z * (A * D + B * C) + B * D) / (C^2 * (z^2 + zr^2) + 2 * C * D * z + D^2);
|
||||
imag_part = (zr * (A * D - B * C)) / (C^2 * (z^2 + zr^2) + 2 * C * D * z + D^2);
|
||||
|
||||
% Output the new position and Rayleigh range
|
||||
z_out = real_part;
|
||||
zr_out = imag_part;
|
||||
end
|
||||
|
||||
function [R, w] = q1_inv_func(z, w0, lam, mat)
|
||||
% Parameters
|
||||
% ----------
|
||||
% z : float or integer
|
||||
% Position of the beam waist in meters.
|
||||
% w0 : float or integer
|
||||
% Radial waist size in meters (of the embedded Gaussian, i.e. W0/M).
|
||||
% lam : float or integer
|
||||
% Wavelength of light in meters.
|
||||
% mat : 2x2 matrix
|
||||
% The ABCD matrix describing the optical system.
|
||||
%
|
||||
% Returns
|
||||
% -------
|
||||
% R : float or integer
|
||||
% Radius of curvature of the wavefront in meters.
|
||||
% w : float or integer
|
||||
% Radius of the beam in meters.
|
||||
|
||||
A = mat(1, 1);
|
||||
B = mat(1, 2);
|
||||
C = mat(2, 1);
|
||||
D = mat(2, 2);
|
||||
|
||||
% Calculate Rayleigh range for the given beam waist and wavelength
|
||||
zr = Zr(w0, lam);
|
||||
|
||||
% Calculate real and imaginary parts
|
||||
real_part = (A * C * (z^2 + zr^2) + z * (A * D + B * C) + B * D) / (A^2 * (z^2 + zr^2) + 2 * A * B * z + B^2);
|
||||
imag_part = -zr * ((A * D) - (B * C)) / (A^2 * (z^2 + zr^2) + (2 * A * B * z) + B^2);
|
||||
|
||||
% Calculate radius of curvature and beam radius
|
||||
R = 1/real_part;
|
||||
w = sqrt(-lam / imag_part / pi);
|
||||
end
|
||||
|
||||
function plotResult(func, var, rang)
|
||||
% Parameters
|
||||
% ----------
|
||||
% func : symfun (symbolic function of one variable)
|
||||
% Symbolic function defining the beam width after the last optical element.
|
||||
% var : symbolic variable
|
||||
% Variable in func that will be plotted.
|
||||
% rang : array (optional)
|
||||
% Array of values along the optical axis to be plotted. Default range is 0 to 3 with step size 0.01.
|
||||
|
||||
% Set default range if not provided
|
||||
if nargin < 3
|
||||
rang = 0:0.01:3;
|
||||
end
|
||||
|
||||
% Create a numeric function from the symbolic function
|
||||
func_handle = matlabFunction(func, 'Vars', var);
|
||||
|
||||
% Compute the values of the function and its negative
|
||||
y_vals = func_handle(rang) * 1E3; % in mm
|
||||
y_neg_vals = -y_vals;
|
||||
|
||||
% Create a figure
|
||||
figure;
|
||||
set(gcf, 'Position', [100 100 950 750]);
|
||||
|
||||
% Plot the function and its negative
|
||||
plot(rang, y_vals, 'b'); hold on;
|
||||
plot(rang, y_neg_vals, 'b');
|
||||
|
||||
% Fill the area between the curves with translucent blue color
|
||||
fill([rang, fliplr(rang)], [y_vals, fliplr(y_neg_vals)], 'b', 'FaceAlpha', 0.2, 'EdgeColor', 'none');
|
||||
|
||||
% Add grid and labels
|
||||
grid on;
|
||||
xlabel('Optical Axis (m)', 'FontSize', 16);
|
||||
ylabel('Beam size (mm)', 'FontSize', 16);
|
||||
end
|
||||
Loading…
Reference in New Issue
Block a user