diff --git a/.gitignore b/.gitignore index 69886cb..ab48627 100644 --- a/.gitignore +++ b/.gitignore @@ -9,5 +9,6 @@ Time-Series-Analyzer/Time-Series-Data *.gif *.mp4 *.bat +*.json .ipynb_checkpoints/ .vscode/ \ No newline at end of file diff --git a/Data-Analyzer/compareSpectralWeights.m b/Data-Analyzer/compareSpectralWeights.m new file mode 100644 index 0000000..e3ffb99 --- /dev/null +++ b/Data-Analyzer/compareSpectralWeights.m @@ -0,0 +1,75 @@ +%% Track spectral weight across the transition + +set(0,'defaulttextInterpreter','latex') +set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex'); + +format long + +% Load data +Data = load('C:/Users/Karthik/Documents/GitRepositories/Calculations/Data-Analyzer/B2.3G/WithoutProcessing/DropletsToStripes.mat', 'unique_theta', 'mean_sf', 'stderr_sf'); + +down_scan_parameter_values = Data.unique_theta; +dts_mean_sf = Data.mean_sf; +down_stderr_sf = Data.stderr_sf; + +Data = load('C:/Users/Karthik/Documents/GitRepositories/Calculations/Data-Analyzer/B2.3G/WithoutProcessing/StripesToDroplets.mat', 'unique_theta', 'mean_sf', 'stderr_sf'); + +std_theta_values = Data.unique_theta; +std_mean_sf = Data.mean_sf; +std_stderr_sf = Data.stderr_sf; + +figure(1); +set(gcf,'Position',[100 100 950 750]) +errorbar(down_scan_parameter_values, dts_mean_sf, down_stderr_sf, 'o--', ... + 'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5, 'DisplayName' , 'Droplets to Stripes'); +hold on +errorbar(std_theta_values, std_mean_sf, std_stderr_sf, 'o--', ... + 'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5, 'DisplayName', 'Stripes to Droplets'); +set(gca, 'FontSize', 14); % For tick labels only +hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex'); +hYLabel = ylabel('Spectral Weight', 'Interpreter', 'tex'); +hTitle = title('B = 2.3 G', 'Interpreter', 'tex'); +legend +set([hXLabel, hYLabel], 'FontName', font) +set([hXLabel, hYLabel], 'FontSize', 14) +set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title +grid on + +%% + +%% Track spectral weight across the transition + +set(0,'defaulttextInterpreter','latex') +set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex'); + +format long + +% Load data +Data = load('C:/Users/Karthik/Documents/GitRepositories/Calculations/Data-Analyzer/RampDownSL.mat', 'unique_scan_parameter_values', 'mean_sf', 'stderr_sf'); + +down_scan_parameter_values = Data.unique_scan_parameter_values; +down_mean_sf = Data.mean_sf; +down_stderr_sf = Data.stderr_sf; + +Data = load('C:/Users/Karthik/Documents/GitRepositories/Calculations/Data-Analyzer/RampUpSL.mat', 'unique_scan_parameter_values', 'mean_sf', 'stderr_sf'); + +up_scan_parameter_values = Data.unique_scan_parameter_values; +up_mean_sf = Data.mean_sf; +up_stderr_sf = Data.stderr_sf; + +figure(1); +set(gcf,'Position',[100 100 950 750]) +errorbar(down_scan_parameter_values, down_mean_sf, down_stderr_sf, 'o--', ... + 'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5, 'DisplayName' , 'BEC to Droplets'); +hold on +errorbar(up_scan_parameter_values, up_mean_sf, up_stderr_sf, 'o--', ... + 'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5, 'DisplayName', 'Droplets to BEC'); +set(gca, 'FontSize', 14); % For tick labels only +hXLabel = xlabel('B_z (G)', 'Interpreter', 'tex'); +hYLabel = ylabel('Spectral Weight', 'Interpreter', 'tex'); +hTitle = title('\alpha = 0', 'Interpreter', 'tex'); +legend +set([hXLabel, hYLabel], 'FontName', font) +set([hXLabel, hYLabel], 'FontSize', 14) +set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title +grid on \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally_compare_numerical_methods.m b/Dipolar-Gas-Simulator/+Scripts/run_locally_compare_numerical_methods.m new file mode 100644 index 0000000..682cec1 --- /dev/null +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally_compare_numerical_methods.m @@ -0,0 +1,127 @@ +% Compare numerical methods + +%% - Polyak's heavy-ball GD + +OptionsStruct = struct; + +OptionsStruct.NumberOfAtoms = 40000; +OptionsStruct.DipolarPolarAngle = deg2rad(0); +OptionsStruct.DipolarAzimuthAngle = 0; +OptionsStruct.ScatteringLength = 95; + +OptionsStruct.TrapFrequencies = [30, 60, 90]; +OptionsStruct.TrapPotentialType = 'Harmonic'; + +OptionsStruct.NumberOfGridPoints = [128, 64, 64]; +OptionsStruct.Dimensions = [30, 20, 20]; +OptionsStruct.UseApproximationForLHY = true; +OptionsStruct.IncludeDDICutOff = true; +OptionsStruct.CutoffType = 'Cylindrical'; +OptionsStruct.SimulationMode = 'EnergyMinimization'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization' +OptionsStruct.GradientDescentMethod = 'HeavyBall'; % 'HeavyBall' | 'NonLinearCGD' +OptionsStruct.MaxIterationsForGD = 1000; +OptionsStruct.NoiseScaleFactor = 0.010; + +OptionsStruct.PlotLive = true; +OptionsStruct.JobNumber = 0; +OptionsStruct.RunOnGPU = false; +OptionsStruct.SaveData = true; +OptionsStruct.SaveDirectory = './Results/Data_3D/GradientDescent'; +options = Helper.convertstruct2cell(OptionsStruct); + +sim = Simulator.DipolarGas(options{:}); +pot = Simulator.Potentials(options{:}); +sim.Potential = pot.trap(); + +%-% Run Simulation %-% +NumberOfOutputs = 5; +[Params, Transf, psi, V, VDk, stats] = Helper.runWithProfiling(@() sim.run(), NumberOfOutputs, OptionsStruct.SaveDirectory); +fprintf('Runtime: %.3f seconds\n', stats.runtime); +fprintf('Memory used: %.2f MB\n', stats.workspaceMemoryMB); + +clear all + +%% - Non-linear CGD + +OptionsStruct = struct; + +OptionsStruct.NumberOfAtoms = 40000; +OptionsStruct.DipolarPolarAngle = deg2rad(0); +OptionsStruct.DipolarAzimuthAngle = 0; +OptionsStruct.ScatteringLength = 95; + +OptionsStruct.TrapFrequencies = [30, 60, 90]; +OptionsStruct.TrapPotentialType = 'Harmonic'; + +OptionsStruct.NumberOfGridPoints = [128, 64, 64]; +OptionsStruct.Dimensions = [30, 20, 20]; +OptionsStruct.UseApproximationForLHY = true; +OptionsStruct.IncludeDDICutOff = true; +OptionsStruct.CutoffType = 'Cylindrical'; +OptionsStruct.SimulationMode = 'EnergyMinimization'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization' +OptionsStruct.GradientDescentMethod = 'NonLinearCGD'; % 'HeavyBall' | 'NonLinearCGD' +OptionsStruct.MaxIterationsForGD = 1000; +OptionsStruct.NoiseScaleFactor = 0.010; + +OptionsStruct.PlotLive = true; +OptionsStruct.JobNumber = 1; +OptionsStruct.RunOnGPU = false; +OptionsStruct.SaveData = true; +OptionsStruct.SaveDirectory = './Results/Data_3D/GradientDescent'; +options = Helper.convertstruct2cell(OptionsStruct); + +sim = Simulator.DipolarGas(options{:}); +pot = Simulator.Potentials(options{:}); +sim.Potential = pot.trap(); + +%-% Run Simulation %-% +NumberOfOutputs = 5; +[Params, Transf, psi, V, VDk, stats] = Helper.runWithProfiling(@() sim.run(), NumberOfOutputs, OptionsStruct.SaveDirectory); +fprintf('Runtime: %.3f seconds\n', stats.runtime); +fprintf('Memory used: %.2f MB\n', stats.workspaceMemoryMB); + +clear all + +%% - Imaginary time propagation + +OptionsStruct = struct; + +OptionsStruct.NumberOfAtoms = 40000; +OptionsStruct.DipolarPolarAngle = deg2rad(0); +OptionsStruct.DipolarAzimuthAngle = 0; +OptionsStruct.ScatteringLength = 95; + +OptionsStruct.TrapFrequencies = [30, 60, 90]; +OptionsStruct.TrapPotentialType = 'Harmonic'; + +OptionsStruct.NumberOfGridPoints = [128, 64, 64]; +OptionsStruct.Dimensions = [30, 20, 20]; +OptionsStruct.UseApproximationForLHY = true; +OptionsStruct.IncludeDDICutOff = true; +OptionsStruct.CutoffType = 'Cylindrical'; +OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization' +OptionsStruct.TimeStepSize = 1E-3; % in s +OptionsStruct.MinimumTimeStepSize = 1E-6; % in s +OptionsStruct.TimeCutOff = 1E5; % in s +OptionsStruct.EnergyTolerance = 5E-10; +OptionsStruct.ResidualTolerance = 1E-08; +OptionsStruct.NoiseScaleFactor = 0.010; + +OptionsStruct.PlotLive = true; +OptionsStruct.JobNumber = 0; +OptionsStruct.RunOnGPU = false; +OptionsStruct.SaveData = true; +OptionsStruct.SaveDirectory = './Results/Data_3D/AnisotropicTrap/Tilted0'; +options = Helper.convertstruct2cell(OptionsStruct); + +sim = Simulator.DipolarGas(options{:}); +pot = Simulator.Potentials(options{:}); +sim.Potential = pot.trap(); + +%-% Run Simulation %-% +NumberOfOutputs = 5; +[Params, Transf, psi, V, VDk, stats] = Helper.runWithProfiling(@() sim.run(), NumberOfOutputs, OptionsStruct.SaveDirectory); +fprintf('Runtime: %.3f seconds\n', stats.runtime); +fprintf('Memory used: %.2f MB\n', stats.workspaceMemoryMB); + +clear all \ No newline at end of file diff --git a/Estimations/2DQHOHeavyBallGradientDescent.m b/Estimations/2DQHOHeavyBallGradientDescent.m index 1fdcbe3..89e63c7 100644 --- a/Estimations/2DQHOHeavyBallGradientDescent.m +++ b/Estimations/2DQHOHeavyBallGradientDescent.m @@ -104,3 +104,7 @@ subplot(2,3,3); imagesc(x, y, abs(psi).^2 - abs(psi_analytic).^2); title('Differ subplot(2,3,[4 5 6]); plot(1:idx, energy, 'b-', 'LineWidth', 1.5); hold on; yline(1, 'r--', 'Analytical E=1'); xlabel('Iteration'); ylabel('Energy'); title('Energy Convergence'); grid on; legend('Numerical', 'Analytical'); + +%% + + diff --git a/Imaging-Response-Function-Extractor/extractIRF.m b/Imaging-Response-Function-Extractor/extractIRF.m index 149e1d3..2a953d8 100644 --- a/Imaging-Response-Function-Extractor/extractIRF.m +++ b/Imaging-Response-Function-Extractor/extractIRF.m @@ -11,9 +11,7 @@ folderPath = strcat(folderPath, run); cam = 5; angle = 0; -% center = [1137, 2023]; center = [1141, 2049]; -% center = [1166, 2055]; span = [255, 255]; fraction = [0.1, 0.1]; @@ -29,7 +27,7 @@ d = 1.22 * (lambda / NA); AbbeLimit = lambda / (2 * NA); % Maximum resolvable spatial frequency for the coherent case -k_cutoff = (NA/lambda) * 1e-6; % (in units of 1/µm) +k_cutoff = 2 * pi * (NA/lambda) * 1e-6; % (in units of 1/µm) removeFringes = false; @@ -38,8 +36,8 @@ removeFringes = false; filePattern = fullfile(folderPath, '*.h5'); files = dir(filePattern); -refimages = zeros(span(1) + 1, span(2) + 1, length(files)); -absimages = zeros(span(1) + 1, span(2) + 1, length(files)); +refimages = zeros(span(1) + 1, span(2) + 1, length(files)); +absimages = zeros(span(1) + 1, span(2) + 1, length(files)); for k = 1 : length(files) baseFileName = files(k).name; @@ -103,14 +101,14 @@ ky = 2*pi*vy; % Wavenumber axis in Y % Create Circular Mask n = 2^8; % size of mask -mask = zeros(n); +fftmask = zeros(n); I = 1:n; x = I-n/2; % mask x-coordinates y = n/2-I; % mask y-coordinates [X,Y] = meshgrid(x,y); % create 2-D mask grid R = 32; % aperture radius -A = (X.^2 + Y.^2 <= R^2); % circular aperture of radius R -mask(A) = 1; % set mask elements inside aperture to 1 +A = (X.^2 + Y.^2 <= R^2); % circular aperture of radius R +fftmask(A) = 1; % set mask elements inside aperture to 1 % Calculate Power Spectrum and plot figure(1) @@ -122,7 +120,7 @@ t = tiledlayout(2, 3, 'TileSpacing', 'compact', 'Padding', 'compact'); for k = 1 : length(od_imgs) mean_subtracted_od_imgs{k} = od_imgs{k} - mean_od_img; - masked_img = mean_subtracted_od_imgs{k} .* mask; + masked_img = mean_subtracted_od_imgs{k} .* fftmask; density_fft{k} = (1/numel(masked_img)) * abs(fftshift(fft2(masked_img))); density_noise_spectrum{k} = density_fft{k}.^2; @@ -161,7 +159,7 @@ for k = 1 : length(od_imgs) % Tile 4: Masked Noise nexttile(4); - imagesc(xvals, yvals, mean_subtracted_od_imgs{k} .* mask) + imagesc(xvals, yvals, mean_subtracted_od_imgs{k} .* fftmask) xlabel('\mum', 'Interpreter', 'tex', 'FontSize', 16) ylabel('\mum', 'Interpreter', 'tex', 'FontSize', 16) axis equal tight; @@ -216,27 +214,27 @@ colormap(flip(jet)); title('Average Density Noise Spectrum', 'FontSize', 16); grid on; -centers = ginput; -radius = 3; -% Plot where clicked. -hVC = viscircles(centers, radius, 'Color', 'r', 'LineWidth', 2); -xc = centers(:,1); -yc = centers(:,2); +centers = ginput; +radius = 3; +% Plot where clicked +hVC = viscircles(centers, radius, 'Color', 'r', 'LineWidth', 2); +xc = centers(:,1); +yc = centers(:,2); [yDim, xDim] = size(averagePowerSpectrum); -[xx,yy] = meshgrid(1:yDim,1:xDim); -mask = false(xDim,yDim); +[xx,yy] = meshgrid(1:yDim,1:xDim); +MaskToRemoveUnwantedPeaks = false(xDim,yDim); for ii = 1:size(centers, 1) - mask = mask | hypot(xx - xc(ii), yy - yc(ii)) <= radius; + MaskToRemoveUnwantedPeaks = MaskToRemoveUnwantedPeaks | hypot(xx - xc(ii), yy - yc(ii)) <= radius; end -mask = not(mask); +MaskToRemoveUnwantedPeaks = not(MaskToRemoveUnwantedPeaks); % Ask user if the circle is acceptable. message = sprintf('Is this acceptable?'); button = questdlg(message, message, 'Accept', 'Reject and Quit', 'Accept'); if contains(button, 'Accept','IgnoreCase',true) - image = mask.*averagePowerSpectrum; + image = MaskToRemoveUnwantedPeaks.*averagePowerSpectrum; image(image==0) = NaN; - imagesc(kx*1E-6, ky*1E-6, mask.*abs(10*log10(averagePowerSpectrum))) + imagesc(kx*1E-6, ky*1E-6, MaskToRemoveUnwantedPeaks.*abs(10*log10(averagePowerSpectrum))) hold on xlabel('k_x (\mum^{-1})', 'Interpreter', 'tex', 'FontSize', 16) ylabel('k_y (\mum^{-1})', 'Interpreter', 'tex', 'FontSize', 16) @@ -386,101 +384,6 @@ xlim([1, xSize]); ylim([1, ySize]); zlim([min(10 * log10(averagePowerSpectrum(:))), max(10 * log10(averagePowerSpectrum(:)))]); % Optional for Z-axis limits - -%% Decompose in Zernike Polynomial basis - -N = size(averagePowerSpectrum, 1); -[X, Y] = meshgrid(linspace(-1, 1, N)); - -max_n = 6; % Adjust based on your needs -basis = []; -orders = []; - -for n = 0:max_n - for m = -n:2:n - % Generate Zernike polynomial for (n, m) - Z = zernike_polynomial(n, m, X, Y); - % Flatten and store valid points - basis = [basis, Z(mask)]; - orders = [orders; [n, m]]; - end -end - -data = 10 * log10(averagePowerSpectrum); -valid_data = data(mask); - -% Solve Ax = b (A = basis matrix, b = data) -coeffs = basis \ valid_data(:); - -% Reconstruct the surface using the coefficients -reconstructed = basis * coeffs; -reconstructed_surface = zeros(size(X)); -reconstructed_surface(mask) = reconstructed; - -figure(5) -clf -set(gcf,'Position',[100, 100, 1500, 700]) - -% Create tiled layout with 2 rows and 3 columns -t = tiledlayout(1, 3, 'TileSpacing', 'compact', 'Padding', 'compact'); - -nexttile(1); -imagesc(data); title('Imaging Response Function', 'FontSize', 16); -axis square; -colorbar -colormap(jet); -grid on; - -nexttile(2); -imagesc(reconstructed_surface); title('Reconstructed with Zernike', 'FontSize', 16); -axis square; -colorbar -colormap(jet); -grid on; - -nexttile(3); -imagesc(data - reconstructed_surface); title('Residuals', 'FontSize', 16); -axis square; -colorbar -colormap(jet); -grid on; - -disp('Zernike Coefficients:'); -disp('---------------------'); -for i = 1:length(coeffs) - fprintf('Order (n=%d, m=%d): Coefficient = %.4f\n', orders(i,1), orders(i,2), coeffs(i)); -end - -% Plot Zernike Coeffecients - -% Find the index of the (n=0, m=0) term -idx_remove = find(orders(:,1) == 0 & orders(:,2) == 0); - -% Remove the Z_0^0 term from coefficients and orders -coeffs_filtered = coeffs; -coeffs_filtered(idx_remove) = []; -orders_filtered = orders; -orders_filtered(idx_remove, :) = []; - -% Generate labels for filtered modes (n, m) -labels_filtered = cell(length(coeffs_filtered), 1); -for i = 1:length(coeffs_filtered) - labels_filtered{i} = sprintf('(%d, %d)', orders_filtered(i,1), orders_filtered(i,2)); -end - -figure(6) -clf -set(gcf,'Position',[100, 100, 1500, 700]) -bar(coeffs_filtered, 'FaceColor', [0.2, 0.6, 0.8]); % Customize bar color -ylim([-2.0, 2.0]) -title('Zernike Coefficients', 'FontSize', 16); -xlabel('Zernike Mode (n, m)', 'FontSize', 16); -ylabel('Coefficient Value', 'FontSize', 16); -xticks(1:length(coeffs_filtered)); % Set x-ticks for all coefficients -xticklabels(labels_filtered); % Assign (n, m) labels -xtickangle(45); % Rotate labels for readability -grid on; - %% Helper Functions function ret = getBkgOffsetFromCorners(img, x_fraction, y_fraction)