From 0647b032d09c60b5a2b847fb79c9e8fd43cee6d7 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Fri, 7 Feb 2025 17:58:19 +0100 Subject: [PATCH] MAJOR UPDATE: Corrected script to extract lattice properties and improved plotting routines. --- .../+Scripts/analyzeGSWavefunction.m | 27 ++- .../+Scripts/extractLatticeProperties.m | 194 +++++++++++++----- Dipolar-Gas-Simulator/+Scripts/run_locally.m | 6 +- .../+Scripts/run_on_cluster.m | 12 +- .../+Scripts/run_on_cluster_in_plane_trap.m | 6 +- 5 files changed, 176 insertions(+), 69 deletions(-) diff --git a/Dipolar-Gas-Simulator/+Scripts/analyzeGSWavefunction.m b/Dipolar-Gas-Simulator/+Scripts/analyzeGSWavefunction.m index b3841e2..b23c0f3 100644 --- a/Dipolar-Gas-Simulator/+Scripts/analyzeGSWavefunction.m +++ b/Dipolar-Gas-Simulator/+Scripts/analyzeGSWavefunction.m @@ -56,7 +56,7 @@ function [contrast, periodX, periodY] = analyzeGSWavefunction(folder_path, run_i %------------------ Lattice Properties ------------------ % - [kx, ky, fftMagnitude, lattice_type, periodX, periodY, freq_x, freq_y, rotation_angle] = Scripts.extractLatticeProperties(nxyScaled, x, y); + [kx, ky, fftMagnitude, lattice_type, periodX, periodY, freq_x, freq_y] = Scripts.extractLatticeProperties(nxyScaled, x, y); %------------------ Plotting ------------------ % @@ -72,10 +72,31 @@ function [contrast, periodX, periodY] = analyzeGSWavefunction(folder_path, run_i cbar1.Label.Interpreter = 'latex'; colormap(gca, Helper.Colormaps.plasma()) % clim(ax1,[0.00,0.3]); + + % Define normalized positions (relative to axis limits) + x_offset = 0.025; % 5% offset from the edges + y_offset = 0.025; % 5% offset from the edges + + % Top-left corner (normalized axis coordinates) + text(0 + x_offset, 1 - y_offset, lattice_type, ... + 'Color', 'white', 'FontWeight', 'bold', 'Interpreter', 'tex', 'FontSize', 30, 'Units', 'normalized', 'HorizontalAlignment', 'left', 'VerticalAlignment', 'top'); + + % Top-right corner (normalized axis coordinates) + text(1 - x_offset, 1 - y_offset, ['C: ', num2str(contrast, '%.3f')], ... + 'Color', 'white', 'FontWeight', 'bold', 'Interpreter', 'tex', 'FontSize', 30, 'Units', 'normalized', 'HorizontalAlignment', 'right', 'VerticalAlignment', 'top'); + + % Bottom-left corner (normalized axis coordinates) + text(0 + x_offset, 0 + y_offset, ['dx: ', num2str(periodX, '%.2f'), ' \mum'], ... + 'Color', 'white', 'FontWeight', 'bold', 'Interpreter', 'tex', 'FontSize', 30, 'Units', 'normalized', 'HorizontalAlignment', 'left', 'VerticalAlignment', 'bottom'); + + % Bottom-right corner (normalized axis coordinates) + text(1 - x_offset, 0 + y_offset, ['dy: ', num2str(periodY, '%.2f'), ' \mum'], ... + 'Color', 'white', 'FontWeight', 'bold', 'Interpreter', 'tex', 'FontSize', 30, 'Units', 'normalized', 'HorizontalAlignment', 'right', 'VerticalAlignment', 'bottom'); + ylabel(cbar1,'$na_{dd}^2$','FontSize',16,'Rotation',270) xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14) ylabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14) - title(['$|\Psi(x,y)|^2$ - Contrast: ', num2str(contrast, '%.3f'), '; Period X = ', num2str(periodX, '%.2f'), '$ \mu$m', '; Period Y = ', num2str(periodY, '%.2f'), '$ \mu$m'], 'Interpreter', 'latex', 'FontSize', 14) + title('$|\Psi(x,y)|^2$ ', 'Interpreter', 'latex', 'FontSize', 14) % Plot 2-D FFT with detected peaks nexttile; % Equivalent to subplot('Position', [0.05, 0.55, 0.28, 0.4]); @@ -90,6 +111,6 @@ function [contrast, periodX, periodY] = analyzeGSWavefunction(folder_path, run_i ylabel('$k_y l_o$', 'Interpreter', 'latex', 'FontSize', 14) title('$\mathcal{F}\{|\Psi(x,y)|^2\}$', 'Interpreter', 'latex', 'FontSize', 16); - sgtitle(['$\omega_z = 2 \pi \times$', num2str(round(Params.wz/(2*pi)), '%.2f'), ' Hz; $\theta = $', num2str(rad2deg(Params.theta), '%.2f'), '$^\circ$; ', sprintf('Detected Lattice Type: %s\n', lattice_type)], 'Interpreter', 'latex', 'FontSize', 18, 'FontWeight', 'bold') + sgtitle(['$\omega_z = 2 \pi \times$', num2str(round(Params.wz/(2*pi)), '%.2f'), ' Hz; $\theta = $', num2str(rad2deg(Params.theta), '%.2f'), '$^\circ$'], 'Interpreter', 'latex', 'FontSize', 18, 'FontWeight', 'bold') end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Scripts/extractLatticeProperties.m b/Dipolar-Gas-Simulator/+Scripts/extractLatticeProperties.m index 69a85e7..a007a6f 100644 --- a/Dipolar-Gas-Simulator/+Scripts/extractLatticeProperties.m +++ b/Dipolar-Gas-Simulator/+Scripts/extractLatticeProperties.m @@ -1,6 +1,6 @@ -function [kx, ky, fftMagnitude, lattice_type, dx_um, dy_um, freq_x, freq_y, rotation_angle] = extractLatticeProperties(I, x, y) +function [kx, ky, fftMagnitude, lattice_type, dx_um, dy_um, freq_x, freq_y] = extractLatticeProperties(I, x, y) % Detects lattice geometry, extracts periodic spacings, and reconstructs real-space lattice vectors. - % Handles arbitrary lattice geometries with rotation correction. + % Handles arbitrary lattice geometries % % Inputs: % I - Grayscale image with periodic structures. @@ -12,7 +12,7 @@ function [kx, ky, fftMagnitude, lattice_type, dx_um, dy_um, freq_x, freq_y, rota % dx_um - Spacing along x-axis in micrometers. % dy_um - Spacing along y-axis in micrometers. % real_lattice_vectors - [2x2] matrix of lattice vectors in micrometers. - % rotation_angle - Angle of rotation of the lattice in degrees. + % angle_in_reciprocal_space - Angle of the reciprocal lattice primitive vectors in degrees. % Compute 2D Fourier Transform F = fft2(double(I)); @@ -64,7 +64,7 @@ function [kx, ky, fftMagnitude, lattice_type, dx_um, dy_um, freq_x, freq_y, rota dy_um = NaN; freq_x = NaN; freq_y = NaN; - rotation_angle = NaN; + angle_in_reciprocal_space = NaN; else % Select two shortest independent lattice vectors - Reciprocal lattice vectors @@ -72,7 +72,7 @@ function [kx, ky, fftMagnitude, lattice_type, dx_um, dy_um, freq_x, freq_y, rota G1 = unique_vectors(1, :); % Reciprocal lattice vector 1 G2 = unique_vectors(3, :); % Reciprocal lattice vector 2 - reciprocal_lattice_vectors = vertcat(G1, G2); + reciprocal_lattice_vectors = abs(vertcat(G1, G2)); % Calculate the angle between the reciprocal lattice vectors using dot product dotProduct = dot(G1, G2); % Dot product of G1 and G2 @@ -80,7 +80,7 @@ function [kx, ky, fftMagnitude, lattice_type, dx_um, dy_um, freq_x, freq_y, rota magnitudeG2 = norm(G2); % Magnitude of G2 % Angle between the reciprocal lattice vectors (in degrees) - rotation_angle = rad2deg(acos(dotProduct / (magnitudeG1 * magnitudeG2))); + angle_in_reciprocal_space = rad2deg(acos(dotProduct / (magnitudeG1 * magnitudeG2))); % Convert to real-space lattice vectors @@ -88,67 +88,153 @@ function [kx, ky, fftMagnitude, lattice_type, dx_um, dy_um, freq_x, freq_y, rota detG = G1(1) * G2(2) - G1(2) * G2(1); % Determinant of the matrix formed by G1 and G2 if detG == 0 - % Handle the case of collinear reciprocal lattice vectors - disp('Reciprocal lattice vectors are collinear.'); + % Handle the case of reciprocal lattice vector matrix being singular - % If both G1, G2 or just G1 are along x-direction - if G1(1) ~= 0 && G1(2) == 0 - a1 = [1 / norm(G1), 0]; % Real-space lattice vector along x - a2 = [0, 0]; % No periodicity in the y-direction (1D lattice) - % If G2 is along x-direction or if G1 is zero. - elseif G2(1) ~= 0 && G2(2) == 0 - a1 = [1 / norm(G2), 0]; % Real-space lattice vector along x - a2 = [0, 0]; % No periodicity in the y-direction - % If G1 and G2 are both in the y-direction (i.e., [0, non-zero]) - elseif G1(1) == 0 && G2(1) == 0 - a1 = [0, 1 / norm(G1)]; % Real-space lattice vector along y - a2 = [0, 0]; % No periodicity in the x-direction - % If both vectors are zero (no periodicity) - else - a1 = [0, 0]; - a2 = [0, 0]; + v1 = reciprocal_lattice_vectors(1, :); + v2 = reciprocal_lattice_vectors(2, :); + + % Check if either vector is the zero vector + if all(v1 == 0) && all(v2 == 0) + % If both vectors are zero, return the same zero matrix + real_lattice_vectors = [0, 0; 0, 0]; end + + % Compute the norms of the vectors + norm_v1 = norm(v1); + norm_v2 = norm(v2); + + % Find the vector with the smallest norm + if norm_v2 < norm_v1 + % If v2 has the smaller norm, set v1 to zero and keep v2 + reciprocal_lattice_vectors = [0, 0; v2]; + if v2(1) == 0 + real_lattice_vectors = [0, 0; 0, 1/v2(2)]; + elseif v2(2) == 0 + real_lattice_vectors = [0, 0; 1/v2(1), 0]; + else + real_lattice_vectors = [1/v2(1), 1/v2(2); 0, 0]; + end + else + % If v1 has the smaller norm, set v2 to zero and keep v1 + reciprocal_lattice_vectors = [v1; 0, 0]; + if v1(1) == 0 + real_lattice_vectors = [0, 1/v1(2); 0, 0]; + elseif v1(2) == 0 + real_lattice_vectors = [1/v1(1), 0; 0, 0]; + else + real_lattice_vectors = [1/v1(1), 1/v1(2); 0, 0]; + end + end + else % If reciprocal vectors are not collinear, compute the 2D real-space lattice - a1 = [G2(2), -G1(2)] / detG; - a2 = [-G2(1), G1(1)] / detG; + real_lattice_vectors = inv(reciprocal_lattice_vectors'); end - real_lattice_vectors = vertcat(a1, a2); - % Ensure correct orientation real_lattice_vectors(isinf(real_lattice_vectors) | isnan(real_lattice_vectors)) = 0; + % Separate the components of real space vectors + real_x = real_lattice_vectors(:, 1); + real_y = real_lattice_vectors(:, 2); + + % Compute projections + proj_real_x = real_lattice_vectors * [1; 0]; % Projections onto x-axis + [smaller_value, idx] = min(proj_real_x); + proj_real_x(idx) = []; % Remove the element + proj_real_y = real_lattice_vectors * [0; 1]; % Projections onto y-axis + [smaller_value, idx] = min(proj_real_y); + proj_real_y(idx) = []; % Remove the element + + % Separate the components of reciprocal space vectors + reciprocal_x = reciprocal_lattice_vectors(:, 1); + reciprocal_y = reciprocal_lattice_vectors(:, 2); + + % Compute projections + proj_reciprocal_x = reciprocal_lattice_vectors * [1; 0]; % Projections onto x-axis + [smaller_value, idx] = min(proj_reciprocal_x); + proj_reciprocal_x(idx) = []; % Remove the element + proj_reciprocal_y = reciprocal_lattice_vectors * [0; 1]; % Projections onto y-axis + [smaller_value, idx] = min(proj_reciprocal_y); + proj_reciprocal_y(idx) = []; % Remove the element + + %{ + % Create a figure for side-by-side subplots + figure('Position', [100, 100, 1600, 800]); + clf + t = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x3 grid + + % Plot real space vectors in the first subplot (2D case) + nexttile; % 1 row, 2 columns, first plot + quiver(0, 0, real_x(1) * 1E6, real_y(1) * 1E6, 'r', 'LineWidth', 2, 'MaxHeadSize', 0.5); + hold on; + quiver(0, 0, real_x(2) * 1E6, real_y(2) * 1E6, 'b', 'LineWidth', 2, 'MaxHeadSize', 0.5); + + % Plot projections of the real space vectors + quiver(0, 0, proj_real_x * 1E6, 0, 'k--', 'LineWidth', 2, 'MaxHeadSize', 0.5); % Projection on X-axis (horizontal dotted line) + quiver(0, 0, 0, proj_real_y * 1E6, 'k--', 'LineWidth', 2, 'MaxHeadSize', 0.5); % Projection on Y-axis (vertical dotted line) + + hold off; + xlabel('X'); + ylabel('Y'); + title('Real Space Primitive Vectors', 'FontSize', 14); + grid on; + axis equal; + + % Plot reciprocal space vectors in the second subplot (2D case) + nexttile; % 1 row, 2 columns, second plot + quiver(0, 0, reciprocal_x(1) * 1E-6, reciprocal_y(1) * 1E-6, 'r', 'LineWidth', 2, 'MaxHeadSize', 0.5); + hold on; + quiver(0, 0, reciprocal_x(2) * 1E-6, reciprocal_y(2) * 1E-6, 'b', 'LineWidth', 2, 'MaxHeadSize', 0.5); + + % Plot projections of the real space vectors + quiver(0, 0, proj_reciprocal_x * 1E-6, 0, 'k--', 'LineWidth', 2, 'MaxHeadSize', 0.5); % Projection on X-axis (horizontal dotted line) + quiver(0, 0, 0, proj_reciprocal_y * 1E-6, 'k--', 'LineWidth', 2, 'MaxHeadSize', 0.5); % Projection on Y-axis (vertical dotted line) + + hold off; + xlabel('X'); + ylabel('Y'); + title('Reciprocal Space Primitive Vectors', 'FontSize', 14); + grid on; + axis equal; + %} + % Compute lattice spacings - dx_um = (1/norm(G2)) * 1E6; - dy_um = (1/norm(G1)) * 1E6; - - % Correct for rotation by applying inverse rotation matrix - theta_rad = deg2rad(-rotation_angle); % Convert to radians - R = [cos(theta_rad), -sin(theta_rad); sin(theta_rad), cos(theta_rad)]; % Rotation matrix - real_lattice_vectors = (R * real_lattice_vectors')'; % Apply rotation correction - + dx_um = proj_real_x * 1E6; + dy_um = proj_real_y * 1E6; + % Classify lattice type based on angle symmetry - angle_between_vectors = atan2d(real_lattice_vectors(2,2), real_lattice_vectors(2,1)) - ... - atan2d(real_lattice_vectors(1,2), real_lattice_vectors(1,1)); - angle_between_vectors = mod(angle_between_vectors, 180); - - if abs(angle_between_vectors - 60) < 5 || abs(angle_between_vectors - 120) < 5 - lattice_type = 'Hexagonal'; - elseif abs(angle_between_vectors - 90) < 5 - lattice_type = 'Square'; - elseif abs(angle_between_vectors - 90) > 5 && abs(angle_between_vectors - 120) > 5 - lattice_type = 'Oblique'; + angle_in_real_space = abs(atan2d(real_lattice_vectors(2,2), real_lattice_vectors(2,1)) - ... + atan2d(real_lattice_vectors(1,2), real_lattice_vectors(1,1))); + angle_in_real_space_mod180 = mod(angle_in_real_space, 180); + + if all(real_lattice_vectors(1, :) == 0) || all(real_lattice_vectors(2, :) == 0) + lattice_type = 'Stripe'; + angle_in_real_space = NaN; + angle_in_reciprocal_space = NaN; else - lattice_type = 'Undetermined'; + if abs(angle_in_real_space_mod180 - 60) < 5 || abs(angle_in_real_space_mod180 - 120) < 5 + if abs(dx_um - dy_um) < 1E-6 + lattice_type = 'Triangular'; + else + lattice_type = 'Sheared Triangular'; + end + elseif abs(angle_in_real_space_mod180 - 90) < 5 + lattice_type = 'Square'; + elseif abs(angle_in_real_space_mod180 - 90) > 5 && abs(angle_in_real_space_mod180 - 120) > 5 + lattice_type = 'Sheared'; + else + lattice_type = 'Undetermined'; + end end - + % Display results - % fprintf('Detected Lattice Type: %s\n', lattice_type); - % fprintf('Estimated Spacing (dx): %.2f µm\n', dx_um); - % fprintf('Estimated Spacing (dy): %.2f µm\n', dy_um); - % fprintf('Rotation Angle: %.2f°\n', rotation_angle); - % fprintf('Lattice Vectors:\n'); - % disp(real_lattice_vectors); + fprintf('Detected Lattice Type: %s\n', lattice_type); + fprintf('Estimated Spacing (dx): %.2f µm\n', dx_um); + fprintf('Estimated Spacing (dy): %.2f µm\n', dy_um); + fprintf('Angle between real space lattice primitve vectors: %.2f°\n', angle_in_real_space); + fprintf('Angle between reciprocal lattice primitve vectors: %.2f°\n', angle_in_reciprocal_space); + fprintf('Real Space Primitive Vectors:\n'); + disp(real_lattice_vectors); end end diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m index 4f362e6..6263b9a 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -251,7 +251,7 @@ Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber) %% - Analysis SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz500'; -JobNumber = 1; +JobNumber = 2; % Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber) [contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction(SaveDirectory, JobNumber); @@ -270,7 +270,7 @@ JobNumber = 2; %% - Analysis SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz2000'; JobNumber = 2; -Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber) +% Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber) [contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction(SaveDirectory, JobNumber); %% - Analysis @@ -298,6 +298,6 @@ SaveOption = 'images'; [contrast_array, periodX_array, periodY_array] = Scripts.analyzeGSWavefunction_multipleruns(SaveDirectory, NumberOfJobs, SaveOption); %% - Analysis SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/Hz500'; -JobNumber = 0; +JobNumber = 1; % Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber) [contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction(SaveDirectory, JobNumber); diff --git a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m index 01a456d..fe38444 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m @@ -1,5 +1,5 @@ %% Tilting of the dipoles -% Atom Number = 1250 ppum +% With an in-plane harmonic trap %% v_z = 500, theta = 0 @@ -17,13 +17,13 @@ OptionsStruct.NumberOfGridPoints = [256, 256]; OptionsStruct.Dimensions = [35, 35]; OptionsStruct.TimeStepSize = 0.005; % in s OptionsStruct.MinimumTimeStepSize = 1E-5; % in s -OptionsStruct.TimeCutOff = 5E5; % in s +OptionsStruct.TimeCutOff = 1E6; % in s OptionsStruct.EnergyTolerance = 5E-10; 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'; @@ -59,13 +59,13 @@ OptionsStruct.NumberOfGridPoints = [256, 256]; OptionsStruct.Dimensions = [35, 35]; OptionsStruct.TimeStepSize = 0.005; % in s OptionsStruct.MinimumTimeStepSize = 1E-5; % in s -OptionsStruct.TimeCutOff = 5E5; % in s +OptionsStruct.TimeCutOff = 1E6; % in s OptionsStruct.EnergyTolerance = 5E-10; 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/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_in_plane_trap.m b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_in_plane_trap.m index 959c0ec..fe38444 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 @@ -1,5 +1,5 @@ %% Tilting of the dipoles -% Atom Number = 1250 ppum +% With an in-plane harmonic trap %% v_z = 500, theta = 0 @@ -17,7 +17,7 @@ OptionsStruct.NumberOfGridPoints = [256, 256]; OptionsStruct.Dimensions = [35, 35]; OptionsStruct.TimeStepSize = 0.005; % in s OptionsStruct.MinimumTimeStepSize = 1E-5; % in s -OptionsStruct.TimeCutOff = 2E6; % in s +OptionsStruct.TimeCutOff = 1E6; % in s OptionsStruct.EnergyTolerance = 5E-10; OptionsStruct.ResidualTolerance = 1E-05; OptionsStruct.NoiseScaleFactor = 0.05; @@ -59,7 +59,7 @@ OptionsStruct.NumberOfGridPoints = [256, 256]; OptionsStruct.Dimensions = [35, 35]; OptionsStruct.TimeStepSize = 0.005; % in s OptionsStruct.MinimumTimeStepSize = 1E-5; % in s -OptionsStruct.TimeCutOff = 2E6; % in s +OptionsStruct.TimeCutOff = 1E6; % in s OptionsStruct.EnergyTolerance = 5E-10; OptionsStruct.ResidualTolerance = 1E-05; OptionsStruct.NoiseScaleFactor = 0.05;