From 5a8b117e3e3345faf866b651cf5fbceda9c1d2c3 Mon Sep 17 00:00:00 2001 From: Karthik Date: Mon, 12 May 2025 11:56:45 +0200 Subject: [PATCH] More minor tweaks --- .../+Plotter/visualizeGSWavefunction.m | 3 +- Dipolar-Gas-Simulator/+Scripts/run_locally.m | 212 +++++++++--------- 2 files changed, 104 insertions(+), 111 deletions(-) diff --git a/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m b/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m index 0016a12..a7529a5 100644 --- a/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m +++ b/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m @@ -34,7 +34,8 @@ function visualizeGSWavefunction(folder_path, run_index) n = abs(psi).^2; %Plotting - figure(1); + fig = figure(1); + fig.WindowState = 'maximized'; clf set(gcf,'Position', [100, 100, 1600, 900]) t = tiledlayout(2, 3, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x3 grid diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m index ef14eef..09c0ff4 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -556,7 +556,32 @@ Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) SaveDirectory = './Results/Data_3D/PhaseDiagram/ImagTimePropagation/'; JobNumber = 0; Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) -%% +%% Generate lists + +% Set display format +format longG + +% Generate the lists +list1 = linspace(1E5, 5E6, 25); +list2 = linspace(80, 105, 25); + +% Convert to strings with no scientific notation +str_list1 = compose('%.0f', list1); +str_list2 = compose('%.2f', list2); + +% Join as space-separated strings +row1 = strjoin(str_list1', ' '); +row2 = strjoin(str_list2', ' '); + +% Display results +disp(row1) +disp(row2) + +%% Use space-separated floating-point/integer values +SCATTERING_LENGTH_RANGE = "[80.00 81.04 82.08 83.12 84.17 85.21 86.25 87.29 88.33 89.38 90.42 91.46 92.50 93.54 94.58 95.62 96.67 97.71 98.75 99.79 100.83 101.88 102.92 103.96 105.00]"; +NUM_ATOMS_LIST = "[100000 304167 508333 712500 916667 1120833 1325000 1529167 1733333 1937500 2141667 2345833 2550000 2754167 2958333 3162500 3366667 3570833 3775000 3979167 4183333 4387500 4591667 4795833 5000000]"; + +%% Explore the phase diagram SCATTERING_LENGTH_RANGE = [91.46 92.50 93.54 94.58 95.62 96.67 97.71]; NUM_ATOMS_LIST = [100000 304167 508333 712500 916667 1120833 1325000 1529167 1733333 1937500 2141667 2345833 2550000 2754167 2958333 3162500 3366667 3570833 3775000 3979167 4183333 4387500 4591667 4795833 5000000]; @@ -567,7 +592,7 @@ for i = 1:length(SCATTERING_LENGTH_RANGE) for j = 1:length(NUM_ATOMS_LIST) N = NUM_ATOMS_LIST(j); - SaveDirectory = sprintf('D:/Results/Data_3D/PhaseDiagram/aS_%.6e_theta_000_phi_000_N_%d', aS, N); + SaveDirectory = sprintf('D:/Results-Numerics/PhaseDiagram/ImagTimePropagation/aS_%.6e_theta_000_phi_000_N_%d', aS, N); fprintf('Processing JobNumber %d: %s\n', JobNumber, SaveDirectory); % Call the plotting function @@ -583,6 +608,81 @@ for i = 1:length(SCATTERING_LENGTH_RANGE) end end +%% Parameters that need to be run again +% Cell array of the input strings +strs = { + 'aS_9.875000e+01_theta_000_phi_000_N_2958333' + 'aS_9.562000e+01_theta_000_phi_000_N_4183333' + 'aS_9.562000e+01_theta_000_phi_000_N_1325000' + 'aS_9.458000e+01_theta_000_phi_000_N_508333' + 'aS_9.458000e+01_theta_000_phi_000_N_5000000' + 'aS_9.458000e+01_theta_000_phi_000_N_1937500' + 'aS_9.354000e+01_theta_000_phi_000_N_4795833' + 'aS_9.250000e+01_theta_000_phi_000_N_5000000' + 'aS_9.250000e+01_theta_000_phi_000_N_4183333' + 'aS_9.146000e+01_theta_000_phi_000_N_1937500' + 'aS_8.833000e+01_theta_000_phi_000_N_2754167' + 'aS_8.833000e+01_theta_000_phi_000_N_1325000' + 'aS_8.729000e+01_theta_000_phi_000_N_3775000' + 'aS_8.625000e+01_theta_000_phi_000_N_712500' + 'aS_8.625000e+01_theta_000_phi_000_N_4795833' + 'aS_8.625000e+01_theta_000_phi_000_N_3979167' + 'aS_8.521000e+01_theta_000_phi_000_N_1937500' + 'aS_8.208000e+01_theta_000_phi_000_N_3979167' + 'aS_105_theta_000_phi_000_N_916667' + 'aS_1.039600e+02_theta_000_phi_000_N_3366667' + 'aS_1.008300e+02_theta_000_phi_000_N_3979167' + 'aS_080_theta_000_phi_000_N_5000000' +}; + +% Extract values using regex +tokens = regexp(strs, 'aS_([0-9.e+-]+)_theta_000_phi_000_N_([0-9]+)', 'tokens'); + +% Convert to numeric array +pairs = cellfun(@(x) [str2double(x{1}), str2double(x{2})], [tokens{:}], 'UniformOutput', false); +pairs = vertcat(pairs{:}); + +% Sort by a_s (first column) +pairs = sortrows(pairs, 1); + +% Display as a table +disp(array2table(pairs, 'VariableNames', {'a_s', 'N'})); + +%% Phase diagram for untilted case +N = [1E5, 3.04E5, 5.08E5, 7.125E5, 9.16E5, 1.12E6, 1.325E6, 1.529E6, ... + 1.733E6, 1.9375E6, 2.141E6, 2.345E6, 2.55E6, 2.75E6, 2.95E6, ... + 3.1625E6, 3.367E6, 3.57E6, 3.775E6, 3.979E6, 4.183E6, 4.3875E6, ... + 4.591E6, 4.795E6, 5E6]; + +as_UB = [88.33, 94.58, 95.62, 96.67, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, ... + 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71]; + +as_LB = [87.29, 93.54, 94.58, 95.62, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, ... + 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67]; + +% Filter only rows with non-NaN data +valid_idx = ~isnan(as_LB) & ~isnan(as_UB); +N_valid = N(valid_idx); +LB_valid = as_LB(valid_idx); +UB_valid = as_UB(valid_idx); + +% Create shaded area between LB and UB +x_fill = [N_valid, fliplr(N_valid)]; +y_fill = [LB_valid, fliplr(UB_valid)]; + +% Plot settings +figure(1); +set(gcf,'Position',[100 100 950 750]) +fill(x_fill, y_fill, [0.8 0.8 1], 'EdgeColor', 'none'); % Light blue shade +% Axes settings +set(gca, 'XScale', 'log'); +xlim([1E4, 1E7]); +ylim([79, 106]); +xlabel('Atom number', 'Interpreter', 'latex', 'FontSize', 16); +ylabel('Scattering length $a_s$ ($a_0$)', 'Interpreter', 'latex', 'FontSize', 16); +grid on; +set(gca,'FontSize',16,'Box','On','Linewidth',2); + %% Visualize phase diagram load('./Results/Data_3D/GradientDescent/phase_diagram_matrix_theta_30.mat') PhaseDiagramMatrix = M; @@ -642,111 +742,3 @@ if ModulationFlag else disp('The state is not modulated'); end - -%% Generate lists - -% Set display format -format longG - -% Generate the lists -list1 = linspace(1E5, 5E6, 25); -list2 = linspace(80, 105, 25); - -% Convert to strings with no scientific notation -str_list1 = compose('%.0f', list1); -str_list2 = compose('%.2f', list2); - -% Join as space-separated strings -row1 = strjoin(str_list1', ' '); -row2 = strjoin(str_list2', ' '); - -% Display results -disp(row1) -disp(row2) - -%% Use space-separated floating-point/integer values -SCATTERING_LENGTH_RANGE = "[80.00 81.04 82.08 83.12 84.17 85.21 86.25 87.29 88.33 89.38 90.42 91.46 92.50 93.54 94.58 95.62 96.67 97.71 98.75 99.79 100.83 101.88 102.92 103.96 105.00]"; -NUM_ATOMS_LIST = "[100000 304167 508333 712500 916667 1120833 1325000 1529167 1733333 1937500 2141667 2345833 2550000 2754167 2958333 3162500 3366667 3570833 3775000 3979167 4183333 4387500 4591667 4795833 5000000]"; - -%% Phase diagram for untilted case -N = [1E5, 3.04E5, 5.08E5, 7.125E5, 9.16E5, 1.12E6, 1.325E6, 1.529E6, ... - 1.733E6, 1.9375E6, 2.141E6, 2.345E6, 2.55E6, 2.75E6, 2.95E6, ... - 3.1625E6, 3.367E6, 3.57E6, 3.775E6, 3.979E6, 4.183E6, 4.3875E6, ... - 4.591E6, 4.795E6, 5E6]; - -as_UB = [88.33, 94.58, 95.62, 96.67, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, ... - 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71]; - -as_LB = [87.29, 93.54, 94.58, 95.62, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, ... - 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67]; - -% Filter only rows with non-NaN data -valid_idx = ~isnan(as_LB) & ~isnan(as_UB); -N_valid = N(valid_idx); -LB_valid = as_LB(valid_idx); -UB_valid = as_UB(valid_idx); - -% Create shaded area between LB and UB -x_fill = [N_valid, fliplr(N_valid)]; -y_fill = [LB_valid, fliplr(UB_valid)]; - -% Plot settings -figure(1); -set(gcf,'Position',[100 100 950 750]) -fill(x_fill, y_fill, [0.8 0.8 1], 'EdgeColor', 'none'); % Light blue shade -% Axes settings -set(gca, 'XScale', 'log'); -xlim([1E4, 1E7]); -ylim([79, 106]); -xlabel('Atom number', 'Interpreter', 'latex', 'FontSize', 16); -ylabel('Scattering length $a_s$ ($a_0$)', 'Interpreter', 'latex', 'FontSize', 16); -grid on; -set(gca,'FontSize',16,'Box','On','Linewidth',2); - -%% - -function ModulationFlag = determineDensityModulation(psi, Params, Transf) - - % Axes scaling and coordinates in micrometers - x = Transf.x * Params.l0 * 1e6; - dx = x(2)-x(1); - % Compute probability density |psi|^2 - n = abs(psi).^2; - nyz = squeeze(trapz(n*dx,1)); - - densityProfile = sum(nyz,2); - - % We need to first smoothen the density profile and subtract the smoothened profile from the - % original density profile to - - % 1. Remove the dominant low-frequency content. - % 2. Isolate the high-frequency components (like a sinusoidal modulation). - % FFT of the residual then highlights the periodic part, making the - % modulation easy to detect. - - % Step 1: Smooth the density profile (Gaussian smoothing) - smoothedProfile = smooth(densityProfile, 10); - - % Step 2: Compute the residual (original - smoothed) - residual = densityProfile - smoothedProfile; - - % Step 3: Compute the Fourier Transform of the residual - N = length(residual); - Y = fft(residual); - P2 = abs(Y/N); % Two-sided spectrum - P1 = P2(1:N/2+1); % Single-sided spectrum - P1(2:end-1) = 2*P1(2:end-1); % Correct for the energy in the negative frequencies - - P1 = P1(15:end); - - % Step 4: Check for significant peaks in the Fourier spectrum - % We check if the peak frequency is above a certain threshold - threshold = 500; % This can be adjusted based on the expected modulation strength - peakValue = max(P1); - - if peakValue > threshold - ModulationFlag = true; % Indicates sinusoidal modulation - else - ModulationFlag = false; % Indicates otherwise - end -end -