From 897743531e3312f917bd8e8b40d70109b6b1cb65 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Thu, 6 Feb 2025 20:57:53 +0100 Subject: [PATCH] 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. --- .../+Scripts/run_on_cluster_in_plane_trap.m | 6 +- Estimations/EstimatesForStaticLattice.m | 36 ++- Estimations/GaussianABCDCalculator.m | 300 ++++++++++++++++++ 3 files changed, 328 insertions(+), 14 deletions(-) create mode 100644 Estimations/GaussianABCDCalculator.m diff --git a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_in_plane_trap.m b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_in_plane_trap.m index 9928276..959c0ec 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_in_plane_trap.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_in_plane_trap.m @@ -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; diff --git a/Estimations/EstimatesForStaticLattice.m b/Estimations/EstimatesForStaticLattice.m index 3b13a8b..da9c868 100644 --- a/Estimations/EstimatesForStaticLattice.m +++ b/Estimations/EstimatesForStaticLattice.m @@ -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); diff --git a/Estimations/GaussianABCDCalculator.m b/Estimations/GaussianABCDCalculator.m new file mode 100644 index 0000000..084104d --- /dev/null +++ b/Estimations/GaussianABCDCalculator.m @@ -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