Merge branch 'master' of https://git.physi.uni-heidelberg.de/karthik/Calculations
This commit is contained in:
Karthik
commit
8ff8fbae58
554
Data-Analyzer/analyzePhaseTransitionSimulation.m
Normal file
554
Data-Analyzer/analyzePhaseTransitionSimulation.m
Normal file
@ -0,0 +1,554 @@
|
||||
%% Extract Images
|
||||
|
||||
baseDir = 'D:/Results - Numerics/Data_Full3D/PhaseTransition/DTS/';
|
||||
JobNumber = 0;
|
||||
runFolder = sprintf('Run_%03d', JobNumber);
|
||||
movieFileName = 'DropletsToStripes.mp4'; % Output file name
|
||||
datafileName = './DropletsToStripes.mat';
|
||||
reverseOrder = false; % Set this to true to reverse the theta ordering
|
||||
|
||||
TitleString = 'Change across transition: Droplets To Stripes';
|
||||
|
||||
%%
|
||||
baseDir = 'D:/Results - Numerics/Data_Full3D/PhaseTransition/STD/';
|
||||
JobNumber = 0;
|
||||
runFolder = sprintf('Run_%03d', JobNumber);
|
||||
movieFileName = 'StripesToDroplets.mp4'; % Output file name
|
||||
datafileName = './StripesToDroplets.mat';
|
||||
reverseOrder = true; % Set this to true to reverse the theta ordering
|
||||
|
||||
|
||||
TitleString = 'Change across transition: Stripes To Droplets';
|
||||
|
||||
%%
|
||||
folderList = dir(baseDir);
|
||||
isValid = [folderList.isdir] & ~ismember({folderList.name}, {'.', '..'});
|
||||
folderNames = {folderList(isValid).name};
|
||||
nimgs = numel(folderNames);
|
||||
|
||||
% Extract theta values from folder names
|
||||
PolarAngleVals = zeros(1, nimgs);
|
||||
for k = 1:nimgs
|
||||
tokens = regexp(folderNames{k}, 'theta_(\d{3})', 'tokens');
|
||||
if isempty(tokens)
|
||||
warning('No theta found in folder name: %s', folderNames{k});
|
||||
PolarAngleVals(k) = NaN;
|
||||
else
|
||||
PolarAngleVals(k) = str2double(tokens{1}{1});
|
||||
end
|
||||
end
|
||||
|
||||
% Choose sort direction
|
||||
sortDirection = 'ascend';
|
||||
if reverseOrder
|
||||
sortDirection = 'descend';
|
||||
end
|
||||
|
||||
% Sort folderNames based on polar angle
|
||||
[~, sortIdx] = sort(PolarAngleVals, sortDirection);
|
||||
folderNames = folderNames(sortIdx);
|
||||
PolarAngleVals = PolarAngleVals(sortIdx); % Optional: if you still want sorted list
|
||||
|
||||
imgs = cell(1, nimgs);
|
||||
alphas = zeros(1, nimgs);
|
||||
|
||||
for k = 1:nimgs
|
||||
folderName = folderNames{k};
|
||||
SaveDirectory = fullfile(baseDir, folderName, runFolder);
|
||||
|
||||
% Extract alpha (theta) again from folder name
|
||||
tokens = regexp(folderName, 'theta_(\d{3})', 'tokens');
|
||||
alpha_val = str2double(tokens{1}{1});
|
||||
alphas(k) = alpha_val;
|
||||
|
||||
matPath = fullfile(SaveDirectory, 'psi_gs.mat');
|
||||
if ~isfile(matPath)
|
||||
warning('Missing psi_gs.mat in %s', SaveDirectory);
|
||||
continue;
|
||||
end
|
||||
|
||||
try
|
||||
Data = load(matPath, 'psi', 'Params', 'Transf', 'Observ');
|
||||
catch ME
|
||||
warning('Failed to load %s: %s', matPath, ME.message);
|
||||
continue;
|
||||
end
|
||||
|
||||
Params = Data.Params;
|
||||
Transf = Data.Transf;
|
||||
Observ = Data.Observ;
|
||||
|
||||
psi = Data.psi;
|
||||
if isgpuarray(psi)
|
||||
psi = gather(psi);
|
||||
end
|
||||
if isgpuarray(Observ.residual)
|
||||
Observ.residual = gather(Observ.residual);
|
||||
end
|
||||
|
||||
% Axes and projection
|
||||
x = Transf.x * Params.l0 * 1e6;
|
||||
y = Transf.y * Params.l0 * 1e6;
|
||||
z = Transf.z * Params.l0 * 1e6;
|
||||
dx = x(2)-x(1); dy = y(2)-y(1); dz = z(2)-z(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 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
|
||||
|
||||
n = abs(psi).^2;
|
||||
nxy = squeeze(trapz(n * dz, 3));
|
||||
|
||||
imgs{k} = nxy;
|
||||
end
|
||||
|
||||
%% Analyze Images
|
||||
|
||||
|
||||
makeMovie = true; % Set to false to disable movie creation
|
||||
font = 'Bahnschrift';
|
||||
|
||||
skipPreprocessing = true;
|
||||
skipMasking = true;
|
||||
skipIntensityThresholding = true;
|
||||
skipBinarization = true;
|
||||
|
||||
% Run Fourier analysis over images
|
||||
|
||||
fft_imgs = cell(1, nimgs);
|
||||
spectral_contrast = zeros(1, nimgs);
|
||||
spectral_weight = zeros(1, nimgs);
|
||||
g2_all = cell(1, nimgs);
|
||||
theta_values_all = cell(1, nimgs);
|
||||
|
||||
N_bins = 180;
|
||||
Threshold = 25;
|
||||
Sigma = 2;
|
||||
|
||||
if makeMovie
|
||||
% Create VideoWriter object for movie
|
||||
videoFile = VideoWriter(movieFileName, 'MPEG-4');
|
||||
videoFile.Quality = 100; % Set quality to maximum (0–100)
|
||||
videoFile.FrameRate = 2; % Set the frame rate (frames per second)
|
||||
open(videoFile); % Open the video file to write
|
||||
end
|
||||
|
||||
% Display the cropped image
|
||||
for k = 1:nimgs
|
||||
IMG = imgs{k};
|
||||
[IMGFFT, IMGPR] = computeFourierTransform(IMG, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization);
|
||||
[theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(IMGFFT, 10, 35, N_bins, Threshold, Sigma);
|
||||
|
||||
g2 = zeros(1, N_bins); % Preallocate
|
||||
|
||||
for dtheta = 0:N_bins-1
|
||||
profile = S_theta;
|
||||
profile_shifted = circshift(profile, -dtheta, 2);
|
||||
|
||||
num = mean(profile .* profile_shifted);
|
||||
denom = mean(profile)^2;
|
||||
|
||||
g2(dtheta+1) = num / denom - 1;
|
||||
end
|
||||
|
||||
g2_all{k} = g2;
|
||||
theta_values_all{k} = theta_vals;
|
||||
|
||||
figure(1);
|
||||
clf
|
||||
set(gcf,'Position',[500 100 1000 800])
|
||||
t = tiledlayout(2, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 1x4 grid
|
||||
|
||||
y_min = min(y);
|
||||
y_max = max(y);
|
||||
x_min = min(x);
|
||||
x_max = max(x);
|
||||
|
||||
% Display the cropped OD image
|
||||
ax1 = nexttile;
|
||||
imagesc(x, y, IMG')
|
||||
% 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-right corner (normalized axis coordinates)
|
||||
hText = text(1 - x_offset, 1 - y_offset, ['Angle = ', num2str(alphas(k), '%.1f')], ...
|
||||
'Color', 'white', 'FontWeight', 'bold', 'Interpreter', 'tex', 'FontSize', 20, 'Units', 'normalized', 'HorizontalAlignment', 'right', 'VerticalAlignment', 'top');
|
||||
axis square;
|
||||
hcb = colorbar;
|
||||
colormap(ax1, 'jet');
|
||||
set(gca, 'FontSize', 14); % For tick labels only
|
||||
hL = ylabel(hcb, 'Optical Density');
|
||||
set(hL,'Rotation',-90);
|
||||
set(gca,'YDir','normal')
|
||||
% set(gca, 'YTick', linspace(y_min, y_max, 5)); % Define y ticks
|
||||
% set(gca, 'YTickLabel', flip(linspace(y_min, y_max, 5))); % Flip only the labels
|
||||
hXLabel = xlabel('x (pixels)', 'Interpreter', 'tex');
|
||||
hYLabel = ylabel('y (pixels)', 'Interpreter', 'tex');
|
||||
hTitle = title('Density', 'Interpreter', 'tex');
|
||||
set([hXLabel, hYLabel, hL, hText], 'FontName', font)
|
||||
set([hXLabel, hYLabel, hL], 'FontSize', 14)
|
||||
set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title
|
||||
|
||||
% Plot the power spectrum
|
||||
ax2 = nexttile;
|
||||
[rows, cols] = size(IMGFFT);
|
||||
zoom_size = 50; % Zoomed-in region around center
|
||||
mid_x = floor(cols/2);
|
||||
mid_y = floor(rows/2);
|
||||
fft_imgs{k} = IMGFFT(mid_y-zoom_size:mid_y+zoom_size, mid_x-zoom_size:mid_x+zoom_size);
|
||||
imagesc(log(1 + abs(fft_imgs{k}).^2));
|
||||
% 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
|
||||
axis square;
|
||||
hcb = colorbar;
|
||||
colormap(ax2, 'jet');
|
||||
set(gca, 'FontSize', 14); % For tick labels only
|
||||
set(gca,'YDir','normal')
|
||||
hXLabel = xlabel('k_x', 'Interpreter', 'tex');
|
||||
hYLabel = ylabel('k_y', 'Interpreter', 'tex');
|
||||
hTitle = title('Power Spectrum - S(k_x,k_y)', 'Interpreter', 'tex');
|
||||
set([hXLabel, hYLabel, hText], 'FontName', font)
|
||||
set([hXLabel, hYLabel], 'FontSize', 14)
|
||||
set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title
|
||||
|
||||
% Plot the angular distribution
|
||||
nexttile
|
||||
spectral_contrast(k) = computeSpectralContrast(fft_imgs{k}, 10, 25, Threshold);
|
||||
[theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(fft_imgs{k}, 10, 25, N_bins, Threshold, Sigma);
|
||||
spectral_weight(k) = trapz(theta_vals, S_theta);
|
||||
plot(theta_vals/pi, S_theta,'Linewidth',2);
|
||||
axis square;
|
||||
set(gca, 'FontSize', 14); % For tick labels only
|
||||
hXLabel = xlabel('\theta/\pi [rad]', 'Interpreter', 'tex');
|
||||
hYLabel = ylabel('Normalized magnitude (a.u.)', 'Interpreter', 'tex');
|
||||
hTitle = title('Angular Spectral Distribution - S(\theta)', 'Interpreter', 'tex');
|
||||
set([hXLabel, hYLabel, hText], 'FontName', font)
|
||||
set([hXLabel, hYLabel], 'FontSize', 14)
|
||||
set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title
|
||||
grid on
|
||||
|
||||
nexttile
|
||||
plot(theta_vals/pi, g2, 'o-', 'LineWidth', 1.2, 'MarkerSize', 5);
|
||||
set(gca, 'FontSize', 14);
|
||||
ylim([-1.5 3.0]); % Set y-axis limits here
|
||||
hXLabel = xlabel('$\delta\theta / \pi$', 'Interpreter', 'latex');
|
||||
hYLabel = ylabel('$g^{(2)}(\delta\theta)$', 'Interpreter', 'latex');
|
||||
hTitle = title('Autocorrelation', 'Interpreter', 'tex');
|
||||
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;
|
||||
|
||||
if makeMovie
|
||||
frame = getframe(gcf);
|
||||
writeVideo(videoFile, frame);
|
||||
else
|
||||
pause(0.5); % Only pause when not recording
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
if makeMovie
|
||||
close(videoFile);
|
||||
disp(['Movie saved to ', movieFileName]);
|
||||
end
|
||||
|
||||
%% Track across the transition
|
||||
|
||||
figure(2);
|
||||
set(gcf,'Position',[100 100 950 750])
|
||||
plot(alphas, spectral_contrast, 'o--', ...
|
||||
'LineWidth', 1.5, 'MarkerSize', 6);
|
||||
set(gca, 'FontSize', 14); % For tick labels only
|
||||
hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex');
|
||||
% hXLabel = xlabel('B_z (G)', 'Interpreter', 'tex');
|
||||
hYLabel = ylabel('Spectral Contrast', 'Interpreter', 'tex');
|
||||
hTitle = title(TitleString, 'Interpreter', 'tex');
|
||||
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
|
||||
|
||||
figure(3);
|
||||
set(gcf,'Position',[100 100 950 750])
|
||||
plot(alphas, spectral_weight, 'o--', ...
|
||||
'LineWidth', 1.5, 'MarkerSize', 6);
|
||||
set(gca, 'FontSize', 14); % For tick labels only
|
||||
hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex');
|
||||
% hXLabel = xlabel('B_z (G)', 'Interpreter', 'tex');
|
||||
hYLabel = ylabel('Spectral Weight', 'Interpreter', 'tex');
|
||||
hTitle = title(TitleString, 'Interpreter', 'tex');
|
||||
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
|
||||
|
||||
save(datafileName, 'alphas', 'spectral_contrast', 'spectral_weight');
|
||||
|
||||
figure(4);
|
||||
clf;
|
||||
set(gcf,'Position',[100 100 950 750])
|
||||
hold on;
|
||||
|
||||
% Reconstruct theta axis from any one of the stored values
|
||||
theta_vals = theta_values_all{1}; % assuming it's in radians
|
||||
|
||||
legend_entries = cell(nimgs, 1);
|
||||
|
||||
% Generate a colormap with enough unique colors
|
||||
cmap = sky(nimgs); % You can also try 'jet', 'turbo', 'hot', etc.
|
||||
|
||||
for i = 1:nimgs
|
||||
plot(theta_vals/pi, g2_all{i}, ...
|
||||
'o-', 'Color', cmap(i,:), 'LineWidth', 1.2, ...
|
||||
'MarkerSize', 5);
|
||||
legend_entries{i} = sprintf('$\\alpha = %g^\\circ$', alphas(i));
|
||||
end
|
||||
|
||||
ylim([-1.5 3.0]); % Set y-axis limits here
|
||||
set(gca, 'FontSize', 14);
|
||||
hXLabel = xlabel('$\delta\theta / \pi$', 'Interpreter', 'latex');
|
||||
hYLabel = ylabel('$g^{(2)}(\delta\theta)$', 'Interpreter', 'latex');
|
||||
hTitle = title(TitleString, 'Interpreter', 'tex');
|
||||
legend(legend_entries, 'Interpreter', 'latex', 'Location', 'bestoutside');
|
||||
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 across the transition
|
||||
|
||||
set(0,'defaulttextInterpreter','latex')
|
||||
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
|
||||
|
||||
format long
|
||||
|
||||
font = 'Bahnschrift';
|
||||
|
||||
% Load data
|
||||
Data = load('./DropletsToStripes.mat', 'alphas', 'spectral_contrast', 'spectral_weight');
|
||||
dts_alphas = Data.alphas;
|
||||
dts_sc = Data.spectral_contrast;
|
||||
dts_sw = Data.spectral_weight;
|
||||
|
||||
Data = load('./StripesToDroplets.mat', 'alphas', 'spectral_contrast', 'spectral_weight');
|
||||
std_alphas = Data.alphas;
|
||||
std_sc = Data.spectral_contrast;
|
||||
std_sw = Data.spectral_weight;
|
||||
|
||||
% Normalize dts data
|
||||
dts_min = min(dts_sw);
|
||||
dts_max = max(dts_sw);
|
||||
dts_range = dts_max - dts_min;
|
||||
dts_sf_norm = (dts_sw - dts_min) / dts_range;
|
||||
|
||||
% Normalize std data
|
||||
std_min = min(std_sw);
|
||||
std_max = max(std_sw);
|
||||
std_range = std_max - std_min;
|
||||
std_sf_norm = (std_sw - std_min) / std_range;
|
||||
|
||||
figure(5);
|
||||
set(gcf,'Position',[100 100 950 750])
|
||||
plot(dts_alphas, dts_sc, 'o--', 'LineWidth', 1.5, 'MarkerSize', 6, 'DisplayName' , 'Droplets to Stripes');
|
||||
hold on
|
||||
plot(std_alphas, std_sc, 'o--', 'LineWidth', 1.5, 'MarkerSize', 6, 'DisplayName' , 'Stripes to Droplets');
|
||||
set(gca, 'FontSize', 14); % For tick labels only
|
||||
hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex');
|
||||
hYLabel = ylabel('Spectral Contrast', 'Interpreter', 'tex');
|
||||
hTitle = title('Change across transition', '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
|
||||
|
||||
figure(6);
|
||||
set(gcf,'Position',[100 100 950 750])
|
||||
plot(dts_alphas, dts_sw, 'o--', 'LineWidth', 1.5, 'MarkerSize', 6, 'DisplayName' , 'Droplets to Stripes');
|
||||
hold on
|
||||
plot(std_alphas, std_sw, 'o--', 'LineWidth', 1.5, 'MarkerSize', 6, '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('Change across transition', '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
|
||||
|
||||
%%
|
||||
function [IMGFFT, IMGPR] = computeFourierTransform(I, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization)
|
||||
% computeFourierSpectrum - Computes the 2D Fourier power spectrum
|
||||
% of binarized and enhanced lattice image features, with optional central mask.
|
||||
%
|
||||
% Inputs:
|
||||
% I - Grayscale or RGB image matrix
|
||||
%
|
||||
% Output:
|
||||
% F_mag - 2D Fourier power spectrum (shifted)
|
||||
|
||||
if ~skipPreprocessing
|
||||
% Preprocessing: Denoise
|
||||
filtered = imgaussfilt(I, 10);
|
||||
IMGPR = I - filtered; % adjust sigma as needed
|
||||
else
|
||||
IMGPR = I;
|
||||
end
|
||||
|
||||
if ~skipMasking
|
||||
[rows, cols] = size(IMGPR);
|
||||
[X, Y] = meshgrid(1:cols, 1:rows);
|
||||
% Elliptical mask parameters
|
||||
cx = cols / 2;
|
||||
cy = rows / 2;
|
||||
|
||||
% Shifted coordinates
|
||||
x = X - cx;
|
||||
y = Y - cy;
|
||||
|
||||
% Ellipse semi-axes
|
||||
rx = 0.4 * cols;
|
||||
ry = 0.2 * rows;
|
||||
|
||||
% Rotation angle in degrees -> radians
|
||||
theta_deg = 30; % Adjust as needed
|
||||
theta = deg2rad(theta_deg);
|
||||
|
||||
% Rotated ellipse equation
|
||||
cos_t = cos(theta);
|
||||
sin_t = sin(theta);
|
||||
|
||||
x_rot = (x * cos_t + y * sin_t);
|
||||
y_rot = (-x * sin_t + y * cos_t);
|
||||
|
||||
ellipseMask = (x_rot.^2) / rx^2 + (y_rot.^2) / ry^2 <= 1;
|
||||
|
||||
% Apply cutout mask
|
||||
IMGPR = IMGPR .* ellipseMask;
|
||||
end
|
||||
|
||||
if ~skipIntensityThresholding
|
||||
% Apply global intensity threshold mask
|
||||
intensity_thresh = 0.20;
|
||||
intensity_mask = IMGPR > intensity_thresh;
|
||||
IMGPR = IMGPR .* intensity_mask;
|
||||
end
|
||||
|
||||
if ~skipBinarization
|
||||
% Adaptive binarization and cleanup
|
||||
IMGPR = imbinarize(IMGPR, 'adaptive', 'Sensitivity', 0.0);
|
||||
IMGPR = imdilate(IMGPR, strel('disk', 2));
|
||||
IMGPR = imerode(IMGPR, strel('disk', 1));
|
||||
IMGPR = imfill(IMGPR, 'holes');
|
||||
F = fft2(double(IMGPR)); % Compute 2D Fourier Transform
|
||||
IMGFFT = abs(fftshift(F))'; % Shift zero frequency to center
|
||||
else
|
||||
F = fft2(double(IMGPR)); % Compute 2D Fourier Transform
|
||||
IMGFFT = abs(fftshift(F))'; % Shift zero frequency to center
|
||||
end
|
||||
end
|
||||
|
||||
function [theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(IMGFFT, r_min, r_max, num_bins, threshold, sigma)
|
||||
% Apply threshold to isolate strong peaks
|
||||
IMGFFT(IMGFFT < threshold) = 0;
|
||||
|
||||
% Prepare polar coordinates
|
||||
[ny, nx] = size(IMGFFT);
|
||||
[X, Y] = meshgrid(1:nx, 1:ny);
|
||||
cx = ceil(nx/2);
|
||||
cy = ceil(ny/2);
|
||||
R = sqrt((X - cx).^2 + (Y - cy).^2);
|
||||
Theta = atan2(Y - cy, X - cx); % range [-pi, pi]
|
||||
|
||||
% Choose radial band
|
||||
radial_mask = (R >= r_min) & (R <= r_max);
|
||||
|
||||
% Initialize the angular structure factor array
|
||||
S_theta = zeros(1, num_bins); % Pre-allocate for 180 angle bins
|
||||
% Define the angle values for the x-axis
|
||||
theta_vals = linspace(0, pi, num_bins);
|
||||
|
||||
% Loop through each angle bin
|
||||
for i = 1:num_bins
|
||||
angle_start = (i-1) * pi / num_bins;
|
||||
angle_end = i * pi / num_bins;
|
||||
|
||||
% Define a mask for the given angle range
|
||||
angle_mask = (Theta >= angle_start & Theta < angle_end);
|
||||
|
||||
bin_mask = radial_mask & angle_mask;
|
||||
|
||||
% Extract the Fourier components for the given angle
|
||||
fft_angle = IMGFFT .* bin_mask;
|
||||
|
||||
% Integrate the Fourier components over the radius at the angle
|
||||
S_theta(i) = sum(sum(abs(fft_angle).^2)); % sum of squared magnitudes
|
||||
end
|
||||
|
||||
% Create a 1D Gaussian kernel
|
||||
half_width = ceil(3 * sigma);
|
||||
x = -half_width:half_width;
|
||||
gauss_kernel = exp(-x.^2 / (2 * sigma^2));
|
||||
gauss_kernel = gauss_kernel / sum(gauss_kernel); % normalize
|
||||
|
||||
% Apply convolution (circular padding to preserve periodicity)
|
||||
S_theta = conv([S_theta(end-half_width+1:end), S_theta, S_theta(1:half_width)], gauss_kernel, 'same');
|
||||
S_theta = S_theta(half_width+1:end-half_width); % crop back to original size
|
||||
|
||||
% Normalize to 1
|
||||
S_theta = S_theta / max(S_theta);
|
||||
end
|
||||
|
||||
function contrast = computeSpectralContrast(IMGFFT, r_min, r_max, threshold)
|
||||
% Apply threshold to isolate strong peaks
|
||||
IMGFFT(IMGFFT < threshold) = 0;
|
||||
|
||||
% Prepare polar coordinates
|
||||
[ny, nx] = size(IMGFFT);
|
||||
[X, Y] = meshgrid(1:nx, 1:ny);
|
||||
cx = ceil(nx/2);
|
||||
cy = ceil(ny/2);
|
||||
R = sqrt((X - cx).^2 + (Y - cy).^2);
|
||||
|
||||
% Ring region (annulus) mask
|
||||
ring_mask = (R >= r_min) & (R <= r_max);
|
||||
|
||||
% Squared magnitude in the ring
|
||||
ring_power = abs(IMGFFT).^2 .* ring_mask;
|
||||
|
||||
% Maximum power in the ring
|
||||
ring_max = max(ring_power(:));
|
||||
|
||||
% Power at the DC component
|
||||
dc_power = abs(IMGFFT(cy, cx))^2;
|
||||
|
||||
% Avoid division by zero
|
||||
if dc_power == 0
|
||||
contrast = Inf; % or NaN or 0, depending on how you want to handle this
|
||||
else
|
||||
contrast = ring_max / dc_power;
|
||||
end
|
||||
end
|
||||
@ -8,42 +8,63 @@ format long
|
||||
font = 'Bahnschrift';
|
||||
|
||||
% Load data
|
||||
Data = load('C:/Users/Karthik/Documents/GitRepositories/Calculations/Data-Analyzer/B2.45G/DropletsToStripes.mat', 'unique_scan_parameter_values', 'mean_sf', 'stderr_sf');
|
||||
Data = load('D:/Results - Experiment/B2.45G/DropletsToStripes.mat', 'unique_scan_parameter_values', 'mean_sc', 'stderr_sc', 'mean_sw', 'stderr_sw');
|
||||
|
||||
dts_scan_parameter_values = Data.unique_scan_parameter_values;
|
||||
dts_mean_sf = Data.mean_sf;
|
||||
dts_stderr_sf = Data.stderr_sf;
|
||||
dts_mean_sc = Data.mean_sc;
|
||||
dts_stderr_sc = Data.stderr_sc;
|
||||
dts_mean_sw = Data.mean_sw;
|
||||
dts_stderr_sw = Data.stderr_sw;
|
||||
|
||||
Data = load('C:/Users/Karthik/Documents/GitRepositories/Calculations/Data-Analyzer/B2.45G/StripesToDroplets.mat', 'unique_scan_parameter_values', 'mean_sf', 'stderr_sf');
|
||||
Data = load('D:/Results - Experiment/B2.45G/StripesToDroplets.mat', 'unique_scan_parameter_values', 'mean_sc', 'stderr_sc', 'mean_sw', 'stderr_sw');
|
||||
|
||||
std_scan_parameter_values = Data.unique_scan_parameter_values;
|
||||
std_mean_sf = Data.mean_sf;
|
||||
std_stderr_sf = Data.stderr_sf;
|
||||
std_mean_sw = Data.mean_sw;
|
||||
std_stderr_sw = Data.stderr_sw;
|
||||
std_mean_sc = Data.mean_sc;
|
||||
std_stderr_sc = Data.stderr_sc;
|
||||
|
||||
% Normalize dts data
|
||||
dts_min = min(dts_mean_sf);
|
||||
dts_max = max(dts_mean_sf);
|
||||
dts_min = min(dts_mean_sw);
|
||||
dts_max = max(dts_mean_sw);
|
||||
dts_range = dts_max - dts_min;
|
||||
dts_mean_sf_norm = (dts_mean_sf - dts_min) / dts_range;
|
||||
dts_stderr_sf_norm = dts_stderr_sf / dts_range;
|
||||
dts_mean_sw_norm = (dts_mean_sw - dts_min) / dts_range;
|
||||
dts_stderr_sw_norm = dts_stderr_sw / dts_range;
|
||||
|
||||
% Normalize std data
|
||||
std_min = min(std_mean_sf);
|
||||
std_max = max(std_mean_sf);
|
||||
std_min = min(std_mean_sw);
|
||||
std_max = max(std_mean_sw);
|
||||
std_range = std_max - std_min;
|
||||
std_mean_sf_norm = (std_mean_sf - std_min) / std_range;
|
||||
std_stderr_sf_norm = std_stderr_sf / std_range;
|
||||
std_mean_sw_norm = (std_mean_sw - std_min) / std_range;
|
||||
std_stderr_sw_norm = std_stderr_sw / std_range;
|
||||
|
||||
figure(1);
|
||||
set(gcf,'Position',[100 100 950 750])
|
||||
errorbar(dts_scan_parameter_values, dts_mean_sf_norm, dts_stderr_sf_norm, 'o--', ...
|
||||
errorbar(dts_scan_parameter_values, dts_mean_sc, dts_stderr_sc, 'o--', ...
|
||||
'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5, 'DisplayName' , 'Droplets to Stripes');
|
||||
hold on
|
||||
errorbar(std_scan_parameter_values, std_mean_sf_norm, std_stderr_sf_norm, 'o--', ...
|
||||
errorbar(std_scan_parameter_values, std_mean_sc, std_stderr_sc, '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('Normalized Spectral Weight', 'Interpreter', 'tex');
|
||||
hYLabel = ylabel('Spectral Contrast', 'Interpreter', 'tex');
|
||||
hTitle = title('B = 2.45 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
|
||||
|
||||
figure(2);
|
||||
set(gcf,'Position',[100 100 950 750])
|
||||
errorbar(dts_scan_parameter_values, dts_mean_sw, dts_stderr_sw, 'o--', ...
|
||||
'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5, 'DisplayName' , 'Droplets to Stripes');
|
||||
hold on
|
||||
errorbar(std_scan_parameter_values, std_mean_sw, std_stderr_sw, '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.45 G', 'Interpreter', 'tex');
|
||||
legend
|
||||
set([hXLabel, hYLabel], 'FontName', font)
|
||||
|
||||
179
Data-Analyzer/computeBondOrderParameters.m
Normal file
179
Data-Analyzer/computeBondOrderParameters.m
Normal file
@ -0,0 +1,179 @@
|
||||
function [psi_global, order_ratio, Npoints] = computeBondOrderParameters(I)
|
||||
% Analyze bond-orientational order in a 2D image, restricted to a circular region
|
||||
%
|
||||
% Input:
|
||||
% I - 2D grayscale image (double or uint8)
|
||||
%
|
||||
% Outputs:
|
||||
% psi_global - struct with fields psi2, psi4, psi6 (⟨|ψₙ|⟩ values)
|
||||
% order_ratio - scalar = ⟨|ψ₆|⟩ / ⟨|ψ₂|⟩
|
||||
% Npoints - number of detected points used
|
||||
|
||||
%% Step 1: Create mask
|
||||
[Ny, Nx] = size(I);
|
||||
% Parameters for elliptical mask
|
||||
a = Nx / 2.25; % Semi-major axis (x-direction)
|
||||
b = Ny / 2.25; % Semi-minor axis (y-direction)
|
||||
theta_deg = 0; % Rotation angle in degrees
|
||||
theta = deg2rad(theta_deg); % Convert to radians
|
||||
|
||||
% Meshgrid coordinates
|
||||
[Ny, Nx] = size(I);
|
||||
[X, Y] = meshgrid(1:Nx, 1:Ny);
|
||||
cx = Nx / 2;
|
||||
cy = Ny / 2;
|
||||
|
||||
% Shifted coordinates
|
||||
Xc = X - cx;
|
||||
Yc = Y - cy;
|
||||
|
||||
% Rotate coordinates
|
||||
Xr = cos(theta) * Xc + sin(theta) * Yc;
|
||||
Yr = -sin(theta) * Xc + cos(theta) * Yc;
|
||||
|
||||
% Elliptical mask equation
|
||||
mask = (Xr / a).^2 + (Yr / b).^2 <= 1;
|
||||
|
||||
% Apply mask to image
|
||||
I = double(I) .* mask;
|
||||
|
||||
%% Step 2: Detect local maxima within mask
|
||||
I_smooth = imgaussfilt(I, 2);
|
||||
peaks = imregionalmax(I_smooth);
|
||||
[y, x] = find(peaks);
|
||||
Npoints = length(x);
|
||||
|
||||
if Npoints < 6
|
||||
warning('Too few points (%d) to compute bond order meaningfully.', Npoints);
|
||||
psi_global = struct('psi2', NaN, 'psi4', NaN, 'psi6', NaN);
|
||||
order_ratio = NaN;
|
||||
return;
|
||||
end
|
||||
|
||||
%% Step 3: Try Delaunay triangulation
|
||||
use_kNN = false;
|
||||
try
|
||||
tri = delaunay(x, y);
|
||||
catch
|
||||
warning('Delaunay triangulation failed - falling back to kNN neighbor search');
|
||||
use_kNN = true;
|
||||
end
|
||||
|
||||
% 3) Find neighbors
|
||||
neighbors = cell(Npoints, 1);
|
||||
if ~use_kNN
|
||||
% Use Delaunay neighbors
|
||||
for i = 1:size(tri,1)
|
||||
for j = 1:3
|
||||
idx = tri(i,j);
|
||||
nbrs = tri(i,[1 2 3] ~= j);
|
||||
neighbors{idx} = unique([neighbors{idx}, nbrs]);
|
||||
end
|
||||
end
|
||||
else
|
||||
% Use kNN neighbors (e.g. k=6)
|
||||
k = min(6, Npoints-1);
|
||||
pts = [x, y];
|
||||
[idx, ~] = knnsearch(pts, pts, 'K', k+1); % +1 for self
|
||||
for j = 1:Npoints
|
||||
neighbors{j} = idx(j, 2:end);
|
||||
end
|
||||
end
|
||||
|
||||
% 4) Compute bond order parameters
|
||||
psi2 = zeros(Npoints,1);
|
||||
psi4 = zeros(Npoints,1);
|
||||
psi6 = zeros(Npoints,1);
|
||||
|
||||
for j = 1:Npoints
|
||||
nbrs = neighbors{j};
|
||||
nbrs(nbrs == j) = [];
|
||||
if isempty(nbrs)
|
||||
continue;
|
||||
end
|
||||
angles = atan2(y(nbrs) - y(j), x(nbrs) - x(j));
|
||||
psi2(j) = abs(mean(exp(1i * 2 * angles)));
|
||||
psi4(j) = abs(mean(exp(1i * 4 * angles)));
|
||||
psi6(j) = abs(mean(exp(1i * 6 * angles)));
|
||||
end
|
||||
|
||||
psi_global.psi2 = mean(psi2);
|
||||
psi_global.psi4 = mean(psi4);
|
||||
psi_global.psi6 = mean(psi6);
|
||||
|
||||
order_ratio = psi_global.psi6 / psi_global.psi2;
|
||||
|
||||
% --- After neighbors are found and (x,y) extracted ---
|
||||
|
||||
% 1) Angular distribution of all bonds pooled from all points
|
||||
all_angles = [];
|
||||
for j = 1:Npoints
|
||||
nbrs = neighbors{j};
|
||||
nbrs(nbrs == j) = [];
|
||||
if isempty(nbrs), continue; end
|
||||
angles = atan2(y(nbrs)-y(j), x(nbrs)-x(j));
|
||||
all_angles = [all_angles; angles(:)];
|
||||
end
|
||||
|
||||
% Histogram of bond angles over [-pi, pi]
|
||||
nbins = 36; % 10 degree bins
|
||||
edges = linspace(-pi, pi, nbins+1);
|
||||
counts = histcounts(all_angles, edges, 'Normalization','probability');
|
||||
|
||||
% Circular variance of bond angles
|
||||
R = abs(mean(exp(1i*all_angles)));
|
||||
circ_variance = 1 - R; % 0 means angles clustered; 1 means uniform
|
||||
|
||||
% 2) Positional correlation function g(r) estimate
|
||||
|
||||
% Compute pairwise distances
|
||||
coords = [x y];
|
||||
D = pdist(coords);
|
||||
max_r = min([Nx, Ny])/2; % max radius for g(r)
|
||||
|
||||
% Compute histogram of distances
|
||||
dr = 1; % bin width in pixels (adjust)
|
||||
r_edges = 0:dr:max_r;
|
||||
[counts_r, r_bins] = histcounts(D, r_edges);
|
||||
|
||||
% Normalize g(r) by ideal gas distribution (circle annulus area)
|
||||
r_centers = r_edges(1:end-1) + dr/2;
|
||||
area_annulus = pi*((r_edges(2:end)).^2 - (r_edges(1:end-1)).^2);
|
||||
density = Npoints / (Nx*Ny);
|
||||
|
||||
g_r = counts_r ./ (density * area_annulus * Npoints);
|
||||
|
||||
figure(1);
|
||||
clf
|
||||
set(gcf,'Position',[500 100 1250 500])
|
||||
t = tiledlayout(1, 3, 'TileSpacing', 'compact', 'Padding', 'compact'); % 1x4 grid
|
||||
|
||||
nexttile
|
||||
imshow(I, []);
|
||||
hold on;
|
||||
scatter(x, y, 30, 'w', 'filled');
|
||||
hold off;
|
||||
colormap(Helper.Colormaps.plasma()); % Example colormap for original image plot
|
||||
cbar1 = colorbar;
|
||||
cbar1.Label.Interpreter = 'latex';
|
||||
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);
|
||||
axis square;
|
||||
|
||||
% --- Plot angular distribution ---
|
||||
nexttile
|
||||
polarhistogram(all_angles, nbins, 'Normalization', 'probability');
|
||||
title(sprintf('Bond angle distribution (circ. variance = %.3f)', circ_variance));
|
||||
|
||||
% --- Plot g(r) ---
|
||||
nexttile
|
||||
plot(r_centers, g_r, 'LineWidth', 2);
|
||||
xlabel('r (pixels)');
|
||||
ylabel('g(r)');
|
||||
title('Radial Pair Correlation Function g(r)');
|
||||
grid on;
|
||||
axis square
|
||||
|
||||
|
||||
end
|
||||
63
Data-Analyzer/execution_scripts.m
Normal file
63
Data-Analyzer/execution_scripts.m
Normal file
@ -0,0 +1,63 @@
|
||||
%%
|
||||
|
||||
Data = load('D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_712500/Run_000/psi_gs.mat','psi','Params','Transf','Observ');
|
||||
|
||||
Params = Data.Params;
|
||||
Transf = Data.Transf;
|
||||
Observ = Data.Observ;
|
||||
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
|
||||
% Axes scaling and coordinates in micrometers
|
||||
x = Transf.x * Params.l0 * 1e6;
|
||||
y = Transf.y * Params.l0 * 1e6;
|
||||
z = Transf.z * Params.l0 * 1e6;
|
||||
dz = z(2)-z(1);
|
||||
% Compute probability density |psi|^2
|
||||
n = abs(psi).^2;
|
||||
nxy = squeeze(trapz(n*dz,3));
|
||||
IMG = nxy;
|
||||
|
||||
extractHexaticOrder(IMG)
|
||||
|
||||
%%
|
||||
|
||||
Data = load('D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_712500/Run_000/psi_gs.mat','psi','Params','Transf','Observ');
|
||||
% Data = load('D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta40/HighN/aS_9.562000e+01_theta_040_phi_000_N_508333/Run_000/psi_gs.mat','psi','Params','Transf','Observ');
|
||||
|
||||
Params = Data.Params;
|
||||
Transf = Data.Transf;
|
||||
Observ = Data.Observ;
|
||||
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
|
||||
% Axes scaling and coordinates in micrometers
|
||||
x = Transf.x * Params.l0 * 1e6;
|
||||
y = Transf.y * Params.l0 * 1e6;
|
||||
z = Transf.z * Params.l0 * 1e6;
|
||||
dz = z(2)-z(1);
|
||||
% Compute probability density |psi|^2
|
||||
n = abs(psi).^2;
|
||||
nxy = squeeze(trapz(n*dz,3));
|
||||
IMG = nxy;
|
||||
|
||||
%
|
||||
|
||||
[psi, ratio, N] = computeBondOrderParameters(IMG);
|
||||
|
||||
fprintf('Points: %d\n⟨|ψ₂|⟩ = %.3f, ⟨|ψ₄|⟩ = %.3f, ⟨|ψ₆|⟩ = %.3f\n', N, psi.psi2, psi.psi4, psi.psi6);
|
||||
fprintf('(⟨|ψ₆|⟩ / ⟨|ψ₂|⟩) = %.3f\n', ratio);
|
||||
@ -198,150 +198,6 @@ set([hXLabel, hYLabel], 'FontSize', 14)
|
||||
set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title
|
||||
grid on;
|
||||
|
||||
%% Extract g2 from simulation data
|
||||
|
||||
Data = load('E:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_712500/Run_000/psi_gs.mat','psi','Params','Transf','Observ');
|
||||
|
||||
% Data = load('E:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta40/HighN/aS_9.562000e+01_theta_040_phi_000_N_508333/Run_000/psi_gs.mat','psi','Params','Transf','Observ');
|
||||
|
||||
Params = Data.Params;
|
||||
Transf = Data.Transf;
|
||||
Observ = Data.Observ;
|
||||
|
||||
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
|
||||
|
||||
% Axes scaling and coordinates in micrometers
|
||||
x = Transf.x * Params.l0 * 1e6;
|
||||
y = Transf.y * Params.l0 * 1e6;
|
||||
z = Transf.z * Params.l0 * 1e6;
|
||||
|
||||
dx = x(2)-x(1); dy = y(2)-y(1); dz = z(2)-z(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 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
|
||||
|
||||
% Compute probability density |psi|^2
|
||||
n = abs(psi).^2;
|
||||
|
||||
nxz = squeeze(trapz(n*dy,2));
|
||||
nyz = squeeze(trapz(n*dx,1));
|
||||
nxy = squeeze(trapz(n*dz,3));
|
||||
|
||||
skipPreprocessing = true;
|
||||
skipMasking = true;
|
||||
skipIntensityThresholding = true;
|
||||
skipBinarization = true;
|
||||
|
||||
font = 'Bahnschrift';
|
||||
% Extract g2
|
||||
|
||||
N_bins = 90;
|
||||
Threshold = 75;
|
||||
Sigma = 2;
|
||||
|
||||
IMG = nxy;
|
||||
|
||||
[IMGFFT, IMGPR] = computeFourierTransform(IMG, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization);
|
||||
|
||||
[theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(IMGFFT, 10, 35, N_bins, Threshold, Sigma);
|
||||
|
||||
g2_all = zeros(1, N_bins); % Preallocate
|
||||
|
||||
for dtheta = 0:N_bins-1
|
||||
profile = S_theta;
|
||||
profile_shifted = circshift(profile, -dtheta, 2);
|
||||
|
||||
num = mean(profile .* profile_shifted);
|
||||
denom = mean(profile)^2;
|
||||
|
||||
g2_all(dtheta+1) = num / denom - 1;
|
||||
end
|
||||
|
||||
figure(2);
|
||||
clf
|
||||
set(gcf,'Position',[500 100 1000 800])
|
||||
t = tiledlayout(2, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 1x4 grid
|
||||
|
||||
% Display the cropped OD image
|
||||
nexttile
|
||||
plotxy = pcolor(x,y,IMG');
|
||||
set(plotxy, 'EdgeColor', 'none');
|
||||
cbar1 = colorbar;
|
||||
cbar1.Label.Interpreter = 'latex';
|
||||
colormap(gca, Helper.Colormaps.plasma())
|
||||
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 the power spectrum
|
||||
nexttile;
|
||||
imagesc(kx, ky, log(1 + abs(IMGFFT).^2));
|
||||
axis square;
|
||||
hcb = colorbar;
|
||||
colormap(gca, Helper.Colormaps.plasma())
|
||||
set(gca, 'FontSize', 14); % For tick labels only
|
||||
set(gca,'YDir','normal')
|
||||
hXLabel = xlabel('k_x', 'Interpreter', 'tex');
|
||||
hYLabel = ylabel('k_y', 'Interpreter', 'tex');
|
||||
hTitle = title('Power Spectrum - S(k_x,k_y)', 'Interpreter', 'tex');
|
||||
set([hXLabel, hYLabel], 'FontName', font)
|
||||
set([hXLabel, hYLabel], 'FontSize', 14)
|
||||
set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title
|
||||
|
||||
% Plot the angular distribution
|
||||
nexttile
|
||||
plot(theta_vals/pi, S_theta,'Linewidth',2);
|
||||
set(gca, 'FontSize', 14); % For tick labels only
|
||||
hXLabel = xlabel('\theta/\pi [rad]', 'Interpreter', 'tex');
|
||||
hYLabel = ylabel('Normalized magnitude (a.u.)', 'Interpreter', 'tex');
|
||||
hTitle = title('Angular Spectral Distribution - S(\theta)', 'Interpreter', 'tex');
|
||||
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
|
||||
|
||||
nexttile
|
||||
plot(theta_vals/pi, g2_all, 'o-', 'LineWidth', 1.2, 'MarkerSize', 5);
|
||||
set(gca, 'FontSize', 14);
|
||||
ylim([-1.5 3.0]); % Set y-axis limits here
|
||||
hXLabel = xlabel('$\delta\theta / \pi$', 'Interpreter', 'latex');
|
||||
hYLabel = ylabel('$g^{(2)}(\delta\theta)$', 'Interpreter', 'latex');
|
||||
hTitle = title('Autocorrelation', 'Interpreter', 'tex');
|
||||
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;
|
||||
|
||||
%% Helper Functions
|
||||
function [IMGFFT, IMGPR] = computeFourierTransform(I, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization)
|
||||
% computeFourierSpectrum - Computes the 2D Fourier power spectrum
|
||||
|
||||
98
Data-Analyzer/extractHexaticOrder.m
Normal file
98
Data-Analyzer/extractHexaticOrder.m
Normal file
@ -0,0 +1,98 @@
|
||||
function extractHexaticOrder(I)
|
||||
% Input:
|
||||
% I = 2D image matrix (grayscale, double)
|
||||
% Output:
|
||||
% Plots original image + hexatic order + histogram
|
||||
|
||||
%% 1) Detect dots (local maxima)
|
||||
I_smooth = imgaussfilt(I, 2);
|
||||
bw = imregionalmax(I_smooth);
|
||||
[y, x] = find(bw);
|
||||
|
||||
if length(x) < 6
|
||||
error('Too few detected dots (%d) for meaningful hexatic order.', length(x));
|
||||
end
|
||||
|
||||
%% 2) Compute hexatic order parameter psi6
|
||||
|
||||
% Delaunay triangulation to find neighbors
|
||||
tri = delaunay(x,y);
|
||||
N = length(x);
|
||||
neighbors = cell(N,1);
|
||||
for i = 1:size(tri,1)
|
||||
neighbors{tri(i,1)} = unique([neighbors{tri(i,1)} tri(i,2) tri(i,3)]);
|
||||
neighbors{tri(i,2)} = unique([neighbors{tri(i,2)} tri(i,1) tri(i,3)]);
|
||||
neighbors{tri(i,3)} = unique([neighbors{tri(i,3)} tri(i,1) tri(i,2)]);
|
||||
end
|
||||
|
||||
psi6_vals = zeros(N,1);
|
||||
for j = 1:N
|
||||
nbrs = neighbors{j};
|
||||
nbrs(nbrs == j) = [];
|
||||
if isempty(nbrs)
|
||||
psi6_vals(j) = 0;
|
||||
continue;
|
||||
end
|
||||
angles = atan2(y(nbrs)-y(j), x(nbrs)-x(j));
|
||||
psi6_vals(j) = abs(mean(exp(1i*6*angles)));
|
||||
end
|
||||
psi6_global = mean(psi6_vals);
|
||||
|
||||
%% 3) Plotting
|
||||
|
||||
figure('Name','Hexatic Order Analysis','NumberTitle','off', 'Units','normalized', 'Position',[0.2 0.2 0.6 0.7]);
|
||||
t = tiledlayout(2,2, 'Padding','none', 'TileSpacing','compact');
|
||||
|
||||
% Top left: Original image with detected dots (square)
|
||||
ax1 = nexttile(1);
|
||||
imshow(I, []);
|
||||
hold on;
|
||||
scatter(x, y, 30, 'w', 'filled');
|
||||
hold off;
|
||||
colormap(ax1, Helper.Colormaps.plasma()); % Example colormap for original image plot
|
||||
cbar1 = colorbar(ax1);
|
||||
cbar1.Label.Interpreter = 'latex';
|
||||
xlabel(ax1, '$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14);
|
||||
ylabel(ax1, '$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14);
|
||||
title(ax1, '$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14);
|
||||
axis(ax1, 'square');
|
||||
|
||||
% Top right: Histogram of local |psi6| (square)
|
||||
ax2 = nexttile(2);
|
||||
histogram(psi6_vals, 20, 'FaceColor', 'b');
|
||||
xlabel(ax2, '$|\psi_6|$', 'Interpreter', 'latex', 'FontSize', 14);
|
||||
ylabel(ax2, 'Count');
|
||||
title(ax2, 'Histogram of Local Hexatic Order');
|
||||
xlim(ax2, [0 1]);
|
||||
grid(ax2, 'on');
|
||||
axis(ax2, 'square');
|
||||
|
||||
% Bottom: Hexatic order magnitude on dots with bond orientation arrows (span 2 tiles)
|
||||
ax3 = nexttile([1 2]);
|
||||
scatter_handle = scatter(x, y, 50, psi6_vals, 'filled');
|
||||
colormap(ax3, 'jet'); % Different colormap for hexatic order magnitude
|
||||
cbar3 = colorbar(ax3);
|
||||
cbar3.Label.String = '|$\psi_6$|';
|
||||
cbar3.Label.Interpreter = 'latex';
|
||||
caxis(ax3, [0 1]);
|
||||
axis(ax3, 'equal');
|
||||
axis(ax3, 'tight');
|
||||
title(ax3, sprintf('Local Hexatic Order |\\psi_6|, Global = %.3f', psi6_global));
|
||||
hold(ax3, 'on');
|
||||
|
||||
for j = 1:N
|
||||
nbrs = neighbors{j};
|
||||
nbrs(nbrs == j) = [];
|
||||
if isempty(nbrs)
|
||||
continue;
|
||||
end
|
||||
angles = atan2(y(nbrs)-y(j), x(nbrs)-x(j));
|
||||
psi6_complex = mean(exp(1i*6*angles));
|
||||
local_angle = angle(psi6_complex)/6;
|
||||
arrow_length = 5;
|
||||
quiver(ax3, x(j), y(j), arrow_length*cos(local_angle), arrow_length*sin(local_angle), ...
|
||||
'k', 'MaxHeadSize', 2, 'AutoScale', 'off');
|
||||
end
|
||||
hold(ax3, 'off');
|
||||
|
||||
end
|
||||
@ -6,7 +6,7 @@ groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axi
|
||||
|
||||
folderPath = "D:/Data - Experiment/2025/05/22/";
|
||||
|
||||
run = '0078';
|
||||
run = '0079';
|
||||
|
||||
folderPath = strcat(folderPath, run);
|
||||
|
||||
@ -25,6 +25,8 @@ scan_parameter = 'rot_mag_fin_pol_angle';
|
||||
scan_parameter_text = 'Angle = ';
|
||||
% scan_parameter_text = 'BField = ';
|
||||
|
||||
savefodlerPath = 'D:/Results - Experiment/B2.45G/';
|
||||
savefileName = 'StripesToDroplets.mat';
|
||||
font = 'Bahnschrift';
|
||||
|
||||
skipPreprocessing = true;
|
||||
@ -96,6 +98,7 @@ end
|
||||
%% Run Fourier analysis over images
|
||||
|
||||
fft_imgs = cell(1, nimgs);
|
||||
spectral_contrast = zeros(1, nimgs);
|
||||
spectral_weight = zeros(1, nimgs);
|
||||
|
||||
N_bins = 180;
|
||||
@ -197,6 +200,7 @@ for k = 1:N_shots
|
||||
|
||||
% Plot the angular distribution
|
||||
nexttile
|
||||
spectral_contrast(k) = computeSpectralContrast(fft_imgs{k}, 10, 25, Threshold);
|
||||
[theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(fft_imgs{k}, 10, 20, N_bins, Threshold, Sigma);
|
||||
spectral_weight(k) = trapz(theta_vals, S_theta);
|
||||
plot(theta_vals/pi, S_theta,'Linewidth',2);
|
||||
@ -220,25 +224,51 @@ end
|
||||
% Close the video file
|
||||
close(videoFile);
|
||||
|
||||
%% Track spectral weight across the transition
|
||||
%% Track across the transition
|
||||
|
||||
% Assuming scan_parameter_values and spectral_weight are column vectors (or row vectors of same length)
|
||||
[unique_scan_parameter_values, ~, idx] = unique(scan_parameter_values);
|
||||
|
||||
% Preallocate arrays
|
||||
mean_sf = zeros(size(unique_scan_parameter_values));
|
||||
stderr_sf = zeros(size(unique_scan_parameter_values));
|
||||
mean_sc = zeros(size(unique_scan_parameter_values));
|
||||
stderr_sc = zeros(size(unique_scan_parameter_values));
|
||||
|
||||
% Loop through each unique theta and compute mean and standard error
|
||||
for i = 1:length(unique_scan_parameter_values)
|
||||
group_vals = spectral_weight(idx == i);
|
||||
mean_sf(i) = mean(group_vals);
|
||||
stderr_sf(i) = std(group_vals) / sqrt(length(group_vals)); % standard error = std / sqrt(N)
|
||||
group_vals = spectral_contrast(idx == i);
|
||||
mean_sc(i) = mean(group_vals);
|
||||
stderr_sc(i) = std(group_vals) / sqrt(length(group_vals)); % standard error = std / sqrt(N)
|
||||
end
|
||||
|
||||
figure(2);
|
||||
set(gcf,'Position',[100 100 950 750])
|
||||
errorbar(unique_scan_parameter_values, mean_sf, stderr_sf, 'o--', ...
|
||||
errorbar(unique_scan_parameter_values, mean_sc, stderr_sc, 'o--', ...
|
||||
'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5);
|
||||
set(gca, 'FontSize', 14); % For tick labels only
|
||||
hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex');
|
||||
% hXLabel = xlabel('B_z (G)', 'Interpreter', 'tex');
|
||||
hYLabel = ylabel('Spectral Contrast', 'Interpreter', 'tex');
|
||||
hTitle = title('Change across transition', 'Interpreter', 'tex');
|
||||
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
|
||||
|
||||
|
||||
% Preallocate arrays
|
||||
mean_sw = zeros(size(unique_scan_parameter_values));
|
||||
stderr_sw = zeros(size(unique_scan_parameter_values));
|
||||
|
||||
% Loop through each unique theta and compute mean and standard error
|
||||
for i = 1:length(unique_scan_parameter_values)
|
||||
group_vals = spectral_weight(idx == i);
|
||||
mean_sw(i) = mean(group_vals);
|
||||
stderr_sw(i) = std(group_vals) / sqrt(length(group_vals)); % standard error = std / sqrt(N)
|
||||
end
|
||||
|
||||
figure(3);
|
||||
set(gcf,'Position',[100 100 950 750])
|
||||
errorbar(unique_scan_parameter_values, mean_sw, stderr_sw, 'o--', ...
|
||||
'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5);
|
||||
set(gca, 'FontSize', 14); % For tick labels only
|
||||
hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex');
|
||||
@ -250,6 +280,8 @@ set([hXLabel, hYLabel], 'FontSize', 14)
|
||||
set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title
|
||||
grid on
|
||||
|
||||
save([savefolderPath savefileName], 'unique_scan_parameter_values', 'mean_sc', 'stderr_sc', 'mean_sw', 'stderr_sw');
|
||||
|
||||
%% k-means Clustering
|
||||
|
||||
% Reshape to column vector
|
||||
@ -433,6 +465,37 @@ function [theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(IM
|
||||
S_theta = S_theta / max(S_theta);
|
||||
end
|
||||
|
||||
function contrast = computeSpectralContrast(IMGFFT, r_min, r_max, threshold)
|
||||
% Apply threshold to isolate strong peaks
|
||||
IMGFFT(IMGFFT < threshold) = 0;
|
||||
|
||||
% Prepare polar coordinates
|
||||
[ny, nx] = size(IMGFFT);
|
||||
[X, Y] = meshgrid(1:nx, 1:ny);
|
||||
cx = ceil(nx/2);
|
||||
cy = ceil(ny/2);
|
||||
R = sqrt((X - cx).^2 + (Y - cy).^2);
|
||||
|
||||
% Ring region (annulus) mask
|
||||
ring_mask = (R >= r_min) & (R <= r_max);
|
||||
|
||||
% Squared magnitude in the ring
|
||||
ring_power = abs(IMGFFT).^2 .* ring_mask;
|
||||
|
||||
% Maximum power in the ring
|
||||
ring_max = max(ring_power(:));
|
||||
|
||||
% Power at the DC component
|
||||
dc_power = abs(IMGFFT(cy, cx))^2;
|
||||
|
||||
% Avoid division by zero
|
||||
if dc_power == 0
|
||||
contrast = Inf; % or NaN or 0, depending on how you want to handle this
|
||||
else
|
||||
contrast = ring_max / dc_power;
|
||||
end
|
||||
end
|
||||
|
||||
function ret = getBkgOffsetFromCorners(img, x_fraction, y_fraction)
|
||||
% image must be a 2D numerical array
|
||||
[dim1, dim2] = size(img);
|
||||
|
||||
164
Data-Analyzer/matchTFRadii.m
Normal file
164
Data-Analyzer/matchTFRadii.m
Normal file
@ -0,0 +1,164 @@
|
||||
% ----------------------------
|
||||
% Experimental data
|
||||
% ----------------------------
|
||||
PixelSize = 5.86; % microns
|
||||
|
||||
AtomNumbers = [9.14950197, 8.6845907 , 8.82521245, 8.5899089 , 8.21675841, ...
|
||||
8.96234044, 8.8636914 , 8.70332154, 8.82930706, 8.91869919, ...
|
||||
8.58553165, 8.73391981, 8.71943552, 8.11717678, 8.59490351, ...
|
||||
8.57514491, 8.81628891, 8.37211343, 8.76077699, 8.71297796, ...
|
||||
9.17634469, 8.81424285, 8.61176745, 8.40555897, 8.97137861, ...
|
||||
8.88393124, 8.66625724, 8.30688943, 9.02338201, 8.57729816, ...
|
||||
8.50333458, 8.67617084, 8.8936879 , 9.02031475, 8.98459233, ...
|
||||
8.76525048, 8.76801503, 8.58302559, 8.4617431 , 8.74479855, ...
|
||||
8.83882896, 8.69091377, 8.79282459, 8.51785483, 8.75629649, ...
|
||||
8.58994308, 8.36816564, 9.2429294 , 8.6583425 , 8.55827961];
|
||||
|
||||
EstimatedAtomNumber = mean(AtomNumbers) * 1E4;
|
||||
|
||||
TF_Radii_X_pixels = [48.44308968, 46.01326593, 46.45950681, 46.41644117, 45.56176919, ...
|
||||
46.60816438, 46.85307478, 47.61086543, 46.66687703, 46.17986721, ...
|
||||
46.67877165, 46.07789481, 46.42285497, 46.22167708, 45.95144492, ...
|
||||
47.05400117, 49.03005788, 45.84588659, 46.85742777, 45.9824117 , ...
|
||||
47.14731188, 47.4984484 , 45.9055646 , 47.31804553, 47.52321888, ...
|
||||
47.76823968, 46.459749 , 45.4498851 , 45.38339308, 46.68736642, ...
|
||||
45.76607233, 48.1796053 , 46.94291541, 47.54092708, 48.26130406, ...
|
||||
47.44092616, 48.73463214, 46.39356452, 48.74120217, 45.57014182, ...
|
||||
47.56467835, 46.62867035, 46.62322802, 46.03032919, 44.78559832, ...
|
||||
46.31282562, 46.83537518, 47.68015029, 47.71093571, 47.34079816];
|
||||
|
||||
TF_Radii_Y_pixels = [113.52610841, 113.68862761, 112.84031747, 114.22062324, ...
|
||||
112.45378285, 114.53863928, 111.39181472, 112.67024271, ...
|
||||
113.65387448, 113.57576769, 110.22576589, 110.45091803, ...
|
||||
109.97966067, 112.84785553, 109.3836049 , 111.22290862, ...
|
||||
111.17028727, 110.71088554, 111.72973603, 112.39623635, ...
|
||||
113.18160954, 112.00016346, 109.66542778, 111.98705097, ...
|
||||
112.35983901, 110.21703075, 112.14565939, 111.2029942 , ...
|
||||
110.74296 , 112.56607054, 112.58015318, 111.93031032, ...
|
||||
111.59774288, 112.30723266, 112.79543793, 111.08288891, ...
|
||||
113.85269603, 111.77349741, 113.58639434, 111.28694353, ...
|
||||
112.1993445 , 111.72215918, 111.93271101, 112.17593036, ...
|
||||
110.82246602, 113.61806907, 114.13693144, 114.27245731, ...
|
||||
114.24223538, 112.61704196];
|
||||
|
||||
% ----------------------------
|
||||
% Load simulation data for fitting
|
||||
% ----------------------------
|
||||
% baseDir = 'D:\Results - Numerics\Data_Full3D\PhaseDiagram\TFRadii\LowN';
|
||||
baseDir = 'C:\Users\Karthik\Documents\GitRepositories\Calculations\Estimations\ThomasFermiRadius';
|
||||
refData = load(fullfile(baseDir, 'TFFermi_Theta0.mat'));
|
||||
|
||||
% ----------------------------
|
||||
% Find simulation values at the estimated atom number
|
||||
% ----------------------------
|
||||
[~, idxClosest] = min(abs(refData.NUM_ATOMS_LIST - EstimatedAtomNumber));
|
||||
|
||||
if ndims(refData.TF_Radii) > 2
|
||||
TF_Radii = squeeze(refData.TF_Radii);
|
||||
TF_X_target = TF_Radii(idxClosest, 1);
|
||||
TF_Y_target = TF_Radii(idxClosest, 2);
|
||||
else
|
||||
TF_X_target = refData.TF_Radii(idxClosest, 1);
|
||||
TF_Y_target = refData.TF_Radii(idxClosest, 2);
|
||||
end
|
||||
|
||||
fprintf('Target radii from simulation at N = %.0f:\n', refData.NUM_ATOMS_LIST(idxClosest));
|
||||
fprintf('TF_X_target = %.2f µm\n', TF_X_target);
|
||||
fprintf('TF_Y_target = %.2f µm\n', TF_Y_target);
|
||||
|
||||
% ----------------------------
|
||||
% Find optimal magnification
|
||||
% ----------------------------
|
||||
errorFunc = @(mag) ...
|
||||
(mean(TF_Radii_X_pixels)*PixelSize/mag - TF_X_target)^2 + ...
|
||||
(mean(TF_Radii_Y_pixels)*PixelSize/mag - TF_Y_target)^2;
|
||||
|
||||
Magnification = fminbnd(errorFunc, 10, 50);
|
||||
fprintf('Best-fit magnification: %.4f\n', Magnification);
|
||||
|
||||
% ----------------------------
|
||||
% Convert to real space and get stats
|
||||
% ----------------------------
|
||||
TF_Radii_X_Real = TF_Radii_X_pixels * PixelSize / Magnification;
|
||||
TF_Radii_Y_Real = TF_Radii_Y_pixels * PixelSize / Magnification;
|
||||
|
||||
Avg_X = mean(TF_Radii_X_Real);
|
||||
Avg_Y = mean(TF_Radii_Y_Real);
|
||||
Std_X = std(TF_Radii_X_Real);
|
||||
Std_Y = std(TF_Radii_Y_Real);
|
||||
|
||||
fprintf('TF Radius X = %.2f ± %.2f µm\n', Avg_X, Std_X);
|
||||
fprintf('TF Radius Y = %.2f ± %.2f µm\n', Avg_Y, Std_Y);
|
||||
|
||||
% ----------------------------
|
||||
% Plotting
|
||||
% ----------------------------
|
||||
fileList = {'TFFermi_Theta0.mat'};
|
||||
thetaLabels = {'\theta = 0^\circ', '\theta = 20^\circ', '\theta = 40^\circ'};
|
||||
|
||||
fig = figure(1); clf;
|
||||
set(gcf,'Position', [100, 100, 1200, 500])
|
||||
t = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact');
|
||||
colors = lines(length(fileList));
|
||||
legendEntries = cell(1, length(fileList));
|
||||
|
||||
for j = 1:length(fileList)
|
||||
data = load(fullfile(baseDir, fileList{j}));
|
||||
|
||||
aS = data.SCATTERING_LENGTH_RANGE;
|
||||
NUM_ATOMS_LIST = data.NUM_ATOMS_LIST;
|
||||
if ndims(data.TF_Radii) > 2
|
||||
TF_Radii = squeeze(data.TF_Radii);
|
||||
else
|
||||
TF_Radii = data.TF_Radii;
|
||||
end
|
||||
|
||||
legendEntries{j} = sprintf('%s, a_s = %.2f a_0', thetaLabels{j}, aS);
|
||||
|
||||
% Rx
|
||||
nexttile(1);
|
||||
plot(NUM_ATOMS_LIST, TF_Radii(:,1), '-', ...
|
||||
'Color', colors(j,:), 'LineWidth', 1.5, ...
|
||||
'DisplayName', legendEntries{j}); hold on;
|
||||
|
||||
% Ry
|
||||
nexttile(2);
|
||||
plot(NUM_ATOMS_LIST, TF_Radii(:,2), '-', ...
|
||||
'Color', colors(j,:), 'LineWidth', 1.5, ...
|
||||
'DisplayName', legendEntries{j}); hold on;
|
||||
end
|
||||
|
||||
% ----------------------------
|
||||
% Add experimental point w/ error bars and annotation
|
||||
% ----------------------------
|
||||
% TF Radius X
|
||||
nexttile(1);
|
||||
errorbar(EstimatedAtomNumber, Avg_X, Std_X, ...
|
||||
'd', 'MarkerSize', 8, 'MarkerFaceColor', [0.2 0.2 0.8], ...
|
||||
'MarkerEdgeColor', 'k', 'Color', 'k', 'LineWidth', 1.2, ...
|
||||
'DisplayName', '\theta = 0^\circ, Experimental Value');
|
||||
|
||||
% TF Radius Y
|
||||
nexttile(2);
|
||||
errorbar(EstimatedAtomNumber, Avg_Y, Std_Y, ...
|
||||
'd', 'MarkerSize', 8, 'MarkerFaceColor', [0.2 0.2 0.8], ...
|
||||
'MarkerEdgeColor', 'k', 'Color', 'k', 'LineWidth', 1.2, ...
|
||||
'DisplayName', '\theta = 0^\circ, Experimental Value');
|
||||
|
||||
% ----------------------------
|
||||
% Finalize
|
||||
% ----------------------------
|
||||
nexttile(1);
|
||||
xlabel('Number of Atoms', 'FontSize', 16);
|
||||
ylabel('TF Radius - X ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16);
|
||||
legend('FontSize', 12, 'Interpreter', 'tex', 'Location', 'bestoutside');
|
||||
axis square; grid on;
|
||||
|
||||
nexttile(2);
|
||||
xlabel('Number of Atoms', 'FontSize', 16);
|
||||
ylabel('TF Radius - Y ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16);
|
||||
legend('FontSize', 12, 'Interpreter', 'tex', 'Location', 'bestoutside');
|
||||
axis square; grid on;
|
||||
|
||||
sgtitle('[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz', ...
|
||||
'Interpreter', 'tex', 'FontSize', 18);
|
||||
296
Dipolar-Gas-Simulator/+Scripts/conductFourierAnalysis.m
Normal file
296
Dipolar-Gas-Simulator/+Scripts/conductFourierAnalysis.m
Normal file
@ -0,0 +1,296 @@
|
||||
function [spectral_weight, g2, theta_vals] = conductFourierAnalysis(folder_path, run_index, N_bins, Threshold, Sigma, SuppressPlotFlag)
|
||||
|
||||
arguments
|
||||
folder_path (1,:) char
|
||||
run_index (1,:) {mustBeNumeric,mustBeReal}
|
||||
N_bins (1,:) {mustBeNumeric,mustBeReal}
|
||||
Threshold (1,:) {mustBeNumeric,mustBeReal}
|
||||
Sigma (1,:) {mustBeNumeric,mustBeReal}
|
||||
SuppressPlotFlag (1,:) logical = true
|
||||
end
|
||||
|
||||
set(0,'defaulttextInterpreter','latex')
|
||||
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
|
||||
|
||||
% Load data
|
||||
Data = load(fullfile(fullfile(folder_path, sprintf('Run_%03i', run_index)), 'psi_gs.mat'), 'psi', 'Transf', 'Observ', 'Params');
|
||||
|
||||
Params = Data.Params;
|
||||
Transf = Data.Transf;
|
||||
Observ = Data.Observ;
|
||||
|
||||
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
|
||||
|
||||
alpha = Params.theta;
|
||||
|
||||
% Axes scaling and coordinates in micrometers
|
||||
x = Transf.x * Params.l0 * 1e6;
|
||||
y = Transf.y * Params.l0 * 1e6;
|
||||
z = Transf.z * Params.l0 * 1e6;
|
||||
|
||||
dz = z(2)-z(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 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
|
||||
|
||||
% Compute probability density |psi|^2
|
||||
n = abs(psi).^2;
|
||||
|
||||
nxy = squeeze(trapz(n*dz,3));
|
||||
|
||||
skipPreprocessing = true;
|
||||
skipMasking = true;
|
||||
skipIntensityThresholding = true;
|
||||
skipBinarization = true;
|
||||
|
||||
%% Extract Spectral Weight and g2
|
||||
|
||||
IMG = nxy;
|
||||
|
||||
[IMGFFT, ~] = computeFourierTransform(IMG, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization);
|
||||
|
||||
[theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(IMGFFT, 10, 35, N_bins, Threshold, Sigma);
|
||||
|
||||
spectral_weight = trapz(theta_vals, S_theta);
|
||||
|
||||
g2 = zeros(1, N_bins); % Preallocate
|
||||
|
||||
for dtheta = 0:N_bins-1
|
||||
profile = S_theta;
|
||||
profile_shifted = circshift(profile, -dtheta, 2);
|
||||
|
||||
num = mean(profile .* profile_shifted);
|
||||
denom = mean(profile)^2;
|
||||
|
||||
g2(dtheta+1) = num / denom - 1;
|
||||
end
|
||||
|
||||
if ~SuppressPlotFlag
|
||||
figure(1);
|
||||
clf
|
||||
set(gcf,'Position',[500 100 1000 800])
|
||||
t = tiledlayout(2, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 1x4 grid
|
||||
font = 'Bahnschrift';
|
||||
% Display the cropped OD image
|
||||
ax1 = nexttile;
|
||||
plotxy = pcolor(x,y,IMG');
|
||||
set(plotxy, 'EdgeColor', 'none');
|
||||
% 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-right corner (normalized axis coordinates)
|
||||
hText = text(1 - x_offset, 1 - y_offset, ...
|
||||
['\alpha = ', num2str(rad2deg(alpha), '%.1f'), '^\circ'], ...
|
||||
'Color', 'white', 'FontWeight', 'bold', ...
|
||||
'Interpreter', 'tex', 'FontSize', 20, ...
|
||||
'Units', 'normalized', ...
|
||||
'HorizontalAlignment', 'right', ...
|
||||
'VerticalAlignment', 'top');
|
||||
axis square;
|
||||
hcb = colorbar;
|
||||
hcb.Label.Interpreter = 'latex';
|
||||
colormap(gca, Helper.Colormaps.plasma())
|
||||
set(gca, 'FontSize', 14); % For tick labels only
|
||||
hL = ylabel(hcb, 'Column Density');
|
||||
hXLabel = xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14);
|
||||
hYLabel = ylabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14);
|
||||
hTitle = title('$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14) ;
|
||||
set(hText, 'FontName', font)
|
||||
set([hXLabel, hYLabel, hL], 'FontSize', 14)
|
||||
set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title
|
||||
|
||||
|
||||
% Plot the power spectrum
|
||||
nexttile;
|
||||
imagesc(kx, ky, log(1 + abs(IMGFFT).^2));
|
||||
axis square;
|
||||
hcb = colorbar;
|
||||
colormap(gca, Helper.Colormaps.plasma())
|
||||
set(gca, 'FontSize', 14); % For tick labels only
|
||||
set(gca,'YDir','normal')
|
||||
hXLabel = xlabel('k_x', 'Interpreter', 'tex');
|
||||
hYLabel = ylabel('k_y', 'Interpreter', 'tex');
|
||||
hTitle = title('Power Spectrum - S(k_x,k_y)', 'Interpreter', 'tex');
|
||||
set([hXLabel, hYLabel], 'FontName', font)
|
||||
set([hXLabel, hYLabel], 'FontSize', 14)
|
||||
set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title
|
||||
|
||||
% Plot the angular distribution
|
||||
nexttile
|
||||
plot(theta_vals/pi, S_theta,'Linewidth',2);
|
||||
set(gca, 'FontSize', 14); % For tick labels only
|
||||
hXLabel = xlabel('\theta/\pi [rad]', 'Interpreter', 'tex');
|
||||
hYLabel = ylabel('Normalized magnitude (a.u.)', 'Interpreter', 'tex');
|
||||
hTitle = title('Angular Spectral Distribution - S(\theta)', 'Interpreter', 'tex');
|
||||
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
|
||||
|
||||
nexttile
|
||||
plot(theta_vals/pi, g2, 'o-', 'LineWidth', 1.2, 'MarkerSize', 5);
|
||||
set(gca, 'FontSize', 14);
|
||||
ylim([-1.5 3.0]); % Set y-axis limits here
|
||||
hXLabel = xlabel('$\delta\theta / \pi$', 'Interpreter', 'latex');
|
||||
hYLabel = ylabel('$g^{(2)}(\delta\theta)$', 'Interpreter', 'latex');
|
||||
hTitle = title('Autocorrelation', 'Interpreter', 'tex');
|
||||
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;
|
||||
end
|
||||
end
|
||||
|
||||
%% Helper Functions
|
||||
function [IMGFFT, IMGPR] = computeFourierTransform(I, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization)
|
||||
% computeFourierSpectrum - Computes the 2D Fourier power spectrum
|
||||
% of binarized and enhanced lattice image features, with optional central mask.
|
||||
%
|
||||
% Inputs:
|
||||
% I - Grayscale or RGB image matrix
|
||||
%
|
||||
% Output:
|
||||
% F_mag - 2D Fourier power spectrum (shifted)
|
||||
|
||||
if ~skipPreprocessing
|
||||
% Preprocessing: Denoise
|
||||
filtered = imgaussfilt(I, 10);
|
||||
IMGPR = I - filtered; % adjust sigma as needed
|
||||
else
|
||||
IMGPR = I;
|
||||
end
|
||||
|
||||
if ~skipMasking
|
||||
[rows, cols] = size(IMGPR);
|
||||
[X, Y] = meshgrid(1:cols, 1:rows);
|
||||
% Elliptical mask parameters
|
||||
cx = cols / 2;
|
||||
cy = rows / 2;
|
||||
|
||||
% Shifted coordinates
|
||||
x = X - cx;
|
||||
y = Y - cy;
|
||||
|
||||
% Ellipse semi-axes
|
||||
rx = 0.4 * cols;
|
||||
ry = 0.2 * rows;
|
||||
|
||||
% Rotation angle in degrees -> radians
|
||||
theta_deg = 30; % Adjust as needed
|
||||
theta = deg2rad(theta_deg);
|
||||
|
||||
% Rotated ellipse equation
|
||||
cos_t = cos(theta);
|
||||
sin_t = sin(theta);
|
||||
|
||||
x_rot = (x * cos_t + y * sin_t);
|
||||
y_rot = (-x * sin_t + y * cos_t);
|
||||
|
||||
ellipseMask = (x_rot.^2) / rx^2 + (y_rot.^2) / ry^2 <= 1;
|
||||
|
||||
% Apply cutout mask
|
||||
IMGPR = IMGPR .* ellipseMask;
|
||||
end
|
||||
|
||||
if ~skipIntensityThresholding
|
||||
% Apply global intensity threshold mask
|
||||
intensity_thresh = 0.20;
|
||||
intensity_mask = IMGPR > intensity_thresh;
|
||||
IMGPR = IMGPR .* intensity_mask;
|
||||
end
|
||||
|
||||
if ~skipBinarization
|
||||
% Adaptive binarization and cleanup
|
||||
IMGPR = imbinarize(IMGPR, 'adaptive', 'Sensitivity', 0.0);
|
||||
IMGPR = imdilate(IMGPR, strel('disk', 2));
|
||||
IMGPR = imerode(IMGPR, strel('disk', 1));
|
||||
IMGPR = imfill(IMGPR, 'holes');
|
||||
F = fft2(double(IMGPR)); % Compute 2D Fourier Transform
|
||||
IMGFFT = abs(fftshift(F))'; % Shift zero frequency to center
|
||||
else
|
||||
F = fft2(double(IMGPR)); % Compute 2D Fourier Transform
|
||||
IMGFFT = abs(fftshift(F))'; % Shift zero frequency to center
|
||||
end
|
||||
end
|
||||
|
||||
function [theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(IMGFFT, r_min, r_max, num_bins, threshold, sigma)
|
||||
% Apply threshold to isolate strong peaks
|
||||
IMGFFT(IMGFFT < threshold) = 0;
|
||||
|
||||
% Prepare polar coordinates
|
||||
[ny, nx] = size(IMGFFT);
|
||||
[X, Y] = meshgrid(1:nx, 1:ny);
|
||||
cx = ceil(nx/2);
|
||||
cy = ceil(ny/2);
|
||||
R = sqrt((X - cx).^2 + (Y - cy).^2);
|
||||
Theta = atan2(Y - cy, X - cx); % range [-pi, pi]
|
||||
|
||||
% Choose radial band
|
||||
radial_mask = (R >= r_min) & (R <= r_max);
|
||||
|
||||
% Initialize the angular structure factor array
|
||||
S_theta = zeros(1, num_bins); % Pre-allocate for 180 angle bins
|
||||
% Define the angle values for the x-axis
|
||||
theta_vals = linspace(0, pi, num_bins);
|
||||
|
||||
% Loop through each angle bin
|
||||
for i = 1:num_bins
|
||||
angle_start = (i-1) * pi / num_bins;
|
||||
angle_end = i * pi / num_bins;
|
||||
|
||||
% Define a mask for the given angle range
|
||||
angle_mask = (Theta >= angle_start & Theta < angle_end);
|
||||
|
||||
bin_mask = radial_mask & angle_mask;
|
||||
|
||||
% Extract the Fourier components for the given angle
|
||||
fft_angle = IMGFFT .* bin_mask;
|
||||
|
||||
% Integrate the Fourier components over the radius at the angle
|
||||
S_theta(i) = sum(sum(abs(fft_angle).^2)); % sum of squared magnitudes
|
||||
end
|
||||
|
||||
% Create a 1D Gaussian kernel
|
||||
half_width = ceil(3 * sigma);
|
||||
x = -half_width:half_width;
|
||||
gauss_kernel = exp(-x.^2 / (2 * sigma^2));
|
||||
gauss_kernel = gauss_kernel / sum(gauss_kernel); % normalize
|
||||
|
||||
% Apply convolution (circular padding to preserve periodicity)
|
||||
S_theta = conv([S_theta(end-half_width+1:end), S_theta, S_theta(1:half_width)], gauss_kernel, 'same');
|
||||
S_theta = S_theta(half_width+1:end-half_width); % crop back to original size
|
||||
|
||||
% Normalize to 1
|
||||
S_theta = S_theta / max(S_theta);
|
||||
end
|
||||
|
||||
@ -5,7 +5,15 @@ function [AveragePCD] = extractAveragePeakColumnDensity(folder_path, run_index,
|
||||
format long
|
||||
|
||||
% Load data
|
||||
Data = load(sprintf(horzcat(folder_path, '/Run_%03i/psi_gs.mat'),run_index),'psi','Params','Transf','Observ');
|
||||
filePath = sprintf(horzcat(folder_path, '/Run_%03i/psi_gs.mat'), run_index);
|
||||
|
||||
try
|
||||
Data = load(filePath, 'psi', 'Params', 'Transf', 'Observ');
|
||||
catch ME
|
||||
warning('Failed to load file: %s\n%s', filePath, ME.message);
|
||||
AveragePCD = NaN;
|
||||
return;
|
||||
end
|
||||
|
||||
Params = Data.Params;
|
||||
Transf = Data.Transf;
|
||||
|
||||
@ -1,85 +1,141 @@
|
||||
function [UCD] = extractAverageUnitCellDensity(folder_path, run_index, radius, minPeakHeight, SuppressPlotFlag)
|
||||
set(0,'defaulttextInterpreter','latex')
|
||||
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
|
||||
|
||||
format long
|
||||
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','Params','Transf','Observ');
|
||||
filePath = sprintf(horzcat(folder_path, '/Run_%03i/psi_gs.mat'), run_index);
|
||||
|
||||
try
|
||||
Data = load(filePath, 'psi', 'Params', 'Transf', 'Observ');
|
||||
catch ME
|
||||
warning('Failed to load file: %s\n%s', filePath, ME.message);
|
||||
UCD = NaN;
|
||||
return;
|
||||
end
|
||||
|
||||
Params = Data.Params;
|
||||
Transf = Data.Transf;
|
||||
Observ = Data.Observ;
|
||||
|
||||
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
|
||||
|
||||
% Axes scaling and coordinates in micrometers
|
||||
% Axes in micrometers
|
||||
x = Transf.x * Params.l0 * 1e6;
|
||||
y = Transf.y * Params.l0 * 1e6;
|
||||
z = Transf.z * Params.l0 * 1e6;
|
||||
|
||||
dz = z(2) - z(1);
|
||||
|
||||
% Compute probability density |psi|^2
|
||||
% Compute integrated density
|
||||
n = abs(psi).^2;
|
||||
nxy = squeeze(trapz(n * dz, 3));
|
||||
|
||||
state = nxy';
|
||||
|
||||
% Fourier spectrum for orientation detection
|
||||
F = fftshift(abs(fft2(state - mean(state(:)))));
|
||||
[nx, ny] = size(state);
|
||||
kx = linspace(-pi / (x(2) - x(1)), pi / (x(2) - x(1)), nx);
|
||||
ky = linspace(-pi / (y(2) - y(1)), pi / (y(2) - y(1)), ny);
|
||||
[KX, KY] = meshgrid(kx, ky);
|
||||
theta = atan2(KY, KX);
|
||||
|
||||
% Remove center/DC
|
||||
R = sqrt(KX.^2 + KY.^2);
|
||||
F(R < 0.1 * max(kx(:))) = 0;
|
||||
|
||||
% Angular histogram of power spectrum
|
||||
nbins = 180;
|
||||
angles = linspace(-pi, pi, nbins);
|
||||
powerAngular = zeros(1, nbins);
|
||||
for k = 1:nbins-1
|
||||
mask = theta >= angles(k) & theta < angles(k+1);
|
||||
powerAngular(k) = sum(F(mask), 'all');
|
||||
end
|
||||
powerAngular(end) = powerAngular(1); % wrap around
|
||||
|
||||
% Anisotropy measure
|
||||
peakRatio = max(powerAngular) / mean(powerAngular);
|
||||
|
||||
% Threshold to distinguish stripe vs droplet
|
||||
isStripe = peakRatio > 20;
|
||||
|
||||
if ~isStripe
|
||||
% === DROPLET MODE: Voronoi cell ===
|
||||
peaks = extractPeaks(state, radius, minPeakHeight);
|
||||
peakHeights = state(peaks);
|
||||
[row, col] = find(peaks);
|
||||
peakHeights = state(peaks);
|
||||
[~, maxIdx] = max(peakHeights);
|
||||
MaxPeakLocation = [row(maxIdx), col(maxIdx)];
|
||||
points = [col, row];
|
||||
|
||||
% Voronoi diagram of peak positions
|
||||
points = [col, row]; % Voronoi uses [x, y]
|
||||
[V, C] = voronoin(points);
|
||||
|
||||
% Voronoi cell of the max peak
|
||||
Vcell_indices = C{maxIdx};
|
||||
|
||||
% Plot the Voronoi cell
|
||||
if all(Vcell_indices > 0) && all(Vcell_indices <= size(V,1)) && ~any(isinf(V(Vcell_indices, 1)))
|
||||
vx = interp1(1:numel(x), x, V(Vcell_indices,1), 'linear', 'extrap');
|
||||
vy = interp1(1:numel(y), y, V(Vcell_indices,2), 'linear', 'extrap');
|
||||
% Close the polygon by repeating the first vertex
|
||||
vx(end+1) = vx(1);
|
||||
vy(end+1) = vy(1);
|
||||
% Compute area of the Voronoi cell polygon using the shoelace formula
|
||||
AreaOfVoronoiCell = polyarea(vx, vy); % Area of Voronoi cell around max peak in um^2
|
||||
AreaOfVoronoiCell = polyarea(vx, vy);
|
||||
|
||||
[X, Y] = meshgrid(x, y);
|
||||
inCellMask = inpolygon(X, Y, vx, vy);
|
||||
NumberOfParticles = sum(state(inCellMask));
|
||||
UCD = NumberOfParticles / AreaOfVoronoiCell;
|
||||
else
|
||||
warning('Voronoi cell for max peak is unbounded or invalid.');
|
||||
warning('Voronoi cell for max peak is invalid.');
|
||||
UCD = NaN;
|
||||
end
|
||||
|
||||
% Create grid points mesh
|
||||
[X, Y] = meshgrid(x, y); % Note: size(X) and size(Y) match size(state)
|
||||
else
|
||||
% === STRIPE MODE: Use rectangular unit cell aligned with stripe direction ===
|
||||
[~, idxMax] = max(F(:));
|
||||
stripeAngle = theta(idxMax); % radians
|
||||
|
||||
% Determine points inside Voronoi polygon
|
||||
inCellMask = inpolygon(X, Y, vx, vy);
|
||||
% Rotate image to align stripes horizontally
|
||||
rotatedState = imrotate(state, -rad2deg(stripeAngle), 'bilinear', 'crop');
|
||||
|
||||
% Sum all state values inside the Voronoi cell polygon
|
||||
NumberOfParticlesInVoronoiCell = sum(state(inCellMask));
|
||||
% Estimate vertical stripe spacing (stripe-normal direction)
|
||||
stripeProfileY = mean(rotatedState, 2) - mean(rotatedState(:));
|
||||
fftY = abs(fft(stripeProfileY));
|
||||
fftY(1:3) = 0;
|
||||
[~, fyIdx] = max(fftY(1:floor(end/2)));
|
||||
spacingY = round(numel(stripeProfileY) / fyIdx);
|
||||
|
||||
UCD = NumberOfParticlesInVoronoiCell / AreaOfVoronoiCell;
|
||||
% Estimate horizontal spacing (along-stripe periodicity)
|
||||
stripeProfileX = mean(rotatedState, 1) - mean(rotatedState(:));
|
||||
fftX = abs(fft(stripeProfileX));
|
||||
fftX(1:3) = 0;
|
||||
[~, fxIdx] = max(fftX(1:floor(end/2)));
|
||||
spacingX = round(numel(stripeProfileX) / fxIdx);
|
||||
|
||||
% Find stripe center (vertical)
|
||||
[~, centerY] = max(stripeProfileY);
|
||||
rowRange = max(1, centerY - round(0.5 * spacingY)) : ...
|
||||
min(size(rotatedState,1), centerY + round(0.5 * spacingY));
|
||||
|
||||
% Find repeat center along stripe direction (horizontal)
|
||||
[~, centerX] = max(stripeProfileX);
|
||||
colRange = max(1, centerX - round(0.5 * spacingX)) : ...
|
||||
min(size(rotatedState,2), centerX + round(0.5 * spacingX));
|
||||
|
||||
% Extract unit cell region
|
||||
unitCellRegion = rotatedState(rowRange, colRange);
|
||||
NumberOfParticles = sum(unitCellRegion(:));
|
||||
AreaOfUnitCell = numel(unitCellRegion) * abs(x(2)-x(1)) * abs(y(2)-y(1));
|
||||
UCD = NumberOfParticles / AreaOfUnitCell;
|
||||
end
|
||||
|
||||
|
||||
% Optional plot
|
||||
if ~SuppressPlotFlag
|
||||
figure(1);
|
||||
clf
|
||||
set(gcf,'Position', [100, 100, 900, 900])
|
||||
if ~isStripe
|
||||
plotxy = pcolor(x,y,state);
|
||||
set(plotxy, 'EdgeColor', 'none');
|
||||
hold on;
|
||||
plot(x(MaxPeakLocation(2)), y(MaxPeakLocation(1)), 'w+', 'MarkerSize', 8, 'LineWidth', 1.5);
|
||||
plot(vx, vy, 'w-', 'LineWidth', 1.5);
|
||||
plot(vx, vy, 'w-', 'LineWidth', 2);
|
||||
cbar1 = colorbar;
|
||||
cbar1.Label.Interpreter = 'latex';
|
||||
% cbar1.Ticks = []; % Disable the ticks
|
||||
@ -89,6 +145,31 @@ function [UCD] = extractAverageUnitCellDensity(folder_path, run_index, radius, m
|
||||
title(['UCD = ' num2str(UCD) ' \mum^{-2}'], ...
|
||||
'Interpreter', 'tex', 'FontSize', 16)
|
||||
set(gca, 'FontSize', 16); % For tick labels only
|
||||
else
|
||||
imagesc(rotatedState);
|
||||
axis image;
|
||||
cbar1 = colorbar;
|
||||
cbar1.Label.Interpreter = 'latex';
|
||||
colormap(gca, Helper.Colormaps.plasma())
|
||||
xlabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16)
|
||||
ylabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16)
|
||||
|
||||
% Title with UCD
|
||||
title(['UCD = ' num2str(UCD, '%.3f') ' $\mu$m$^{-2}$'], ...
|
||||
'Interpreter', 'latex', 'FontSize', 16)
|
||||
|
||||
% Highlight unit cell region
|
||||
rectX = colRange(1) - 0.5;
|
||||
rectY = rowRange(1) - 0.5;
|
||||
rectW = length(colRange);
|
||||
rectH = length(rowRange);
|
||||
rectangle('Position', [rectX, rectY, rectW, rectH], ...
|
||||
'EdgeColor', 'w', 'LineWidth', 2);
|
||||
title(sprintf('UCD = %.4f $\\mu$m$^{-2}$', UCD), ...
|
||||
'Interpreter', 'latex', 'FontSize', 16);
|
||||
|
||||
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
@ -48,7 +48,7 @@ function run_hybrid_worker(batchParams, batchIdx)
|
||||
OptionsStruct.MaxIterationsForGD = 15000;
|
||||
OptionsStruct.TimeStepSize = 1E-3; % in s
|
||||
OptionsStruct.MinimumTimeStepSize = 1E-6; % in s
|
||||
OptionsStruct.TimeCutOff = 2E5; % in s
|
||||
OptionsStruct.TimeCutOff = 1E6; % in s
|
||||
OptionsStruct.EnergyTolerance = 5E-10;
|
||||
OptionsStruct.ResidualTolerance = 1E-05;
|
||||
OptionsStruct.NoiseScaleFactor = 0.010;
|
||||
|
||||
@ -573,9 +573,14 @@ JobNumber = 0;
|
||||
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
|
||||
|
||||
%%
|
||||
SaveDirectory = 'D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_712500';
|
||||
SaveDirectory = 'D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_1733333';
|
||||
JobNumber = 0;
|
||||
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
|
||||
%%
|
||||
SaveDirectory = 'D:/Results - Numerics/Data_Full3D/PhaseTransition/STD/aS_9.562000e+01_theta_025_phi_000_N_500000';
|
||||
JobNumber = 0;
|
||||
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
|
||||
|
||||
|
||||
%% Identify and count droplets
|
||||
Radius = 2; % The radius within which peaks will be considered duplicates
|
||||
@ -593,7 +598,7 @@ JobNumber = 0;
|
||||
SuppressPlotFlag = false;
|
||||
AveragePCD = Scripts.extractAveragePeakColumnDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag);
|
||||
|
||||
%% Extract average unit cell density
|
||||
%% Extract average unit cell density - Droplets
|
||||
Radius = 2; % The radius within which peaks will be considered duplicates
|
||||
PeakThreshold = 3E3;
|
||||
SaveDirectory = 'D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_712500';
|
||||
@ -601,6 +606,19 @@ JobNumber = 0;
|
||||
SuppressPlotFlag = false;
|
||||
UCD = Scripts.extractAverageUnitCellDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag);
|
||||
|
||||
%% Extract average unit cell density - Stripes
|
||||
Radius = 2; % The radius within which peaks will be considered duplicates
|
||||
PeakThreshold = 3E3;
|
||||
SaveDirectory = 'D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_1325000';
|
||||
JobNumber = 0;
|
||||
SuppressPlotFlag = false;
|
||||
UCD = Scripts.extractAverageUnitCellDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag);
|
||||
|
||||
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];
|
||||
|
||||
%% Plot number of droplets
|
||||
% Parameters
|
||||
Radius = 2;
|
||||
@ -676,60 +694,166 @@ SuppressPlotFlag = true;
|
||||
|
||||
SCATTERING_LENGTH_RANGE = [95.62];
|
||||
|
||||
NUM_ATOMS_LIST = [50000 54545 59091 63636 68182 72727 77273 81818 86364 90909 95455 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];
|
||||
NUM_ATOMS_LIST_FULL = [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];
|
||||
|
||||
NUM_ATOMS_LIST_INSET = [50000 54545 59091 63636 68182 72727 77273 81818 86364 90909 95455];
|
||||
|
||||
% Prepare figure
|
||||
figure(1);
|
||||
clf;
|
||||
set(gcf,'Position', [100, 100, 1000, 700]);
|
||||
hold on;
|
||||
|
||||
% Color order for better visibility
|
||||
% Color order
|
||||
colors = lines(length(SCATTERING_LENGTH_RANGE));
|
||||
|
||||
% Store legend labels
|
||||
legendEntries = cell(1, length(SCATTERING_LENGTH_RANGE));
|
||||
% Store data for both sets
|
||||
AverageCDs_full = zeros(length(SCATTERING_LENGTH_RANGE), length(NUM_ATOMS_LIST_FULL));
|
||||
AverageCDs_inset = zeros(length(SCATTERING_LENGTH_RANGE), length(NUM_ATOMS_LIST_INSET));
|
||||
|
||||
% Main plot
|
||||
main_ax = axes;
|
||||
hold(main_ax, 'on');
|
||||
|
||||
% Loop over scattering lengths
|
||||
for j = 1:length(SCATTERING_LENGTH_RANGE)
|
||||
aS = SCATTERING_LENGTH_RANGE(j);
|
||||
|
||||
% Format scattering length in scientific notation with 6 decimal places
|
||||
aS_string = sprintf('%.6e', aS);
|
||||
|
||||
% Construct base directory for this aS
|
||||
baseDir = ['D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/aS_' ...
|
||||
baseDir = ['D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_' ...
|
||||
aS_string '_theta_000_phi_000_N_'];
|
||||
|
||||
% Preallocate results for this curve
|
||||
AverageCDs = zeros(size(NUM_ATOMS_LIST));
|
||||
|
||||
% Loop over all atom numbers
|
||||
for i = 1:length(NUM_ATOMS_LIST)
|
||||
N = NUM_ATOMS_LIST(i);
|
||||
% Full list
|
||||
for i = 1:length(NUM_ATOMS_LIST_FULL)
|
||||
N = NUM_ATOMS_LIST_FULL(i);
|
||||
SaveDirectory = [baseDir num2str(N)];
|
||||
|
||||
% Call your droplet counting function
|
||||
AverageCDs(i) = Scripts.extractAveragePeakColumnDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag);
|
||||
AverageCDs_full(j,i) = Scripts.extractAveragePeakColumnDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag);
|
||||
end
|
||||
|
||||
% Plot curve
|
||||
plot(NUM_ATOMS_LIST, AverageCDs, 'o-', ...
|
||||
baseDir = ['D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/LowN/aS_' ...
|
||||
aS_string '_theta_000_phi_000_N_'];
|
||||
|
||||
% Inset list
|
||||
for i = 1:length(NUM_ATOMS_LIST_INSET)
|
||||
N = NUM_ATOMS_LIST_INSET(i);
|
||||
SaveDirectory = [baseDir num2str(N)];
|
||||
AverageCDs_inset(j,i) = Scripts.extractAveragePeakColumnDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag);
|
||||
end
|
||||
|
||||
% Plot main curve
|
||||
x = NUM_ATOMS_LIST_FULL;
|
||||
y = AverageCDs_full(j,:);
|
||||
valid = ~isnan(y); % logical index of valid points
|
||||
plot(main_ax,x(valid), y(valid), 'o-', ...
|
||||
'Color', colors(j,:), 'LineWidth', 1.5);
|
||||
|
||||
% Store legend entry
|
||||
legendEntries{j} = ['a_s = ' num2str(aS) 'a_o'];
|
||||
end
|
||||
|
||||
% Finalize plot
|
||||
xlabel('Number of Atoms', 'Interpreter', 'tex', 'FontSize', 16);
|
||||
ylabel('Average Peak Column Density', 'Interpreter', 'tex', 'FontSize', 16);
|
||||
title(TitleString, 'Interpreter', 'tex', 'FontSize', 18);
|
||||
legend(legendEntries, 'Interpreter', 'tex', 'FontSize', 12, 'Location', 'bestoutside');
|
||||
legend(arrayfun(@(aS) sprintf('a_s = %.2f a_0', aS), SCATTERING_LENGTH_RANGE, ...
|
||||
'UniformOutput', false), ...
|
||||
'Interpreter', 'tex', 'FontSize', 12, 'Location', 'bestoutside');
|
||||
grid(main_ax, 'on');
|
||||
|
||||
% Inset plot
|
||||
inset_ax = axes('Position', [0.45 0.18 0.28 0.28]); % Normalized position [x y w h]
|
||||
box(inset_ax, 'on');
|
||||
hold(inset_ax, 'on');
|
||||
|
||||
for j = 1:length(SCATTERING_LENGTH_RANGE)
|
||||
plot(inset_ax, NUM_ATOMS_LIST_INSET, AverageCDs_inset(j,:), 'o-', ...
|
||||
'Color', colors(j,:), 'LineWidth', 1.2);
|
||||
end
|
||||
|
||||
set(inset_ax, 'FontSize', 8);
|
||||
title(inset_ax, 'Low-N', 'FontSize', 9);
|
||||
grid(inset_ax, 'on');
|
||||
xlabel(inset_ax, 'N', 'FontSize', 9);
|
||||
ylabel(inset_ax, 'CD', 'FontSize', 9);
|
||||
|
||||
%% Plot average unit cell density
|
||||
Radius = 2;
|
||||
PeakThreshold = 3E3;
|
||||
JobNumber = 0;
|
||||
SuppressPlotFlag = true; % Suppress plots during batch processing
|
||||
TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 0^\circ";
|
||||
|
||||
|
||||
SCATTERING_LENGTH_RANGE = [95.62];
|
||||
|
||||
NUM_ATOMS_LIST = [712500 916667 1120833 1325000 1529167 1733333 1937500 2141667 2345833 2550000 2754167 2958333 3162500 3366667 3570833];
|
||||
|
||||
UCD_values = zeros(length(SCATTERING_LENGTH_RANGE), length(NUM_ATOMS_LIST));
|
||||
|
||||
% Prepare figure
|
||||
figure(1);
|
||||
clf;
|
||||
set(gcf,'Position', [100, 100, 1000, 700]);
|
||||
hold on
|
||||
|
||||
% Color order
|
||||
colors = lines(length(SCATTERING_LENGTH_RANGE));
|
||||
|
||||
for j = 1:length(SCATTERING_LENGTH_RANGE)
|
||||
aS = SCATTERING_LENGTH_RANGE(j);
|
||||
aS_string = sprintf('%.6e', aS);
|
||||
|
||||
baseDir = ['D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_' ...
|
||||
aS_string '_theta_000_phi_000_N_'];
|
||||
|
||||
for i = 1:length(NUM_ATOMS_LIST)
|
||||
N = NUM_ATOMS_LIST(i);
|
||||
% Construct folder path for this N
|
||||
SaveDirectory = sprintf('%s%d', baseDir, N);
|
||||
try
|
||||
UCD_values(j,i) = Scripts.extractAverageUnitCellDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag);
|
||||
catch ME
|
||||
warning('Error processing N=%d: %s', N, ME.message);
|
||||
UCD_values(j,i) = NaN; % mark as NaN on error
|
||||
end
|
||||
end
|
||||
|
||||
x = NUM_ATOMS_LIST;
|
||||
y = UCD_values(j,:);
|
||||
|
||||
valid = ~isnan(y); % logical index of valid points
|
||||
plot(x(valid), y(valid), 'o-', 'Color', colors(j,:), 'LineWidth', 1.5);
|
||||
|
||||
end
|
||||
|
||||
xlabel('Number of Atoms', 'Interpreter', 'latex', 'FontSize', 16);
|
||||
ylabel('Unit Cell Density (UCD) [$\mu m^{-2}$]', 'Interpreter', 'latex', 'FontSize', 16);
|
||||
title(TitleString, 'Interpreter', 'tex', 'FontSize', 18);
|
||||
legend(arrayfun(@(aS) sprintf('a_s = %.2f a_0', aS), SCATTERING_LENGTH_RANGE, ...
|
||||
'UniformOutput', false), ...
|
||||
'Interpreter', 'tex', 'FontSize', 12, 'Location', 'bestoutside');
|
||||
set(gca, 'FontSize', 14);
|
||||
grid on;
|
||||
|
||||
% Physical constants
|
||||
PlanckConstant = 6.62607015E-34;
|
||||
PlanckConstantReduced = 6.62607015E-34/(2*pi);
|
||||
AtomicMassUnit = 1.660539066E-27;
|
||||
BohrMagneton = 9.274009994E-24;
|
||||
|
||||
% Dy specific constants
|
||||
Dy164Mass = 163.929174751*AtomicMassUnit;
|
||||
Dy164IsotopicAbundance = 0.2826;
|
||||
DyMagneticMoment = 9.93*BohrMagneton;
|
||||
|
||||
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
|
||||
nadd2s = 0.01:0.01:0.25;
|
||||
ppmum = nadd2s.*(1E12*add^2)^-1;
|
||||
%
|
||||
figure(2);
|
||||
clf;
|
||||
set(gcf,'Position', [100, 100, 850, 700]);
|
||||
hold on
|
||||
plot(nadd2s, ppmum, 'o-', 'LineWidth', 1.5)
|
||||
xlabel('na_{dd}^2', 'Interpreter', 'tex', 'FontSize', 16);
|
||||
ylabel('Unit Cell Density (UCD) [$\mu m^{-2}$]', 'Interpreter', 'latex', 'FontSize', 16);
|
||||
title("[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 0, 0, 72.4 ] Hz; \theta = 0^\circ", 'Interpreter', 'tex', 'FontSize', 18);
|
||||
set(gca, 'FontSize', 14);
|
||||
grid on;
|
||||
hold off;
|
||||
|
||||
%% Plot TF radii of unmodulated states
|
||||
% Parameters
|
||||
@ -1066,3 +1190,69 @@ else
|
||||
disp('The state is not modulated');
|
||||
end
|
||||
|
||||
%% Plot Spectral weight and g2 across transition for simulated data
|
||||
N_atoms = 5E5;
|
||||
a_s = 95.62;
|
||||
phi_deg = 0.0;
|
||||
alpha_vals = [0.0 5.0 10.0 15.0 20.0 25.0 30.0 35.0 40.0 45.0];
|
||||
|
||||
N_alpha = length(alpha_vals);
|
||||
N_bins = 90;
|
||||
Threshold = 75;
|
||||
Sigma = 2;
|
||||
|
||||
spectral_weights = zeros(size(alpha_vals));
|
||||
g2_all = zeros(N_alpha, N_bins);
|
||||
theta_vals_all = zeros(N_alpha, N_bins);
|
||||
|
||||
JobNumber = 1;
|
||||
SuppressPlotFlag = true;
|
||||
|
||||
for i = 1:N_alpha
|
||||
folderName = sprintf('aS_%03d_theta_%03d_phi_%03d_N_%d', a_s, alpha_vals(i), phi_deg, N_atoms);
|
||||
SaveDirectory = fullfile('D:/Results - Numerics/Data_Full3D/PhaseTransition/DTS/', folderName);
|
||||
[spectral_weight, g2, theta_vals] = Scripts.conductFourierAnalysis(SaveDirectory, JobNumber, N_bins, Threshold, Sigma, SuppressPlotFlag);
|
||||
spectral_weights(i) = spectral_weight; % Store the spectral weight for the current alpha value
|
||||
g2_all(i, :) = g2; % Store the g2 results for the current alpha value
|
||||
theta_vals_all(i, :) = theta_vals;
|
||||
end
|
||||
|
||||
figure(2);
|
||||
set(gcf,'Position',[100 100 950 750])
|
||||
font = 'Bahnschrift';
|
||||
plot(alpha_vals, spectral_weights, 'o--', ...
|
||||
'LineWidth', 1.5, 'MarkerSize', 6);
|
||||
set(gca, 'FontSize', 14); % For tick labels only
|
||||
hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex');
|
||||
hYLabel = ylabel('Spectral Weight', 'Interpreter', 'tex');
|
||||
hTitle = title('Change across transition', 'Interpreter', 'tex');
|
||||
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
|
||||
|
||||
figure(3);
|
||||
set(gcf,'Position',[100 100 950 750])
|
||||
font = 'Bahnschrift';
|
||||
legend_entries = cell(N_alpha, 1);
|
||||
theta_vals = theta_vals_all(1, :);
|
||||
cmap = sky(N_alpha); % Generate a colormap with enough unique colors
|
||||
hold on
|
||||
|
||||
for i = 1:N_alpha
|
||||
plot(theta_vals/pi, g2_all(i, :), ...
|
||||
'o-', 'Color', cmap(i,:), 'LineWidth', 1.2, ...
|
||||
'MarkerSize', 5);
|
||||
legend_entries{i} = sprintf('$\\alpha = %g^\\circ$', alpha_vals(i));
|
||||
end
|
||||
set(gca, 'FontSize', 14);
|
||||
ylim([-1.5 10.0]); % Set y-axis limits here
|
||||
hXLabel = xlabel('$\delta\theta / \pi$', 'Interpreter', 'latex');
|
||||
hYLabel = ylabel('$g^{(2)}(\delta\theta)$', 'Interpreter', 'latex');
|
||||
hTitle = title('Transition from Droplets to Stripes', 'Interpreter', 'tex');
|
||||
legend(legend_entries, 'Interpreter', 'latex', 'Location', 'bestoutside');
|
||||
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;
|
||||
|
||||
|
||||
@ -46,13 +46,13 @@ function run_on_cluster(batchParams, batchIdx)
|
||||
OptionsStruct.MaxIterationsForGD = 15000;
|
||||
OptionsStruct.TimeStepSize = 1E-3;
|
||||
OptionsStruct.MinimumTimeStepSize = 1E-6;
|
||||
OptionsStruct.TimeCutOff = 5E5;
|
||||
OptionsStruct.TimeCutOff = 2E5;
|
||||
OptionsStruct.EnergyTolerance = 5E-08;
|
||||
OptionsStruct.ResidualTolerance = 1E-05;
|
||||
OptionsStruct.NoiseScaleFactor = 0.010;
|
||||
|
||||
OptionsStruct.PlotLive = false;
|
||||
OptionsStruct.JobNumber = k;
|
||||
OptionsStruct.JobNumber = 0;
|
||||
OptionsStruct.RunOnGPU = true;
|
||||
OptionsStruct.SaveData = true;
|
||||
OptionsStruct.SaveDirectory = saveDir;
|
||||
|
||||
Loading…
Reference in New Issue
Block a user