Latest - Plotting routines for results added.

Estimations/DipolarDispersionAndRotonInstabilityBoundary/AnalyzeResults.m

@@ -0,0 +1,135 @@
+load('ExtractingKRoton_Result.mat')
+
+PlanckConstantReduced = 6.62607015E-34/(2*pi);
+AtomicMassUnit        = 1.660539066E-27; 
+Dy164Mass             = 163.929174751*AtomicMassUnit;
+VacuumPermeability    = 1.25663706212E-6; 
+BohrMagneton          = 9.274009994E-24; 
+BohrRadius            = 5.2917721067E-11; 
+DyMagneticMoment      = 9.93*BohrMagneton;
+
+add          = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
+
+%{
+figure(13)
+clf
+set(gcf,'Position',[50 50 950 750])
+for idx = 1:length(data_struct)
+    theta_values   = data_struct(idx).theta_values;
+    eps_dd_values  = data_struct(idx).eps_dd_values;
+    plot(theta_values, eps_dd_values, '-o', LineWidth=2.0, DisplayName=['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
+    hold on
+end
+xlabel('$\theta$','fontsize',16,'interpreter','latex');
+ylabel('$\epsilon_{dd}$','fontsize',16,'interpreter','latex');
+% title([''],'fontsize',16,'interpreter','latex')
+grid on
+legend('location', 'northeast','fontsize', 16, 'Interpreter','latex')
+
+figure(14)
+clf
+set(gcf,'Position',[50 50 950 750])
+for idx = 1:length(data_struct)
+    theta_values   = data_struct(idx).theta_values;
+    eps_dd_values  = data_struct(idx).eps_dd_values;
+    plot(theta_values, (1./eps_dd_values) * (add/BohrRadius), '-o', LineWidth=2.0, DisplayName=['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
+    hold on
+end
+xlabel('$\theta$','fontsize',16,'interpreter','latex');
+ylabel('$a_s (\times a_o)$','fontsize',16,'interpreter','latex');
+% title([''],'fontsize',16,'interpreter','latex')
+grid on
+legend('location', 'southeast','fontsize', 16, 'Interpreter','latex')
+
+figure(15)
+clf
+set(gcf,'Position',[50 50 950 750])
+for idx = 1:length(data_struct)
+    theta_values   = data_struct(idx).theta_values;
+    n_values       = data_struct(idx).n_values;
+    plot(theta_values, n_values * 1E-15, '-o', LineWidth=2.0, DisplayName=['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
+    hold on
+end
+xlabel('$\theta$','fontsize',16,'interpreter','latex');
+ylabel('$n (\times 10^{3} \mu m^{-2})$','fontsize',16,'interpreter','latex');
+% title([''],'fontsize',16,'interpreter','latex')
+grid on
+legend('location', 'northeast','fontsize', 16, 'Interpreter','latex')
+
+figure(16)
+clf
+set(gcf,'Position',[50 50 950 750])
+for idx = 1:length(data_struct)
+    theta_values   = data_struct(idx).theta_values;
+    k_roton_values = data_struct(idx).k_roton_values;
+    plot(theta_values, k_roton_values * 1E-6, '-o', LineWidth=2.0, DisplayName=['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
+    hold on
+end
+xlabel('$\theta$','fontsize',16,'interpreter','latex');
+ylabel('$k_{roton} (\mu m^{-1})$','fontsize',16,'interpreter','latex');
+% title([''],'fontsize',16,'interpreter','latex')
+grid on
+legend('location', 'northeast','fontsize', 16, 'Interpreter','latex')
+%}
+
+% Create a tiled layout with tighter spacing
+figure(17)
+clf
+set(gcf,'Position',[50 50 1200 900])
+t = tiledlayout(2, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x2 grid
+
+% First subplot
+nexttile;  % Equivalent to subplot(2, 2, 1)
+for idx = 1:length(data_struct)
+    theta_values   = data_struct(idx).theta_values;
+    eps_dd_values  = data_struct(idx).eps_dd_values;
+    plot(theta_values, eps_dd_values, '-o', 'LineWidth', 2.0, 'DisplayName', ['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
+    hold on
+end
+xlabel('$\theta$', 'fontsize', 16, 'interpreter', 'latex');
+ylabel('$\epsilon_{dd}$', 'fontsize', 16, 'interpreter', 'latex');
+grid on
+legend('location', 'northeast', 'fontsize', 10, 'Interpreter', 'latex'); % Reduced font size
+
+% Second subplot
+nexttile;  % Equivalent to subplot(2, 2, 2)
+for idx = 1:length(data_struct)
+    theta_values   = data_struct(idx).theta_values;
+    eps_dd_values  = data_struct(idx).eps_dd_values;
+    plot(theta_values, (1 ./ eps_dd_values) * (add / BohrRadius), '-o', 'LineWidth', 2.0, 'DisplayName', ['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
+    hold on
+end
+xlabel('$\theta$', 'fontsize', 16, 'interpreter', 'latex');
+ylabel('$a_s (\times a_o)$', 'fontsize', 16, 'interpreter', 'latex');
+grid on
+legend('location', 'southeast', 'fontsize', 10, 'Interpreter', 'latex'); % Reduced font size
+
+% Third subplot
+nexttile;  % Equivalent to subplot(2, 2, 3)
+for idx = 1:length(data_struct)
+    theta_values   = data_struct(idx).theta_values;
+    n_values       = data_struct(idx).n_values;
+    plot(theta_values, n_values * 1E-15, '-o', 'LineWidth', 2.0, 'DisplayName', ['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
+    hold on
+end
+xlabel('$\theta$', 'fontsize', 16, 'interpreter', 'latex');
+ylabel('$n (\times 10^{3} \mu m^{-2})$', 'fontsize', 16, 'interpreter', 'latex');
+grid on
+legend('location', 'northeast', 'fontsize', 10, 'Interpreter', 'latex'); % Reduced font size
+
+% Fourth subplot
+nexttile;  % Equivalent to subplot(2, 2, 4)
+for idx = 1:length(data_struct)
+    theta_values   = data_struct(idx).theta_values;
+    k_roton_values = data_struct(idx).k_roton_values;
+    plot(theta_values, k_roton_values * 1E-6, '-o', LineWidth=2.0, DisplayName=['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
+    hold on
+end
+xlabel('$\theta$','fontsize',16,'interpreter','latex');
+ylabel('$k_{roton} (\mu m^{-1})$','fontsize',16,'interpreter','latex');
+grid on
+legend('location', 'northeast','fontsize', 10, 'Interpreter','latex')
+
+% Adjust layout to minimize space
+t.TileSpacing = 'compact'; % Minimize space between tiles
+t.Padding = 'compact';     % Minimize padding around the layout

Estimations/DipolarDispersionAndRotonInstabilityBoundary/ExtractingKRoton.m

@@ -32,17 +32,17 @@ DyMagneticMoment             = 9.93*BohrMagneton;
 
 %-------TEST-------%
 % nadd2s       = 0.05:0.005:0.25;
-% as_to_add    = 0.76:0.0001:0.81;
+% as_to_add    = 0.76:0.001:0.81;
 
 %-------DEPLOY-------%
 nadd2s       = 0.005:0.005:0.5;
-as_to_add    = 0.35:0.0001:1.15;
+as_to_add    = 0.35:0.001:1.15;
 
 data_struct  = struct;
 wz_values    = [150, 300, 500, 750];
-theta_values = 0:5:45;  % Range of theta values (you can modify this)
+theta_values = 0:5:45;                                                                       % Range of theta values
-phi          = 0;                                                                               % Azimuthal angle of momentum vector
+phi          = 0;                                                                            % Azimuthal angle of momentum vector
-kvec         = linspace(0, 2.25e6, 1000);                                                       % Vector of magnitudes of k vector
+kvec         = linspace(0, 5e6, 1000);                                                       % Vector of magnitudes of k vector
 
 for mainloop_idx = 1:length(wz_values)
     format long
@@ -280,6 +280,7 @@ function [eps_dd, AtomNumberDensity, k_roton] = extractFromBoundaryCurve(theta,
         val                = rounded_val;
     else
         val                = nearest_val;
+        idx                = nearest_idx;
     end
 
     AtomNumberDensity      = xvals(idx) / add^2;                                                              % Areal density of atoms
@@ -298,7 +299,7 @@ function [eps_dd, AtomNumberDensity, k_roton] = extractFromBoundaryCurve(theta,
     TotalEnergyPerParticle = @(x) computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, gs, add, gdd, PlanckConstantReduced);
     sigma                  = fmincon(TotalEnergyPerParticle, x0, Aineq, Bineq, Aeq, Beq, lb, ub, nonlcon, fminconopts);
     MeanWidth              = sigma * lz;                                                                      % Mean width of Gaussian ansatz
-    [Go,gamma4,Fka,Ukk]    = computePotentialInMomentumSpace(kvec, gs, gdd, MeanWidth, theta, phi);              % DDI potential in k-space
+    [Go,gamma4,Fka,Ukk]    = computePotentialInMomentumSpace(kvec, gs, gdd, MeanWidth, theta, phi);           % DDI potential in k-space
     
     % == Quantum Fluctuations term == %
     gammaQF                = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2));