PlanckConstant=6.62607015E-34;
PlanckConstantReduced=6.62607015E-34/(2*pi);
FineStructureConstant=7.2973525698E-3;
ElectronMass=9.10938291E-31;
GravitationalConstant=6.67384E-11;
ProtonMass=1.672621777E-27;
AtomicMassUnit=1.660539066E-27; 
BohrRadius=5.2917721067E-11; 
BohrMagneton=9.274009994E-24; 
BoltzmannConstant=1.38064852E-23;
StandardGravityAcceleration=9.80665;
SpeedOfLight=299792458;
StefanBoltzmannConstant=5.670373E-8;
ElectronCharge=1.602176634E-19;
VacuumPermeability=1.25663706212E-6; 
DielectricConstant=8.8541878128E-12;
ElectronGyromagneticFactor=-2.00231930436153;
AvogadroConstant=6.02214076E23;
ZeroKelvin = 273.15;
GravitationalAcceleration = 9.80553;

% Dy specific constants
Dy164Mass                    = 163.929174751*1.660539066E-27;
Dy164IsotopicAbundance       = 0.2826;
DyMagneticMoment             = 9.93*9.274009994E-24;

hbar           = PlanckConstantReduced;           % [J.s]
kbol           = BoltzmannConstant;               % [J/K]
mu0            = VacuumPermeability;              % [N/A^2]
muB            = BohrMagneton;                    % [J/T]
a0             = BohrRadius;                      % [m]
m0             = AtomicMassUnit;                  % [kg]
% w0             = 2*pi*61.658214297935530;       % Angular frequency unit [s^-1]
mu0factor      = 0.3049584233607396;              % =(m0/me)*pi*alpha^2 -- me=mass of electron, alpha=fine struct. const.
m              = Dy164Mass;                       % Mass of Dysprosium
mu             = DyMagneticMoment;                % Dipole lengths (units of muB)
add            = mu0*mu^2*m/(12*pi*hbar^2);
%%

N              = linspace(0.5, 10, 100) * 1E5;
as             = 85 * a0;
w0             = 2*pi*(linspace(50, 1000, 100));         % Angular frequency unit [s^-1]
xs             = sqrt(hbar./(m.*w0));

Q              = zeros(length(xs), length(N));

for i = 1:numel(xs)
    for j = 1:numel(N)
        % Contact interaction strength (units of l0/m)
        C       = 4*pi*as *  N(j)/xs(i);
        
        % DDI strength
        D       = 4*pi*add * N(j)/xs(i);
        
        epsilondd = D/C;
        
        Q(i,j)   = ((4/(3 * pi^2)) * (C^(5/2)/N(j)) * (1 + ((3/2)*epsilondd^2))) * 1E-7;
    end
end

% Plot the resulting matrix
figure(1)
clf
set(gcf,'Position',[50 50 950 750])
set(gca,'FontSize',16,'Box','On','Linewidth',2);
imagesc(w0/(2*pi), N*1E-5, Q);  % Specify x and y data for axes
hold on

% Contours of constant Q
num_levels = 10; % Increase this to get more contours
contour_levels = round(linspace(min(Q(:)*1E-6), max(Q(:)), num_levels),2);
[contour_matrix, contour_handle] = contour(w0/(2*pi), N*1E-5, Q, contour_levels, 'LineColor', 'white', 'LineWidth', 1.5);
clabel(contour_matrix, contour_handle, 'FontSize', 12, 'Color', 'black', 'LabelSpacing', 600);

% Contours of constant N*sqrt(w0)
Nsqrtw0 = N' .* sqrt(w0/(2*pi)) * 1E-7; % Transpose N to match dimensions
% Overlay the contour plot for N*sqrt(w0) in red and dotted
num_levels = 10; % Increase this to get more contours
contour_levels_overlay = round(linspace(min(Nsqrtw0(:)), max(Nsqrtw0(:)), num_levels),2); % Adjust number of contours as needed
[contour_matrix_overlay, contour_handle_overlay] = contour(w0/(2*pi), N*1E-5, Nsqrtw0, contour_levels_overlay, ...
    'LineColor', 'white', 'LineStyle', '--', 'LineWidth', 1.5); % Red and dotted
clabel(contour_matrix_overlay, contour_handle_overlay, 'FontSize', 12, 'Color', 'black', 'LabelSpacing', 500); % Adjust label spacing

% Additional plot formatting
% Add the legend for the two contour types
    legend([contour_handle, contour_handle_overlay], {'Constant Q', 'Constant N \cdot \surd{\omega_0}'},'Interpreter', 'tex', 'Location', 'Best', 'TextColor', 'white', 'FontSize', 16, 'Box', 'off');
set(gca, 'YDir', 'normal');              % Correct the y-axis direction
cbar = colorbar;                                % Add a colorbar
set(gca,'ColorScale','log')
xlabel('Trap Frequency (Hz)','FontSize',16);
ylabel('Atom Numbers (x 1E5)', 'Interpreter', 'tex','FontSize',16);
ylabel(cbar,'Q(x 1E7)', 'Interpreter', 'tex','FontSize',16,'Rotation',270)
title('\bf Scaling of Quantum Fluctuations','fontsize',16, 'Interpreter', 'tex');

%%
N              = linspace(0.5, 10, 100) * 1E5;
  
% Lengths
as             = linspace(100, 70, 100) * a0;
w0             = 2*pi*125;         % Angular frequency unit [s^-1]
xs            = sqrt(hbar/(m*w0));

Q              = zeros(length(as), length(N));

for i = 1:numel(as)
    for j = 1:numel(N)
        % Contact interaction strength (units of l0/m)
        C       = 4*pi*as(i) *  N(j)/xs;
        
        % DDI strength
        D       = 4*pi*add * N(j)/xs;
        
        epsilondd = D/C;
        
        Q(i,j)   = (4/(3 * pi^2)) * (C^(5/2)/N(j)) * (1 + ((3/2)*epsilondd^2));
    end
end

% Plot the resulting matrix
figure(2)
clf
set(gcf,'Position',[50 50 950 750])
set(gca,'FontSize',16,'Box','On','Linewidth',2);
imagesc(as/a0, N*1E-5, Q*1E-5);  % Specify x and y data for axes
hold on
[contour_matrix, contour_handle] = contour(as/a0, N*1E-5, Q*1E-5, 'LineColor', 'white', 'LineWidth', 1.5);
clabel(contour_matrix, contour_handle, 'FontSize', 12, 'Color', 'black', 'LabelSpacing', 600);  % Adjust LabelSpacing for better placement
set(gca, 'YDir', 'normal');              % Correct the y-axis direction
cbar = colorbar;                                % Add a colorbar
xlabel('Scattering Length (a0)','FontSize',16);
ylabel('Atom Numbers (x 1E5)', 'Interpreter', 'tex','FontSize',16);
ylabel(cbar,'Q', 'Interpreter', 'tex','FontSize',16,'Rotation',270)
title('\bf Scaling of Quantum Fluctuations','fontsize',16, 'Interpreter', 'tex');

%%