%% Across range of a_s, n

% load('.\Results\ExtractingParameters_Result_Below1000.mat')
% load('.\Results\ExtractingParameters_Result_Above1000.mat')
load('.\Results\ExtractingParameters_Result_Above10000.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

% Create a tiled layout with tighter spacing
figure(17)
clf
set(gcf,'Position',[50 50 1800 500])
t = tiledlayout(1, 3, '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;
    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

% Third subplot
nexttile;  % Equivalent to subplot(2, 2, 3)
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

% Convert to units relevant to experiment

% Create a tiled layout with tighter spacing
figure(18)
clf
set(gcf,'Position',[50 50 1800 500])
t = tiledlayout(1, 3, '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, (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

% Second subplot
nexttile;  % Equivalent to subplot(2, 2, 2)
for idx = 1:length(data_struct)
    theta_values   = data_struct(idx).theta_values;
    n_values       = data_struct(idx).n_values;
    Lx             = 10e-6;
    Ly             = 10e-6;
    AtomNumber     = n_values .* Lx * Ly;
    plot(theta_values, AtomNumber * 1e-5, '-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('Atom number in a trap of area 100 $\mu m^2 (\times 10^{5})$', 'fontsize', 16, 'interpreter', 'latex');
grid on
legend('location', 'northeast', '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;
    lambda_roton_values = (2 * pi) ./ data_struct(idx).k_roton_values;
    plot(theta_values, lambda_roton_values * 1E6, '-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('$\lambda_{roton} (\mu m)$','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

%% Fixed Density results

load('.\Results\ExtractingParameters_Result_FixedDensity_phi0.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

% Create a tiled layout with tighter spacing
figure(19)
clf
set(gcf,'Position',[50 50 1200 500])
t = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x2 grid

% First subplot
nexttile;
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;
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

% Create a tiled layout with tighter spacing
figure(20)
clf
set(gcf,'Position',[50 50 1200 500])
t = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x2 grid

% First subplot
nexttile;
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', 'northwest', 'fontsize', 10, 'Interpreter', 'latex'); % Reduced font size

% Second subplot
nexttile;
for idx = 1:length(data_struct)
    theta_values        = data_struct(idx).theta_values;
    lambda_roton_values = (2 * pi) ./ data_struct(idx).k_roton_values;
    semilogy(theta_values, lambda_roton_values * 1E6, '-o', LineWidth=2.0, DisplayName=['$w_z = 2 \pi \times $', num2str(data_struct(idx).wz_value), ' Hz']);
    hold on
end
% ylim([0 2])
xlabel('$\theta$','fontsize',16,'interpreter','latex');
ylabel('$\lambda_{roton} (\mu m)$','fontsize',16,'interpreter','latex');
grid on
legend('location', 'southeast','fontsize', 10, 'Interpreter','latex')

% Adjust layout to minimize space
t.TileSpacing = 'compact'; % Minimize space between tiles
t.Padding = 'compact';     % Minimize padding around the layout

%% Fixed Density results - compare two orthogonal directions

data0 = load('.\Results\ExtractingParameters_Result_FixedDensity_phi0.mat');
data90 = load('.\Results\ExtractingParameters_Result_FixedDensity_phi90.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

% Create a tiled layout with tighter spacing
figure(21)
clf
set(gcf,'Position',[50 50 1200 500])
t = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x2 grid

idx = 4;

% First subplot
nexttile;
theta_values   = data0.data_struct(idx).theta_values;
eps_dd_values  = data0.data_struct(idx).eps_dd_values;
plot(theta_values, eps_dd_values, '-o', 'LineWidth', 2.0, 'DisplayName', ['$w_z = 2 \pi \times $', num2str(data0.data_struct(idx).wz_value), ' Hz; $\phi = 0^\circ$']);
hold on
theta_values   = data90.data_struct(idx).theta_values;
eps_dd_values  = data90.data_struct(idx).eps_dd_values;
plot(theta_values, eps_dd_values, '-o', 'LineWidth', 2.0, 'DisplayName', ['$w_z = 2 \pi \times $', num2str(data90.data_struct(idx).wz_value), ' Hz; $\phi = 90^\circ$']);
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;
theta_values   = data0.data_struct(idx).theta_values;
k_roton_values  = data0.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(data0.data_struct(idx).wz_value), ' Hz; $\phi = 0^\circ$']);
hold on
theta_values   = data90.data_struct(idx).theta_values;
k_roton_values  = data90.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(data90.data_struct(idx).wz_value), ' Hz; $\phi = 90^\circ$']);
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


% Create a tiled layout with tighter spacing
figure(22)
clf
set(gcf,'Position',[50 50 1200 500])
t = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x2 grid

% First subplot
nexttile;
theta_values   = data0.data_struct(idx).theta_values;
eps_dd_values  = data0.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(data0.data_struct(idx).wz_value), ' Hz; $\phi = 0^\circ$']);
hold on
theta_values   = data90.data_struct(idx).theta_values;
eps_dd_values  = data90.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(data90.data_struct(idx).wz_value), ' Hz; $\phi = 90^\circ$']);
xlabel('$\theta$', 'fontsize', 16, 'interpreter', 'latex');
ylabel('$a_s (\times a_o)$', 'fontsize', 16, 'interpreter', 'latex');
grid on
legend('location', 'northwest', 'fontsize', 10, 'Interpreter', 'latex'); % Reduced font size

% Second subplot
nexttile;
theta_values   = data0.data_struct(idx).theta_values;
k_roton_values  = data0.data_struct(idx).k_roton_values;
lambda_roton_values = (2 * pi) ./ k_roton_values;
semilogy(theta_values, lambda_roton_values * 1E6, '-o', LineWidth=2.0, DisplayName=['$w_z = 2 \pi \times $', num2str(data0.data_struct(idx).wz_value), ' Hz; $\phi = 0^\circ$']);
hold on
theta_values   = data90.data_struct(idx).theta_values;
k_roton_values  = data90.data_struct(idx).k_roton_values;
lambda_roton_values = (2 * pi) ./ k_roton_values;
semilogy(theta_values, lambda_roton_values * 1E6, '-o', LineWidth=2.0, DisplayName=['$w_z = 2 \pi \times $', num2str(data90.data_struct(idx).wz_value), ' Hz; $\phi = 90^\circ$']);
xlabel('$\theta$','fontsize',16,'interpreter','latex');
ylabel('$\lambda_{roton} (\mu m)$','fontsize',16,'interpreter','latex');
grid on
legend('location', 'northwest','fontsize', 10, 'Interpreter','latex')

% Adjust layout to minimize space
t.TileSpacing = 'compact'; % Minimize space between tiles
t.Padding = 'compact';     % Minimize padding around the layout

%%
%{
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')
%}