From e7a4d3990670add67b21b6953c18dbdeb136ea1c Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Mon, 10 Feb 2025 22:46:33 +0100 Subject: [PATCH] Correct way to determine the wavenumber axes and calculate the lattice vectors now established, appropriate changes made elsewhere. --- .../+Scripts/analyzeGSWavefunction.m | 2 +- .../analyzeGSWavefunction_in_plane_trap.m | 118 ++++++++ .../+Scripts/extractLatticeProperties.m | 38 ++- Dipolar-Gas-Simulator/+Scripts/run_locally.m | 23 +- .../+Scripts/run_on_cluster_in_plane_trap.m | 264 +++++++++++++++++- .../run_on_cluster_transition_angle.m | 208 +++++++++++++- 6 files changed, 619 insertions(+), 34 deletions(-) create mode 100644 Dipolar-Gas-Simulator/+Scripts/analyzeGSWavefunction_in_plane_trap.m diff --git a/Dipolar-Gas-Simulator/+Scripts/analyzeGSWavefunction.m b/Dipolar-Gas-Simulator/+Scripts/analyzeGSWavefunction.m index b23c0f3..af1ba73 100644 --- a/Dipolar-Gas-Simulator/+Scripts/analyzeGSWavefunction.m +++ b/Dipolar-Gas-Simulator/+Scripts/analyzeGSWavefunction.m @@ -105,7 +105,7 @@ function [contrast, periodX, periodY] = analyzeGSWavefunction(folder_path, run_i cbar1.Label.Interpreter = 'latex'; colormap(gca, 'jet') hold on; - plot(freq_x*Params.l0, freq_y*Params.l0, 'ro', 'MarkerSize', 10, 'LineWidth', 2); + plot(2*pi*freq_x*Params.l0, 2*pi*freq_y*Params.l0, 'ro', 'MarkerSize', 10, 'LineWidth', 2); hold off; xlabel('$k_x l_o$', 'Interpreter', 'latex', 'FontSize', 14) ylabel('$k_y l_o$', 'Interpreter', 'latex', 'FontSize', 14) diff --git a/Dipolar-Gas-Simulator/+Scripts/analyzeGSWavefunction_in_plane_trap.m b/Dipolar-Gas-Simulator/+Scripts/analyzeGSWavefunction_in_plane_trap.m new file mode 100644 index 0000000..8eab17d --- /dev/null +++ b/Dipolar-Gas-Simulator/+Scripts/analyzeGSWavefunction_in_plane_trap.m @@ -0,0 +1,118 @@ +function [contrast, periodX, periodY] = analyzeGSWavefunction(folder_path, run_index) + + set(0,'defaulttextInterpreter','latex') + set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex'); + + % Load data + Data = load(sprintf(horzcat(folder_path, '/Run_%03i/psi_gs.mat'),run_index),'psi','Transf','Observ','Params','VParams'); + + Transf = Data.Transf; + Observ = Data.Observ; + Params = Data.Params; + VParams = Data.VParams; + + if isgpuarray(Data.psi) + psi = gather(Data.psi); + else + psi = Data.psi; + end + + if isgpuarray(Data.Observ.residual) + Observ.residual = gather(Data.Observ.residual); + else + Observ.residual = Data.Observ.residual; + end + + % Set long format for output + format long + + % Coordinates in micrometers + x = Transf.x * Params.l0; + y = Transf.y * Params.l0; + + % Compute probability density |psi|^2 + nxy = abs(psi).^2; + nxyScaled = nxy*(Params.add*10^6)^2; + + %------------------ Calculate contrast ------------------ % + + % Find unique maximum intensity values + uniqueMaxValues = unique(nxyScaled(nxyScaled == max(nxyScaled(:)))); + if length(uniqueMaxValues) > 1 + maxIntensity = median(uniqueMaxValues); % Choose the median of max values + else + maxIntensity = uniqueMaxValues; % If only one, take the unique value + end + + % Find unique minimum intensity values + uniqueMinValues = unique(nxyScaled(nxyScaled == min(nxyScaled(:)))); + if length(uniqueMinValues) > 1 + minIntensity = median(uniqueMinValues); % Choose the median of min values + else + minIntensity = uniqueMinValues; % If only one, take the unique value + end + + contrast = (maxIntensity - minIntensity) / (maxIntensity + minIntensity); + + %------------------ Lattice Properties ------------------ % + + [kx, ky, fftMagnitude, lattice_type, periodX, periodY, freq_x, freq_y] = Scripts.extractLatticeProperties(nxyScaled, x, y); + + %------------------ Plotting ------------------ % + + figure('Position', [100, 100, 1600, 800]); + clf + t = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x3 grid + + % Plot |psi(x,y)|^2 (density in x and y plane) + nexttile; % Equivalent to subplot('Position', [0.05, 0.55, 0.28, 0.4]); + plotxy = pcolor(x/Params.l0,y/Params.l0,nxyScaled'); + set(plotxy, 'EdgeColor', 'none'); + cbar1 = colorbar; + 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$ ', 'Interpreter', 'latex', 'FontSize', 14) + + % Plot 2-D FFT with detected peaks + nexttile; % Equivalent to subplot('Position', [0.05, 0.55, 0.28, 0.4]); + imagesc(kx*Params.l0,ky*Params.l0, log(1 + fftMagnitude)); + cbar1 = colorbar; + cbar1.Label.Interpreter = 'latex'; + colormap(gca, 'jet') + %{ + hold on; + plot(freq_x*Params.l0, freq_y*Params.l0, 'ro', 'MarkerSize', 10, 'LineWidth', 2); + hold off; + %} + xlabel('$k_x l_o$', 'Interpreter', 'latex', 'FontSize', 14) + ylabel('$k_y l_o$', 'Interpreter', 'latex', 'FontSize', 14) + title('$\mathcal{F}\{|\Psi(x,y)|^2\}$', 'Interpreter', 'latex', 'FontSize', 16); + + sgtitle(['$[\omega_x, \omega_y, \omega_z] = 2 \pi~\times$ [', num2str(round(Params.wx/(2*pi)), '%.2f'), ', ', num2str(round(Params.wy/(2*pi)), '%.2f'), ', ', 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 a007a6f..74c90c3 100644 --- a/Dipolar-Gas-Simulator/+Scripts/extractLatticeProperties.m +++ b/Dipolar-Gas-Simulator/+Scripts/extractLatticeProperties.m @@ -18,18 +18,30 @@ function [kx, ky, fftMagnitude, lattice_type, dx_um, dy_um, freq_x, freq_y] = ex F = fft2(double(I)); fftMagnitude = abs(fftshift(F))'; % Shift zero frequency to center - % Calculate wavenumber increment (frequency axes) - Nx = length(x); - Ny = length(y); - dx = mean(diff(x)); - dy = mean(diff(y)); - dkx = 1 / (Nx * dx); % Increment in the X direction (in micrometers^-1) - dky = 1 / (Ny * dy); % Increment in the Y direction (in micrometers^-1) + % Calculate frequency increment (frequency axes) + Nx = length(x); % grid size along X + Ny = length(y); % grid size along Y + dx = mean(diff(x)); % real space increment in the X direction (in micrometers) + dy = mean(diff(y)); % real space increment in the Y direction (in micrometers) + dvx = 1 / (Nx * dx); % reciprocal space increment in the X direction (in micrometers^-1) + dvy = 1 / (Ny * dy); % reciprocal space increment in the Y direction (in micrometers^-1) - % Create the wavenumber axes - kx = (-Nx/2:Nx/2-1) * dkx; % Wavenumber axis in X (micrometers^-1) - ky = (-Ny/2:Ny/2-1) * dky; % Wavenumber axis in Y (micrometers^-1) + % Create the frequency axes + vx = (-Nx/2:Nx/2-1) * dvx; % Frequency axis in X (micrometers^-1) + vy = (-Ny/2:Ny/2-1) * dvy; % Frequency axis in Y (micrometers^-1) + + % Calculate maximum frequencies + % kx_max = pi / dx; + % ky_max = pi / dy; + + % Generate reciprocal axes + % kx = linspace(-kx_max, kx_max * (Nx-2)/Nx, Nx); + % ky = linspace(-ky_max, ky_max * (Ny-2)/Ny, Ny); + % Create the Wavenumber axes + kx = 2*pi*vx; % Wavenumber axis in X + ky = 2*pi*vy; % Wavenumber axis in Y + % Find peaks in Fourier domain peak_mask = imregionalmax(fftMagnitude); % Detect local maxima threshold = max(fftMagnitude(:)) * 0.1; % Adaptive threshold @@ -39,8 +51,8 @@ function [kx, ky, fftMagnitude, lattice_type, dx_um, dy_um, freq_x, freq_y] = ex [rows, cols] = find(peak_mask); % Convert indices to frequency values - freq_x = kx(cols); - freq_y = ky(rows); + freq_x = vx(cols); + freq_y = vy(rows); distances = sqrt(freq_x.^2 + freq_y.^2); % Remove DC component (center peak) @@ -128,7 +140,7 @@ function [kx, ky, fftMagnitude, lattice_type, dx_um, dy_um, freq_x, freq_y] = ex else % If reciprocal vectors are not collinear, compute the 2D real-space lattice - real_lattice_vectors = inv(reciprocal_lattice_vectors'); + real_lattice_vectors = inv(reciprocal_lattice_vectors'); % If vx, vy are used, no need to multiply by 2*pi (crystallographic convention); If kx, ky are used, need to multiply by 2*pi (physics convention) end % Ensure correct orientation diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m index 6263b9a..a2b4a5b 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -251,53 +251,58 @@ Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber) %% - Analysis SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz500'; -JobNumber = 2; +JobNumber = 1; % Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber) [contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction(SaveDirectory, JobNumber); %% - Analysis SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz750'; -JobNumber = 2; +JobNumber = 1; % Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber) [contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction(SaveDirectory, JobNumber); %% - Analysis SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz1000'; -JobNumber = 2; +JobNumber = 1; % Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber) [contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction(SaveDirectory, JobNumber); %% - Analysis SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz2000'; -JobNumber = 2; +JobNumber = 1; % Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber) [contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction(SaveDirectory, JobNumber); %% - Analysis SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/Hz500'; NumberOfJobs = 16; -SaveOption = 'images'; +SaveOption = 'video'; [contrast_array, periodX_array, periodY_array] = Scripts.analyzeGSWavefunction_multipleruns(SaveDirectory, NumberOfJobs, SaveOption); %% - Analysis SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/Hz750'; NumberOfJobs = 16; -SaveOption = 'images'; +SaveOption = 'video'; [contrast_array, periodX_array, periodY_array] = Scripts.analyzeGSWavefunction_multipleruns(SaveDirectory, NumberOfJobs, SaveOption); %% - Analysis SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/Hz1000'; NumberOfJobs = 16; -SaveOption = 'images'; +SaveOption = 'video'; [contrast_array, periodX_array, periodY_array] = Scripts.analyzeGSWavefunction_multipleruns(SaveDirectory, NumberOfJobs, SaveOption); %% - Analysis SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/Hz2000'; NumberOfJobs = 16; -SaveOption = 'images'; +SaveOption = 'video'; [contrast_array, periodX_array, periodY_array] = Scripts.analyzeGSWavefunction_multipleruns(SaveDirectory, NumberOfJobs, SaveOption); %% - Analysis +SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle_newruns/Hz1000'; +JobNumber = 0; +% Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber) +[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction(SaveDirectory, JobNumber); +%% - Analysis SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/Hz500'; JobNumber = 1; % Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber) -[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction(SaveDirectory, JobNumber); +[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber); 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 fe38444..25df0c5 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 @@ -10,20 +10,20 @@ OptionsStruct.DipolarPolarAngle = 0; OptionsStruct.DipolarAzimuthAngle = 0; OptionsStruct.ScatteringLength = 70.00; -OptionsStruct.TrapFrequencies = [50, 50, 500]; +OptionsStruct.TrapFrequencies = [100, 100, 500]; OptionsStruct.TrapPotentialType = 'Harmonic'; OptionsStruct.NumberOfGridPoints = [256, 256]; OptionsStruct.Dimensions = [35, 35]; OptionsStruct.TimeStepSize = 0.005; % in s OptionsStruct.MinimumTimeStepSize = 1E-5; % in s -OptionsStruct.TimeCutOff = 1E6; % in s +OptionsStruct.TimeCutOff = 2E6; % in s OptionsStruct.EnergyTolerance = 5E-10; OptionsStruct.ResidualTolerance = 1E-05; OptionsStruct.NoiseScaleFactor = 0.05; OptionsStruct.MaxIterations = 10; -OptionsStruct.VariationalWidth = 1.5; +OptionsStruct.VariationalWidth = 1.2; OptionsStruct.WidthLowerBound = 0.01; OptionsStruct.WidthUpperBound = 12; OptionsStruct.WidthCutoff = 5e-3; @@ -52,20 +52,20 @@ OptionsStruct.DipolarPolarAngle = deg2rad(15); OptionsStruct.DipolarAzimuthAngle = 0; OptionsStruct.ScatteringLength = 71.00; -OptionsStruct.TrapFrequencies = [50, 50, 500]; +OptionsStruct.TrapFrequencies = [100, 100, 500]; OptionsStruct.TrapPotentialType = 'Harmonic'; OptionsStruct.NumberOfGridPoints = [256, 256]; OptionsStruct.Dimensions = [35, 35]; OptionsStruct.TimeStepSize = 0.005; % in s OptionsStruct.MinimumTimeStepSize = 1E-5; % in s -OptionsStruct.TimeCutOff = 1E6; % in s +OptionsStruct.TimeCutOff = 2E6; % in s OptionsStruct.EnergyTolerance = 5E-10; OptionsStruct.ResidualTolerance = 1E-05; OptionsStruct.NoiseScaleFactor = 0.05; OptionsStruct.MaxIterations = 10; -OptionsStruct.VariationalWidth = 1.5; +OptionsStruct.VariationalWidth = 1.2; OptionsStruct.WidthLowerBound = 0.01; OptionsStruct.WidthUpperBound = 12; OptionsStruct.WidthCutoff = 5e-3; @@ -82,5 +82,257 @@ solver = VariationalSolver2D.DipolarGas(options{: pot = VariationalSolver2D.Potentials(options{:}); solver.Potential = pot.trap(); +%-% Run Solver %-% +[Params, Transf, psi, V, VDk] = solver.run(); + +%% v_z = 750, theta = 0 + +OptionsStruct = struct; + +OptionsStruct.NumberOfAtoms = 61250; +OptionsStruct.DipolarPolarAngle = 0; +OptionsStruct.DipolarAzimuthAngle = 0; +OptionsStruct.ScatteringLength = 64.0; + +OptionsStruct.TrapFrequencies = [100, 100, 750]; +OptionsStruct.TrapPotentialType = 'Harmonic'; + +OptionsStruct.NumberOfGridPoints = [256, 256]; +OptionsStruct.Dimensions = [35, 35]; +OptionsStruct.TimeStepSize = 0.005; % in s +OptionsStruct.MinimumTimeStepSize = 1E-5; % in s +OptionsStruct.TimeCutOff = 1E5; % in s +OptionsStruct.EnergyTolerance = 5E-10; +OptionsStruct.ResidualTolerance = 1E-05; +OptionsStruct.NoiseScaleFactor = 0.05; + +OptionsStruct.MaxIterations = 10; +OptionsStruct.VariationalWidth = 0.8; +OptionsStruct.WidthLowerBound = 0.01; +OptionsStruct.WidthUpperBound = 12; +OptionsStruct.WidthCutoff = 5e-3; + +OptionsStruct.PlotLive = false; +OptionsStruct.JobNumber = 0; +OptionsStruct.RunOnGPU = true; +OptionsStruct.SaveData = true; +OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/Hz750'; +options = Helper.convertstruct2cell(OptionsStruct); +clear OptionsStruct + +solver = VariationalSolver2D.DipolarGas(options{:}); +pot = VariationalSolver2D.Potentials(options{:}); +solver.Potential = pot.trap(); + +%-% Run Solver %-% +[Params, Transf, psi, V, VDk] = solver.run(); + +%% v_z = 750, theta = 15 + +OptionsStruct = struct; + +OptionsStruct.NumberOfAtoms = 61250; +OptionsStruct.DipolarPolarAngle = deg2rad(15); +OptionsStruct.DipolarAzimuthAngle = 0; +OptionsStruct.ScatteringLength = 64.0; + +OptionsStruct.TrapFrequencies = [100, 100, 750]; +OptionsStruct.TrapPotentialType = 'Harmonic'; + +OptionsStruct.NumberOfGridPoints = [256, 256]; +OptionsStruct.Dimensions = [35, 35]; +OptionsStruct.TimeStepSize = 0.005; % in s +OptionsStruct.MinimumTimeStepSize = 1E-5; % in s +OptionsStruct.TimeCutOff = 1E5; % in s +OptionsStruct.EnergyTolerance = 5E-10; +OptionsStruct.ResidualTolerance = 1E-05; +OptionsStruct.NoiseScaleFactor = 0.05; + +OptionsStruct.MaxIterations = 10; +OptionsStruct.VariationalWidth = 0.8; +OptionsStruct.WidthLowerBound = 0.01; +OptionsStruct.WidthUpperBound = 12; +OptionsStruct.WidthCutoff = 5e-3; + +OptionsStruct.PlotLive = false; +OptionsStruct.JobNumber = 1; +OptionsStruct.RunOnGPU = true; +OptionsStruct.SaveData = true; +OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/Hz750'; +options = Helper.convertstruct2cell(OptionsStruct); +clear OptionsStruct + +solver = VariationalSolver2D.DipolarGas(options{:}); +pot = VariationalSolver2D.Potentials(options{:}); +solver.Potential = pot.trap(); + +%-% Run Solver %-% +[Params, Transf, psi, V, VDk] = solver.run(); + +%% v_z = 1000, theta = 0 + +OptionsStruct = struct; + +OptionsStruct.NumberOfAtoms = 45000; +OptionsStruct.DipolarPolarAngle = 0; +OptionsStruct.DipolarAzimuthAngle = 0; +OptionsStruct.ScatteringLength = 59.0; + +OptionsStruct.TrapFrequencies = [100, 100, 1000]; +OptionsStruct.TrapPotentialType = 'Harmonic'; + +OptionsStruct.NumberOfGridPoints = [256, 256]; +OptionsStruct.Dimensions = [35, 35]; +OptionsStruct.TimeStepSize = 0.005; % in s +OptionsStruct.MinimumTimeStepSize = 1E-5; % in s +OptionsStruct.TimeCutOff = 1E5; % in s +OptionsStruct.EnergyTolerance = 5E-10; +OptionsStruct.ResidualTolerance = 1E-05; +OptionsStruct.NoiseScaleFactor = 0.05; + +OptionsStruct.MaxIterations = 10; +OptionsStruct.VariationalWidth = 0.7; +OptionsStruct.WidthLowerBound = 0.01; +OptionsStruct.WidthUpperBound = 12; +OptionsStruct.WidthCutoff = 5e-3; + +OptionsStruct.PlotLive = false; +OptionsStruct.JobNumber = 0; +OptionsStruct.RunOnGPU = true; +OptionsStruct.SaveData = true; +OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/Hz1000'; +options = Helper.convertstruct2cell(OptionsStruct); +clear OptionsStruct + +solver = VariationalSolver2D.DipolarGas(options{:}); +pot = VariationalSolver2D.Potentials(options{:}); +solver.Potential = pot.trap(); + +%-% Run Solver %-% +[Params, Transf, psi, V, VDk] = solver.run(); + +%% v_z = 1000, theta = 15 + +OptionsStruct = struct; + +OptionsStruct.NumberOfAtoms = 45000; +OptionsStruct.DipolarPolarAngle = deg2rad(15); +OptionsStruct.DipolarAzimuthAngle = 0; +OptionsStruct.ScatteringLength = 60.0; + +OptionsStruct.TrapFrequencies = [100, 100, 1000]; +OptionsStruct.TrapPotentialType = 'Harmonic'; + +OptionsStruct.NumberOfGridPoints = [256, 256]; +OptionsStruct.Dimensions = [35, 35]; +OptionsStruct.TimeStepSize = 0.005; % in s +OptionsStruct.MinimumTimeStepSize = 1E-5; % in s +OptionsStruct.TimeCutOff = 1E5; % in s +OptionsStruct.EnergyTolerance = 5E-10; +OptionsStruct.ResidualTolerance = 1E-05; +OptionsStruct.NoiseScaleFactor = 0.05; + +OptionsStruct.MaxIterations = 10; +OptionsStruct.VariationalWidth = 0.7; +OptionsStruct.WidthLowerBound = 0.01; +OptionsStruct.WidthUpperBound = 12; +OptionsStruct.WidthCutoff = 5e-3; + +OptionsStruct.PlotLive = false; +OptionsStruct.JobNumber = 1; +OptionsStruct.RunOnGPU = true; +OptionsStruct.SaveData = true; +OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/Hz1000'; +options = Helper.convertstruct2cell(OptionsStruct); +clear OptionsStruct + +solver = VariationalSolver2D.DipolarGas(options{:}); +pot = VariationalSolver2D.Potentials(options{:}); +solver.Potential = pot.trap(); + +%-% Run Solver %-% +[Params, Transf, psi, V, VDk] = solver.run(); + +%% v_z = 2000, theta = 0 + +OptionsStruct = struct; + +OptionsStruct.NumberOfAtoms = 31250; +OptionsStruct.DipolarPolarAngle = 0; +OptionsStruct.DipolarAzimuthAngle = 0; +OptionsStruct.ScatteringLength = 47; + +OptionsStruct.TrapFrequencies = [100, 100, 2000]; +OptionsStruct.TrapPotentialType = 'Harmonic'; + +OptionsStruct.NumberOfGridPoints = [256, 256]; +OptionsStruct.Dimensions = [35, 35]; +OptionsStruct.TimeStepSize = 0.005; % in s +OptionsStruct.MinimumTimeStepSize = 1E-5; % in s +OptionsStruct.TimeCutOff = 1E5; % in s +OptionsStruct.EnergyTolerance = 5E-10; +OptionsStruct.ResidualTolerance = 1E-05; +OptionsStruct.NoiseScaleFactor = 0.05; + +OptionsStruct.MaxIterations = 10; +OptionsStruct.VariationalWidth = 0.5; +OptionsStruct.WidthLowerBound = 0.01; +OptionsStruct.WidthUpperBound = 12; +OptionsStruct.WidthCutoff = 5e-3; + +OptionsStruct.PlotLive = false; +OptionsStruct.JobNumber = 0; +OptionsStruct.RunOnGPU = true; +OptionsStruct.SaveData = true; +OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/Hz2000'; +options = Helper.convertstruct2cell(OptionsStruct); +clear OptionsStruct + +solver = VariationalSolver2D.DipolarGas(options{:}); +pot = VariationalSolver2D.Potentials(options{:}); +solver.Potential = pot.trap(); + +%-% Run Solver %-% +[Params, Transf, psi, V, VDk] = solver.run(); + +%% v_z = 2000, theta = 15 + +OptionsStruct = struct; + +OptionsStruct.NumberOfAtoms = 31250; +OptionsStruct.DipolarPolarAngle = deg2rad(15); +OptionsStruct.DipolarAzimuthAngle = 0; +OptionsStruct.ScatteringLength = 48; + +OptionsStruct.TrapFrequencies = [100, 100, 2000]; +OptionsStruct.TrapPotentialType = 'Harmonic'; + +OptionsStruct.NumberOfGridPoints = [256, 256]; +OptionsStruct.Dimensions = [35, 35]; +OptionsStruct.TimeStepSize = 0.005; % in s +OptionsStruct.MinimumTimeStepSize = 1E-5; % in s +OptionsStruct.TimeCutOff = 1E5; % in s +OptionsStruct.EnergyTolerance = 5E-10; +OptionsStruct.ResidualTolerance = 1E-05; +OptionsStruct.NoiseScaleFactor = 0.05; + +OptionsStruct.MaxIterations = 10; +OptionsStruct.VariationalWidth = 0.5; +OptionsStruct.WidthLowerBound = 0.01; +OptionsStruct.WidthUpperBound = 12; +OptionsStruct.WidthCutoff = 5e-3; + +OptionsStruct.PlotLive = false; +OptionsStruct.JobNumber = 1; +OptionsStruct.RunOnGPU = true; +OptionsStruct.SaveData = true; +OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/Hz2000'; +options = Helper.convertstruct2cell(OptionsStruct); +clear OptionsStruct + +solver = VariationalSolver2D.DipolarGas(options{:}); +pot = VariationalSolver2D.Potentials(options{:}); +solver.Potential = pot.trap(); + %-% Run Solver %-% [Params, Transf, psi, V, VDk] = solver.run(); \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_transition_angle.m b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_transition_angle.m index 086a46a..25f4c30 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_transition_angle.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_transition_angle.m @@ -2,11 +2,11 @@ % Atom Number = 1250 ppum % System size = [5 * l_rot, 5 * l_rot] -theta_values = 1:1:14; +theta_values = 0:1:7; %% v_z = 500 -as_values_500 = [76.41, 76.54, 76.54, 76.54, 76.54, 76.67, 76.67, 76.80, 76.80, 76.93, 77.06, 77.06, 77.19, 77.32]; +as_values_500 = [76.408020791433103, 76.408020791433103, 76.538777583310534, 76.538777583310534, 76.538777583310534, 76.538777583310534, 76.669534375187951, 76.669534375187951]; num_iterations = length(theta_values); for i = 1:num_iterations @@ -54,7 +54,7 @@ end %% v_z = 750 -as_values_750 = [70.52, 70.52, 70.52, 70.52, 70.52, 70.65, 70.65, 70.78, 70.79, 70.92, 71.05, 71.05, 71.18, 71.31, 71.44, 71.70]; +as_values_750 = [70.523965156949174, 70.523965156949174, 70.523965156949174, 70.523965156949174, 70.523965156949174, 70.654721948826605, 70.654721948826605, 70.785478740704022]; num_iterations = length(theta_values); for i = 1:num_iterations @@ -102,7 +102,7 @@ end %% v_z = 1000 -as_values_1000 = [65.95, 65.95, 65.95, 65.95, 66.08, 66.08, 66.08, 66.21, 66.34, 66.34, 66.47, 66.60, 66.73, 66.86, 66.99, 67.26]; +as_values_1000 = [65.947477441239442, 65.947477441239442, 65.947477441239442, 65.947477441239442, 66.078234233116873, 66.078234233116873, 66.078234233116873, 66.208991024994290]; num_iterations = length(theta_values); for i = 1:num_iterations @@ -150,7 +150,7 @@ end %% v_z = 2000 -as_values_2000 = [53.92, 53.92, 53.92, 54.05, 54.05, 54.05, 54.18, 54.31, 54.31, 54.44, 54.57, 54.70, 54.96, 55.1, 55.23, 55.49]; +as_values_2000 = [53.917852588516730, 53.917852588516730, 53.917852588516730, 54.048609380394161, 54.048609380394161, 54.048609380394161, 54.179366172271585, 54.310122964149002]; num_iterations = length(theta_values); for i = 1:num_iterations @@ -192,6 +192,204 @@ for i = 1:num_iterations pot = VariationalSolver2D.Potentials(options{:}); solver.Potential = pot.trap(); + % Run Solver + [Params, Transf, psi, V, VDk] = solver.run(); +end + +%% Tilting of the dipoles +% Atom Number = 1250 ppum +% System size = [5 * l_rot, 5 * l_rot] + +theta_values = 8:1:15; + +%% v_z = 500 + +as_values_500 = [76.800291167065367, 76.800291167065367, 76.931047958942784, 77.061804750820215, 77.061804750820215, 77.192561542697632, 77.323318334575049, 77.454075126452480]; +num_iterations = length(theta_values); + +for i = 1:num_iterations + OptionsStruct = struct; + + OptionsStruct.NumberOfAtoms = 101250; + OptionsStruct.DipolarPolarAngle = deg2rad(theta_values(i)); + OptionsStruct.DipolarAzimuthAngle = 0; + OptionsStruct.ScatteringLength = as_values_500(i); + + OptionsStruct.TrapFrequencies = [0, 0, 500]; + OptionsStruct.TrapPotentialType = 'None'; + + OptionsStruct.NumberOfGridPoints = [128, 128]; + OptionsStruct.Dimensions = [9, 9]; + OptionsStruct.TimeStepSize = 0.005; % in s + OptionsStruct.MinimumTimeStepSize = 1E-5; % in s + OptionsStruct.TimeCutOff = 3E6; % in s + OptionsStruct.EnergyTolerance = 5E-10; + OptionsStruct.ResidualTolerance = 1E-05; + OptionsStruct.NoiseScaleFactor = 0.05; + + OptionsStruct.MaxIterations = 10; + OptionsStruct.VariationalWidth = 1.2; + OptionsStruct.WidthLowerBound = 0.01; + OptionsStruct.WidthUpperBound = 12; + OptionsStruct.WidthCutoff = 5e-3; + + OptionsStruct.PlotLive = false; + OptionsStruct.JobNumber = i-1; % Assign a unique JobNumber per iteration + OptionsStruct.RunOnGPU = true; + OptionsStruct.SaveData = true; + OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/Hz500'; + + options = Helper.convertstruct2cell(OptionsStruct); + clear OptionsStruct + + solver = VariationalSolver2D.DipolarGas(options{:}); + pot = VariationalSolver2D.Potentials(options{:}); + solver.Potential = pot.trap(); + + % Run Solver + [Params, Transf, psi, V, VDk] = solver.run(); +end + +%% v_z = 750 + +as_values_750 = [70.785478740704022, 70.916235532581439, 71.046992324458870, 71.046992324458870, 71.177749116336287, 71.308505908213704, 71.439262700091120, 71.700776283845968]; +num_iterations = length(theta_values); + +for i = 1:num_iterations + OptionsStruct = struct; + + OptionsStruct.NumberOfAtoms = 61250; + OptionsStruct.DipolarPolarAngle = deg2rad(theta_values(i)); + OptionsStruct.DipolarAzimuthAngle = 0; + OptionsStruct.ScatteringLength = as_values_750(i); + + OptionsStruct.TrapFrequencies = [0, 0, 750]; + OptionsStruct.TrapPotentialType = 'None'; + + OptionsStruct.NumberOfGridPoints = [128, 128]; + OptionsStruct.Dimensions = [7, 7]; + OptionsStruct.TimeStepSize = 0.005; % in s + OptionsStruct.MinimumTimeStepSize = 1E-5; % in s + OptionsStruct.TimeCutOff = 3E6; % in s + OptionsStruct.EnergyTolerance = 5E-10; + OptionsStruct.ResidualTolerance = 1E-05; + OptionsStruct.NoiseScaleFactor = 0.05; + + OptionsStruct.MaxIterations = 10; + OptionsStruct.VariationalWidth = 0.85; + OptionsStruct.WidthLowerBound = 0.01; + OptionsStruct.WidthUpperBound = 12; + OptionsStruct.WidthCutoff = 5e-3; + + OptionsStruct.PlotLive = false; + OptionsStruct.JobNumber = i-1; % Assign a unique JobNumber per iteration + OptionsStruct.RunOnGPU = true; + OptionsStruct.SaveData = true; + OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/Hz750'; + + options = Helper.convertstruct2cell(OptionsStruct); + clear OptionsStruct + + solver = VariationalSolver2D.DipolarGas(options{:}); + pot = VariationalSolver2D.Potentials(options{:}); + solver.Potential = pot.trap(); + + % Run Solver + [Params, Transf, psi, V, VDk] = solver.run(); +end + +%% v_z = 1000 + +as_values_1000 = [66.339747816871707, 66.339747816871707, 66.470504608749138, 66.601261400626555, 66.732018192503972, 66.862774984381389, 66.993531776258820, 67.255045360013654]; +num_iterations = length(theta_values); + +for i = 1:num_iterations + OptionsStruct = struct; + + OptionsStruct.NumberOfAtoms = 45000; + OptionsStruct.DipolarPolarAngle = deg2rad(theta_values(i)); + OptionsStruct.DipolarAzimuthAngle = 0; + OptionsStruct.ScatteringLength = as_values_1000(i); + + OptionsStruct.TrapFrequencies = [0, 0, 1000]; + OptionsStruct.TrapPotentialType = 'None'; + + OptionsStruct.NumberOfGridPoints = [128, 128]; + OptionsStruct.Dimensions = [6, 6]; + OptionsStruct.TimeStepSize = 0.005; % in s + OptionsStruct.MinimumTimeStepSize = 1E-5; % in s + OptionsStruct.TimeCutOff = 3E6; % in s + OptionsStruct.EnergyTolerance = 5E-10; + OptionsStruct.ResidualTolerance = 1E-05; + OptionsStruct.NoiseScaleFactor = 0.05; + + OptionsStruct.MaxIterations = 10; + OptionsStruct.VariationalWidth = 0.7; + OptionsStruct.WidthLowerBound = 0.01; + OptionsStruct.WidthUpperBound = 12; + OptionsStruct.WidthCutoff = 5e-3; + + OptionsStruct.PlotLive = false; + OptionsStruct.JobNumber = i-1; % Assign a unique JobNumber per iteration + OptionsStruct.RunOnGPU = true; + OptionsStruct.SaveData = true; + OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/Hz1000'; + + options = Helper.convertstruct2cell(OptionsStruct); + clear OptionsStruct + + solver = VariationalSolver2D.DipolarGas(options{:}); + pot = VariationalSolver2D.Potentials(options{:}); + solver.Potential = pot.trap(); + + % Run Solver + [Params, Transf, psi, V, VDk] = solver.run(); +end + +%% v_z = 2000 + +as_values_2000 = [54.310122964149002, 54.440879756026426, 54.571636547903843, 54.702393339781267, 54.963906923536108, 55.094663715413525, 55.225420507290949, 55.486934091045789]; +num_iterations = length(theta_values); + +for i = 1:num_iterations + OptionsStruct = struct; + + OptionsStruct.NumberOfAtoms = 31250; + OptionsStruct.DipolarPolarAngle = deg2rad(theta_values(i)); + OptionsStruct.DipolarAzimuthAngle = 0; + OptionsStruct.ScatteringLength = as_values_2000(i); + + OptionsStruct.TrapFrequencies = [0, 0, 2000]; + OptionsStruct.TrapPotentialType = 'None'; + + OptionsStruct.NumberOfGridPoints = [128, 128]; + OptionsStruct.Dimensions = [5, 5]; + OptionsStruct.TimeStepSize = 0.005; % in s + OptionsStruct.MinimumTimeStepSize = 1E-5; % in s + OptionsStruct.TimeCutOff = 3E6; % in s + OptionsStruct.EnergyTolerance = 5E-10; + OptionsStruct.ResidualTolerance = 1E-05; + OptionsStruct.NoiseScaleFactor = 0.05; + + OptionsStruct.MaxIterations = 10; + OptionsStruct.VariationalWidth = 0.5; + OptionsStruct.WidthLowerBound = 0.01; + OptionsStruct.WidthUpperBound = 12; + OptionsStruct.WidthCutoff = 5e-3; + + OptionsStruct.PlotLive = false; + OptionsStruct.JobNumber = i-1; % Assign a unique JobNumber per iteration + OptionsStruct.RunOnGPU = true; + OptionsStruct.SaveData = true; + OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/Hz2000'; + + options = Helper.convertstruct2cell(OptionsStruct); + clear OptionsStruct + + solver = VariationalSolver2D.DipolarGas(options{:}); + pot = VariationalSolver2D.Potentials(options{:}); + solver.Potential = pot.trap(); + % Run Solver [Params, Transf, psi, V, VDk] = solver.run(); end \ No newline at end of file