Latest scripts - added one to compute TF radii of the unmodulated BEC as extracted from the simulations
This commit is contained in:
parent
c90c04249f
commit
99d9413530
170
Dipolar-Gas-Simulator/+Scripts/extractThomasFermiRadii.m
Normal file
170
Dipolar-Gas-Simulator/+Scripts/extractThomasFermiRadii.m
Normal file
@ -0,0 +1,170 @@
|
||||
function [Rx, Ry] = extractThomasFermiRadii(folder_path, run_index, SuppressPlotFlag)
|
||||
set(0,'defaulttextInterpreter','latex')
|
||||
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
|
||||
|
||||
format long
|
||||
|
||||
% Load data
|
||||
Data = load(sprintf(horzcat(folder_path, '/Run_%03i/psi_gs.mat'),run_index),'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);
|
||||
|
||||
% 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));
|
||||
|
||||
state = nxy';
|
||||
|
||||
% Fit
|
||||
fitResult = fitThomasFermiProfile2D(state, x, y);
|
||||
Rx = fitResult.Rx;
|
||||
Ry = fitResult.Ry;
|
||||
|
||||
if ~SuppressPlotFlag
|
||||
fig = figure(1);
|
||||
clf
|
||||
set(gcf,'Position', [100, 100, 1800, 500])
|
||||
t = tiledlayout(1, 3, 'TileSpacing', 'compact', 'Padding', 'compact');
|
||||
|
||||
% Original density plot
|
||||
nexttile;
|
||||
plot1 = pcolor(x, y, state);
|
||||
set(plot1, 'EdgeColor', 'none');
|
||||
colormap(gca, Helper.Colormaps.plasma())
|
||||
colorbar;
|
||||
xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
|
||||
ylabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
|
||||
title('Original Density', 'Interpreter', 'latex', 'FontSize', 16)
|
||||
axis square;
|
||||
set(gca, 'FontSize', 14);
|
||||
|
||||
% Fitted TF profile plot
|
||||
nexttile;
|
||||
plot2 = pcolor(x, y, fitResult.fittedDensity);
|
||||
set(plot2, 'EdgeColor', 'none');
|
||||
colormap(gca, Helper.Colormaps.plasma())
|
||||
colorbar;
|
||||
xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
|
||||
ylabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
|
||||
title(sprintf('TF Fit: $R_x = %.2f\\,\\mu\\mathrm{m}$, $R_y = %.2f\\,\\mu\\mathrm{m}$', Rx, Ry), ...
|
||||
'Interpreter', 'latex', 'FontSize', 16)
|
||||
axis square;
|
||||
set(gca, 'FontSize', 14);
|
||||
|
||||
% 1D cuts through center
|
||||
nexttile;
|
||||
% Find index of center
|
||||
[~, ix0] = min(abs(x - fitResult.x0));
|
||||
[~, iy0] = min(abs(y - fitResult.y0));
|
||||
|
||||
% X-cut: density along y = y0
|
||||
orig_cut_x = state(iy0, :);
|
||||
fit_cut_x = fitResult.fittedDensity(iy0, :);
|
||||
|
||||
% Y-cut: density along x = x0
|
||||
orig_cut_y = state(:, ix0);
|
||||
fit_cut_y = fitResult.fittedDensity(:, ix0);
|
||||
|
||||
% Plot both cuts
|
||||
plot(x, orig_cut_x, 'w-', 'LineWidth', 1.5); hold on
|
||||
plot(x, fit_cut_x, 'r--', 'LineWidth', 1.5)
|
||||
plot(y, orig_cut_y, 'w-', 'LineWidth', 1.5)
|
||||
plot(y, fit_cut_y, 'r--', 'LineWidth', 1.5)
|
||||
hold off
|
||||
axis square;
|
||||
legend({'$n(x)$', '$n_{\mathrm{TF}}(x)$', '$n(y)$', '$n_{\mathrm{TF}}(y)$'}, ...
|
||||
'Interpreter', 'latex', 'FontSize', 12, 'Location', 'northeast')
|
||||
|
||||
xlabel('$x$ or $y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
|
||||
ylabel('Density', 'Interpreter', 'latex', 'FontSize', 14)
|
||||
title('1-D Cut', 'Interpreter', 'latex', 'FontSize', 16)
|
||||
grid on
|
||||
set(gca, 'FontSize', 14)
|
||||
end
|
||||
end
|
||||
|
||||
function fitResult = fitThomasFermiProfile2D(state, xgrid, ygrid)
|
||||
% Fits a 2D Thomas-Fermi profile to the input BEC state
|
||||
%
|
||||
% Input:
|
||||
% state - 2D array (real density or complex wavefunction)
|
||||
% xgrid - vector of x-coordinates (size must match size(state,2))
|
||||
% ygrid - vector of y-coordinates (size must match size(state,1))
|
||||
%
|
||||
% Output (struct):
|
||||
% fitResult.fittedDensity - 2D array of the fitted TF profile
|
||||
% fitResult.n0 - peak density
|
||||
% fitResult.x0, y0 - center coordinates
|
||||
% fitResult.Rx, Ry - Thomas-Fermi radii
|
||||
% fitResult.params - parameter vector [n0, x0, y0, Rx, Ry]
|
||||
|
||||
% Compute density if input is a wavefunction
|
||||
if ~isreal(state)
|
||||
density = abs(state).^2;
|
||||
else
|
||||
density = state;
|
||||
end
|
||||
|
||||
% Set up grid
|
||||
[X, Y] = meshgrid(xgrid, ygrid);
|
||||
xy = [X(:), Y(:)];
|
||||
n = density(:);
|
||||
|
||||
% Mask: ignore near-zero values
|
||||
mask = n > max(n)*1e-3;
|
||||
xy = xy(mask, :);
|
||||
n = n(mask);
|
||||
|
||||
% Initial guess: [n0, x0, y0, Rx, Ry]
|
||||
n0_guess = max(n);
|
||||
x0_guess = mean(xgrid(:));
|
||||
y0_guess = mean(ygrid(:));
|
||||
Rx_guess = (max(xgrid)-min(xgrid)) / 2;
|
||||
Ry_guess = (max(ygrid)-min(ygrid)) / 2;
|
||||
p0 = [n0_guess, x0_guess, y0_guess, Rx_guess, Ry_guess];
|
||||
|
||||
% Bounds
|
||||
lb = [0, min(xgrid), min(ygrid), 0, 0];
|
||||
ub = [Inf, max(xgrid), max(ygrid), Inf, Inf];
|
||||
|
||||
% TF 2D function
|
||||
tf2d = @(p, xy) max(0, p(1) * (1 - ((xy(:,1)-p(2)).^2 / p(4)^2 + (xy(:,2)-p(3)).^2 / p(5)^2)));
|
||||
|
||||
% Fit
|
||||
options = optimoptions('lsqcurvefit','Display','off');
|
||||
pfit = lsqcurvefit(tf2d, p0, xy, n, lb, ub, options);
|
||||
|
||||
% Evaluate fitted profile over full grid
|
||||
fittedDensity = reshape(tf2d(pfit, [X(:), Y(:)]), size(X));
|
||||
|
||||
% Return result
|
||||
fitResult.fittedDensity = fittedDensity;
|
||||
fitResult.n0 = pfit(1);
|
||||
fitResult.x0 = pfit(2);
|
||||
fitResult.y0 = pfit(3);
|
||||
fitResult.Rx = pfit(4);
|
||||
fitResult.Ry = pfit(5);
|
||||
fitResult.params = pfit;
|
||||
end
|
||||
@ -572,7 +572,7 @@ DropletNumber = Scripts.identifyAndCountDroplets(SaveDirectory, JobNum
|
||||
Radius = 2;
|
||||
PeakThreshold = 6E3;
|
||||
JobNumber = 0;
|
||||
TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 0";
|
||||
TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 0^\circ";
|
||||
SuppressPlotFlag = true;
|
||||
|
||||
SCATTERING_LENGTH_RANGE = [85.21 86.25 87.29 88.33 89.38];
|
||||
@ -632,6 +632,144 @@ legend(legendEntries, 'Interpreter', 'tex', 'FontSize', 12, 'Location', 'bestout
|
||||
grid on;
|
||||
hold off;
|
||||
|
||||
%% Plot TF radii of unmodulated states
|
||||
% Parameters
|
||||
JobNumber = 0;
|
||||
TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 40^\circ";
|
||||
SuppressPlotFlag = true;
|
||||
|
||||
SCATTERING_LENGTH_RANGE = [97.71];
|
||||
|
||||
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];
|
||||
|
||||
% 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/Theta40/aS_' ...
|
||||
aS_string '_theta_040_phi_000_N_'];
|
||||
|
||||
% Preallocate results for this curve: rows = N, cols = [Rx, Ry]
|
||||
TF_Radii = zeros(length(NUM_ATOMS_LIST), 2);
|
||||
|
||||
% Loop over all atom numbers
|
||||
for i = 1:length(NUM_ATOMS_LIST)
|
||||
N = NUM_ATOMS_LIST(i);
|
||||
SaveDirectory = [baseDir num2str(N)];
|
||||
|
||||
% Extract both Rx and Ry
|
||||
try
|
||||
[Rx, Ry] = Scripts.extractThomasFermiRadii(SaveDirectory, JobNumber, SuppressPlotFlag);
|
||||
catch ME
|
||||
warning(ME.message)
|
||||
Rx = NaN;
|
||||
Ry = NaN;
|
||||
end
|
||||
|
||||
% Store as row [Rx, Ry]
|
||||
TF_Radii(i, :) = [Rx, Ry];
|
||||
end
|
||||
end
|
||||
|
||||
% Plot curve
|
||||
% Prepare figure
|
||||
fig = figure(1);
|
||||
clf
|
||||
set(gcf,'Position', [100, 100, 1200, 500])
|
||||
t = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact');
|
||||
% Color order for better visibility
|
||||
colors = lines(length(SCATTERING_LENGTH_RANGE));
|
||||
% Original density plot
|
||||
nexttile;
|
||||
plot(NUM_ATOMS_LIST, TF_Radii(:, 1), 'o-', ...
|
||||
'Color', colors(j, :), 'LineWidth', 1.5);
|
||||
% Store legend entry
|
||||
legendEntries{j} = ['a_s = ' num2str(aS) 'a_o'];
|
||||
|
||||
% Finalize plot
|
||||
xlabel('Number of Atoms', 'Interpreter', 'tex', 'FontSize', 16);
|
||||
ylabel('TF Radius - X ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16);
|
||||
legend(legendEntries,'Interpreter', 'tex', 'FontSize', 12, 'Location', 'bestoutside');
|
||||
axis square
|
||||
grid on;
|
||||
|
||||
nexttile
|
||||
plot(NUM_ATOMS_LIST, TF_Radii(:, 2), 'o-', ...
|
||||
'Color', colors(j, :), 'LineWidth', 1.5);
|
||||
% Store legend entry
|
||||
legendEntries{j} = ['a_s = ' num2str(aS) 'a_o'];
|
||||
|
||||
% Finalize plot
|
||||
xlabel('Number of Atoms', 'Interpreter', 'tex', 'FontSize', 16);
|
||||
ylabel('TF Radius - Y ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16);
|
||||
legend(legendEntries,'Interpreter', 'tex', 'FontSize', 12, 'Location', 'bestoutside');
|
||||
axis square
|
||||
grid on;
|
||||
|
||||
sgtitle(TitleString, 'Interpreter', 'tex', 'FontSize', 18);
|
||||
%% Overlay curves
|
||||
|
||||
% File list and associated theta labels
|
||||
fileList = {'TFFermi_Theta0.mat', 'TFFermi_Theta20.mat', 'TFFermi_Theta40.mat'};
|
||||
thetaLabels = {'\theta = 0^\circ', '\theta = 20^\circ', '\theta = 40^\circ'};
|
||||
|
||||
% Set up figure
|
||||
fig = figure(1);
|
||||
clf
|
||||
set(gcf,'Position', [100, 100, 1200, 500])
|
||||
t = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact');
|
||||
|
||||
% Color order for visibility
|
||||
colors = lines(length(fileList));
|
||||
|
||||
legendEntries = cell(1, length(fileList));
|
||||
|
||||
% Loop over files
|
||||
for j = 1:length(fileList)
|
||||
data = load(fileList{j});
|
||||
|
||||
aS = data.SCATTERING_LENGTH_RANGE; % scalar
|
||||
NUM_ATOMS_LIST = data.NUM_ATOMS_LIST;
|
||||
TF_Radii = data.TF_Radii; % [N x 2]
|
||||
|
||||
% Store legend entry (can choose to include theta or a_s)
|
||||
legendEntries{j} = sprintf('%s, a_s = %.2f a_0', thetaLabels{j}, aS);
|
||||
|
||||
% Plot Rx (TF_Radii(:,1))
|
||||
nexttile(1);
|
||||
plot(NUM_ATOMS_LIST, TF_Radii(:, 1), 'o-', ...
|
||||
'Color', colors(j, :), 'LineWidth', 1.5); hold on;
|
||||
|
||||
% Plot Ry (TF_Radii(:,2))
|
||||
nexttile(2);
|
||||
plot(NUM_ATOMS_LIST, TF_Radii(:, 2), 'o-', ...
|
||||
'Color', colors(j, :), 'LineWidth', 1.5); hold on;
|
||||
end
|
||||
|
||||
% Finalize Rx plot
|
||||
nexttile(1);
|
||||
xlabel('Number of Atoms', 'Interpreter', 'tex', 'FontSize', 16);
|
||||
ylabel('TF Radius - X ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16);
|
||||
legend(legendEntries, 'Interpreter', 'tex', 'FontSize', 12, 'Location', 'bestoutside');
|
||||
axis square
|
||||
grid on;
|
||||
|
||||
% Finalize Ry plot
|
||||
nexttile(2);
|
||||
xlabel('Number of Atoms', 'Interpreter', 'tex', 'FontSize', 16);
|
||||
ylabel('TF Radius - Y ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16);
|
||||
legend(legendEntries, 'Interpreter', 'tex', 'FontSize', 12, 'Location', 'bestoutside');
|
||||
axis square
|
||||
grid on;
|
||||
|
||||
% Add super title
|
||||
sgtitle('[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz', ...
|
||||
'Interpreter', 'tex', 'FontSize', 18);
|
||||
|
||||
%% Generate lists
|
||||
|
||||
% Set display format
|
||||
@ -695,6 +833,7 @@ data = {
|
||||
'aS_8.729000e+01_theta_000_phi_000_N_95455'
|
||||
'aS_8.729000e+01_theta_000_phi_000_N_3775000'
|
||||
'aS_9.250000e+01_theta_000_phi_000_N_5000000'
|
||||
'aS_9.771000e+01_theta_020_phi_000_N_3775000'
|
||||
};
|
||||
|
||||
% Preallocate arrays
|
||||
@ -788,12 +927,12 @@ Scripts.plotSmoothedContourOnPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_
|
||||
'interpMethod', 'pchip', 'showPoints', true);
|
||||
|
||||
%% Plot with interpolated phase boundary
|
||||
load('./Results/Data_Full3D/PhaseDiagramTilted_Theta_20.mat')
|
||||
load('./Results/Data_Full3D/PhaseDiagramTilted_Theta_40.mat')
|
||||
PhaseDiagramMatrix = M;
|
||||
ScatteringLengths = SCATTERING_LENGTH_RANGE;
|
||||
NumberOfAtoms = NUM_ATOMS_LIST;
|
||||
PhaseBoundary = load("./Results/Data_Full3D/BoundaryPoints/PhaseBoundary_Tilted_Theta20.mat");
|
||||
TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 20^\circ";
|
||||
PhaseBoundary = load("./Results/Data_Full3D/BoundaryPoints/PhaseBoundary_Tilted_Theta40.mat");
|
||||
TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 40^\circ";
|
||||
|
||||
Scripts.plotPhaseDiagramWithBoundaries(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, PhaseBoundary, TitleString);
|
||||
|
||||
@ -831,5 +970,3 @@ else
|
||||
disp('The state is not modulated');
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user