karthik/Calculations
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

New additions, minor modifications to existing files

Browse Source
This commit is contained in:
Karthik 2025-04-30 12:45:20 +02:00
parent f9bc6d032c
commit 305335f519
5 changed files with 227 additions and 117 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -9,5 +9,6 @@ Time-Series-Analyzer/Time-Series-Data
*.gif
*.mp4
*.bat
*.json
.ipynb_checkpoints/
.vscode/

75
Data-Analyzer/compareSpectralWeights.m Normal file
View File

@ -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

127
Dipolar-Gas-Simulator/+Scripts/run_locally_compare_numerical_methods.m Normal file
View File

@ -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

4
Estimations/2DQHOHeavyBallGradientDescent.m
View File

@ -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');
%%

137
Imaging-Response-Function-Extractor/extractIRF.m
View File

@ -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)