%% Physical constants
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;
VacuumPermittivity           = 1 / (SpeedOfLight^2 * VacuumPermeability);
HartreeEnergy                = ElectronCharge^2 / (4 * pi * VacuumPermittivity * BohrRadius);
AtomicUnitOfPolarizability   = (ElectronCharge^2 * BohrRadius^2) / HartreeEnergy; % Or simply 4*pi*VacuumPermittivity*BohrRadius^3

% Dy specific constants
Dy164Mass                    = 163.929174751*AtomicMassUnit;
Dy164IsotopicAbundance       = 0.2826;
DyMagneticMoment             = 9.93*BohrMagneton;

%% Bogoliubov excitation spectrum for quasi-2D dipolar gas with QF correction
AtomNumber        = 1E5;                                                                             % Total atom number in the system
wz                = 2 * pi * 72.4;                                                                   % Trap frequency in the tight confinement direction
lz                = sqrt(PlanckConstantReduced/(Dy164Mass * wz));                                    % Defining a harmonic oscillator length
as                = 102.515 * BohrRadius;                                                            % Scattering length
Trapsize          = 7.5815 * lz;                                                                     % Trap is assumed to be a box of finite extent , given here in units of the harmonic oscillator length
alpha             = 0;                                                                               % Polar angle of dipole moment
phi               = 0;                                                                               % Azimuthal angle of momentum vector
MeanWidth         = 5.7304888515 * lz;                                                               % Mean width of Gaussian ansatz
k                 = linspace(0, 2e6, 1000);                                                          % Vector of magnitudes of k vector

%  no = 2.0429e+15, eps_dd = 1.2755, as = 5.4249e-09

AtomNumberDensity = AtomNumber / Trapsize^2;                                                         % Areal density of atoms
add               = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
eps_dd            = add/as;                                                                          % Relative interaction strength
gs                = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as;                                 % Contact interaction strength
gdd               = VacuumPermeability*DyMagneticMoment^2/3;

[Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, alpha, phi);            % DDI potential in k-space

% == Quantum Fluctuations term == %
gammaQF           = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2));
gamma5            = sqrt(2/5) / (sqrt(pi) * MeanWidth)^(3/2);
gQF               = gamma5 * gammaQF;

% == Dispersion relation == %
DeltaK            = ((PlanckConstantReduced^2 .* k.^2) ./ (2 * Dy164Mass)) + ((2 * AtomNumberDensity) .* Ukk) + (3 * gQF * AtomNumberDensity^(3/2));
EpsilonK          = sqrt(((PlanckConstantReduced^2 .* k.^2) ./ (2 * Dy164Mass)) .* DeltaK); 

figure(1)
set(gcf,'Position',[50 50 950 750])
xvals = (k .* add);
yvals = EpsilonK ./ PlanckConstant;
plot(xvals, yvals,LineWidth=2.0)
title(horzcat(['$a_s = ',num2str(round(1/eps_dd,3)),'a_{dd}, '], ['na_{dd}^2 = ',num2str(round(AtomNumberDensity * add^2,4)),'$']),'fontsize',16,'interpreter','latex')
xlabel('$k_{\rho}a_{dd}$','fontsize',16,'interpreter','latex')
ylabel('$\epsilon(k_{\rho})/h$ (Hz)','fontsize',16,'interpreter','latex')
grid on

%% For different interaction strengths

AtomNumber        = 1E5;                                                                             % Total atom number in the system
wz                = 2 * pi * 72.4;                                                                   % Trap frequency in the tight confinement direction
lz                = sqrt(PlanckConstantReduced/(Dy164Mass * wz));                                    % Defining a harmonic oscillator length
Trapsize          = 7.5815 * lz;                                                                     % Trap is assumed to be a box of finite extent , given here in units of the harmonic oscillator length
alpha             = 0;                                                                               % Polar angle of dipole moment
phi               = 0;                                                                               % Azimuthal angle of momentum vector
MeanWidth         = 5.7304888515 * lz;                                                               % Mean width of Gaussian ansatz
k                 = linspace(0, 2e6, 1000);                                                          % Vector of magnitudes of k vector

AtomNumberDensity = AtomNumber / Trapsize^2;                                                         % Areal density of atoms
add               = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length

ScatteringLengths = [108.5, 105.9, 103.3, 102.5150];
eps_dds           = zeros(1, length(ScatteringLengths));
EpsilonKs         = zeros(length(k), length(ScatteringLengths));
for idx = 1:length(ScatteringLengths)

    as                = ScatteringLengths(idx) * BohrRadius;                                             % Scattering length
    eps_dd            = add/as;                                                                          % Relative interaction strength
    gs                = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as;                                 % Contact interaction strength
    gdd               = VacuumPermeability*DyMagneticMoment^2/3;
    
    [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, alpha, phi);            % DDI potential in k-space
    
    % == Quantum Fluctuations term == %
    gammaQF           = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2));
    gamma5            = sqrt(2/5) / (sqrt(pi) * MeanWidth)^(3/2);
    gQF               = gamma5 * gammaQF;
    
    % == Dispersion relation == %
    DeltaK            = ((PlanckConstantReduced^2 .* k.^2) ./ (2 * Dy164Mass)) + ((2 * AtomNumberDensity) .* Ukk) + (3 * gQF * AtomNumberDensity^(3/2));
    EpsilonK          = sqrt(((PlanckConstantReduced^2 .* k.^2) ./ (2 * Dy164Mass)) .* DeltaK); 
    
    eps_dds(idx)      = eps_dd;
    EpsilonKs(:,idx)  = EpsilonK;
end

figure(2)
clf
set(gcf,'Position',[50 50 950 750])
xvals = (k .* add);
yvals = EpsilonKs(:, 1) ./ PlanckConstant;
plot(xvals, yvals,LineWidth=2.0, DisplayName=['$a_s = ',num2str(round(1/eps_dds(1),3)),'a_{dd}$'])
hold on
for idx = 2:length(ScatteringLengths)
    yvals = EpsilonKs(:, idx) ./ PlanckConstant;
    plot(xvals, yvals,LineWidth=2.0, DisplayName=['$a_s = ',num2str(round(1/eps_dds(idx),3)),'a_{dd}$'])
end
title(['$na_{dd}^2 = ',num2str(round(AtomNumberDensity * add^2,4)),'$'],'fontsize',16,'interpreter','latex')
xlabel('$k_{\rho}a_{dd}$','fontsize',16,'interpreter','latex')
ylabel('$\epsilon(k_{\rho})/h$ (Hz)','fontsize',16,'interpreter','latex')
grid on
legend('location', 'northwest','fontsize',16, 'Interpreter','latex')

%% For 3 points on the roton instability boundary

wz                = 2 * pi * 72.4;                                                                   % Trap frequency in the tight confinement direction
lz                = sqrt(PlanckConstantReduced/(Dy164Mass * wz));                                    % Defining a harmonic oscillator length
alpha             = 0;                                                                               % Polar angle of dipole moment
phi               = 0;                                                                               % Azimuthal angle of momentum vector
k                 = linspace(0, 2.25e6, 1000);                                                       % Vector of magnitudes of k vector

nadd2s            = [0.0844, 0.0978, 0.123];
as_to_add         = [0.7730, 0.7840, 0.7819];
var_widths        = [4.97165, 5.7296048721, 5.93178];

add               = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
EpsilonKs         = zeros(length(k), length(nadd2s));

for idx = 1:length(nadd2s)
    
    AtomNumberDensity = nadd2s(idx) / add^2;                                                             % Areal density of atoms
    as                = (as_to_add(idx) * add);                                                          % Scattering length
    eps_dd            = add/as;                                                                          % Relative interaction strength
    gs                = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as;                                 % Contact interaction strength
    gdd               = VacuumPermeability*DyMagneticMoment^2/3;
    MeanWidth         = var_widths(idx) * lz;                                                            % Mean width of Gaussian ansatz

    [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, alpha, phi);            % DDI potential in k-space
    
    % == Quantum Fluctuations term == %
    gammaQF           = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2));
    gamma5            = sqrt(2/5) / (sqrt(pi) * MeanWidth)^(3/2);
    gQF               = gamma5 * gammaQF;
    
    % == Dispersion relation == %
    DeltaK            = ((PlanckConstantReduced^2 .* k.^2) ./ (2 * Dy164Mass)) + ((2 * AtomNumberDensity) .* Ukk) + (3 * gQF * AtomNumberDensity^(3/2));
    EpsilonK          = sqrt(((PlanckConstantReduced^2 .* k.^2) ./ (2 * Dy164Mass)) .* DeltaK); 
    EpsilonKs(:,idx)  = EpsilonK;
end

figure(3)
clf
set(gcf,'Position',[50 50 950 750])
xvals = (k .* add);
yvals = EpsilonKs(:, 1) ./ PlanckConstant;
plot(xvals, yvals,LineWidth=2.0, DisplayName=['$a_s = ',num2str(round(as_to_add(1),4)),'a_{dd}, na_{dd}^2 = ',num2str(round(nadd2s(1),4)),'$'])
hold on
for idx = 2:length(nadd2s)
    yvals = EpsilonKs(:, idx) ./ PlanckConstant;
    plot(xvals, yvals,LineWidth=2.0, DisplayName=['$a_s = ',num2str(round(as_to_add(idx),4)),'a_{dd}, na_{dd}^2 = ',num2str(round(nadd2s(idx),4)),'$'])
end
xlabel('$k_{\rho}a_{dd}$','fontsize',16,'interpreter','latex')
ylabel('$\epsilon(k_{\rho})/h$ (Hz)','fontsize',16,'interpreter','latex')
grid on
legend('location', 'northwest','fontsize',16, 'Interpreter','latex')

%% Mean widths of the variational Gaussian ansatz - extremize the total mean field energy per particle wrt to the variational parameter

wz                       = 2 * pi * 72.4;                                                                   % Trap frequency in the tight confinement direction
lz                       = sqrt(PlanckConstantReduced/(Dy164Mass * wz));                                    % Defining a harmonic oscillator length
gs                       = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as;                                 % Contact interaction strength
add                      = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
gdd                      = VacuumPermeability*DyMagneticMoment^2/3;
AtomNumberDensity        = 0.0978 / add^2;
as                       = 0.784  * add;                                                                    % Scattering length
TotalEnergyPerParticle   = @(x) computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, gs, add, gdd, PlanckConstantReduced);

x0                       = 5;
Aineq                    = []; 
Bineq                    = []; 
Aeq                      = []; 
Beq                      = [];
lb                       = [1]; 
ub                       = [7];
nonlcon                  = [];
fminconopts              = optimoptions(@fmincon,'Display','off', 'StepTolerance', 1.0000e-11, 'MaxIterations',1500);
sigma                    = fmincon(TotalEnergyPerParticle, x0, Aineq, Bineq, Aeq, Beq, lb, ub, nonlcon, fminconopts);
fprintf(['Variational width of Gaussian ansatz = ' num2str(sigma) ' * lz \n'])

%% Mean widths of the variational Gaussian ansatz - extremize the total mean field energy per particle wrt to the variational parameter

wz           = 2 * pi * 72.4;                                                                   % Trap frequency in the tight confinement direction
lz           = sqrt(PlanckConstantReduced/(Dy164Mass * wz));                                    % Defining a harmonic oscillator length
gs           = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as;                                 % Contact interaction strength
add          = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
gdd          = VacuumPermeability*DyMagneticMoment^2/3;

nadd2s       = 0.05:0.001:0.25;
as_to_add    = 0.74:0.001:0.79;
var_widths   = zeros(length(as_to_add), length(nadd2s));

x0           = 5;
Aineq        = []; 
Bineq        = []; 
Aeq          = []; 
Beq          = [];
lb           = [1]; 
ub           = [10];
nonlcon      = [];
fminconopts  = optimoptions(@fmincon,'Display','off', 'StepTolerance', 1.0000e-11, 'MaxIterations',1500);

for idx = 1:length(nadd2s)
    for jdx = 1:length(as_to_add)
        AtomNumberDensity      = nadd2s(idx) / add^2;                                                             % Areal density of atoms
        as                     = (as_to_add(jdx) * add);                                                          % Scattering length
        TotalEnergyPerParticle = @(x) computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, gs, add, gdd, PlanckConstantReduced);
        sigma                  = fmincon(TotalEnergyPerParticle, x0, Aineq, Bineq, Aeq, Beq, lb, ub, nonlcon, fminconopts);
        var_widths(jdx, idx)   = sigma;
    end
end

figure(4)
clf
set(gcf,'Position',[50 50 950 750])
imagesc(nadd2s, as_to_add, var_widths);  % Specify x and y data for axes
set(gca, 'YDir', 'normal');              % Correct the y-axis direction
colorbar;                                % Add a colorbar
xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex');
ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex');

% ====================================================================================================================================================== %

alpha                = 0;                                                                               % Polar angle of dipole moment
phi                  = 0;                                                                               % Azimuthal angle of momentum vector
k                    = linspace(0, 2.25e6, 1000);                                                       % Vector of magnitudes of k vector
instability_boundary = zeros(length(as_to_add), length(nadd2s));

for idx = 1:length(nadd2s)
    for jdx = 1:length(as_to_add)
        AtomNumberDensity = nadd2s(idx) / add^2;                                                             % Areal density of atoms
        as                = (as_to_add(jdx) * add);                                                          % Scattering length
        eps_dd            = add/as;                                                                          % Relative interaction strength
        gs                = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as;                                 % Contact interaction strength
        gdd               = VacuumPermeability*DyMagneticMoment^2/3;
        MeanWidth         = var_widths(jdx, idx) * lz;                                                       % Mean width of Gaussian ansatz
    
        [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, alpha, phi);            % DDI potential in k-space
        
        % == Quantum Fluctuations term == %
        gammaQF           = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2));
        gamma5            = sqrt(2/5) / (sqrt(pi) * MeanWidth)^(3/2);
        gQF               = gamma5 * gammaQF;
        
        % == Dispersion relation == %
        DeltaK            = ((PlanckConstantReduced^2 .* k.^2) ./ (2 * Dy164Mass)) + ((2 * AtomNumberDensity) .* Ukk) + (3 * gQF * AtomNumberDensity^(3/2));
        EpsilonK          = sqrt(((PlanckConstantReduced^2 .* k.^2) ./ (2 * Dy164Mass)) .* DeltaK); 
        instability_boundary(jdx, idx)  = ~isreal(EpsilonK);
    end
end

nadd2s_from_figure     = [0.04974, 0.05383, 0.05655, 0.06609, 0.06916, 0.07291, 0.07836, 0.08517, 0.09063, 0.0978, 0.10459, 0.11345, 0.11822, 0.12231, 0.12674, 0.13117, 0.13560, 0.14003, 0.14548, 0.15127, 0.15775, 0.16660, 0.17546, 0.18364, 0.19557, 0.20579, 0.21839, 0.23850, 0.25144];
as_to_add_from_figure  = [0.76383, 0.76766, 0.76974, 0.77543, 0.77675, 0.77828, 0.78003, 0.78178, 0.78288, 0.7840, 0.78474, 0.78540, 0.78562, 0.78572, 0.78583, 0.78583, 0.78583, 0.78583, 0.78567, 0.78551, 0.78529, 0.78485, 0.78441, 0.78386, 0.78310, 0.78233, 0.78135, 0.77970, 0.77861];

figure(5)
clf
set(gcf,'Position',[50 50 950 750])
imagesc(nadd2s, as_to_add, instability_boundary);  % Specify x and y data for axes
hold on 
plot(nadd2s_from_figure, as_to_add_from_figure, 'r*-', 'LineWidth', 2);  % Plot the curve (red line)
set(gca, 'YDir', 'normal');              % Correct the y-axis direction
colorbar;                                % Add a colorbar
xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex');
ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex');
title('Roton instability boundary','fontsize',16,'interpreter','latex')

%% Roton instability boundary for tilted dipoles

wz           = 2 * pi * 72.4;                                                                   % Trap frequency in the tight confinement direction
lz           = sqrt(PlanckConstantReduced/(Dy164Mass * wz));                                    % Defining a harmonic oscillator length
gs           = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as;                                 % Contact interaction strength
add          = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
gdd          = VacuumPermeability*DyMagneticMoment^2/3;

nadd2s       = 0.05:0.001:0.25;
as_to_add    = 0.76:0.001:0.81;
var_widths   = zeros(length(as_to_add), length(nadd2s));

x0           = 5;
Aineq        = []; 
Bineq        = []; 
Aeq          = []; 
Beq          = [];
lb           = [1]; 
ub           = [10];
nonlcon      = [];
fminconopts  = optimoptions(@fmincon,'Display','off', 'StepTolerance', 1.0000e-11, 'MaxIterations',1500);

for idx = 1:length(nadd2s)
    for jdx = 1:length(as_to_add)
        AtomNumberDensity      = nadd2s(idx) / add^2;                                                             % Areal density of atoms
        as                     = (as_to_add(jdx) * add);                                                          % Scattering length
        TotalEnergyPerParticle = @(x) computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, gs, add, gdd, PlanckConstantReduced);
        sigma                  = fmincon(TotalEnergyPerParticle, x0, Aineq, Bineq, Aeq, Beq, lb, ub, nonlcon, fminconopts);
        var_widths(jdx, idx)   = sigma;
    end
end

% ====================================================================================================================================================== %

alpha                = 45;                                                                               % Polar angle of dipole moment
phi                  = 0;                                                                               % Azimuthal angle of momentum vector
k                    = linspace(0, 2.25e6, 1000);                                                       % Vector of magnitudes of k vector
instability_boundary = zeros(length(as_to_add), length(nadd2s));

for idx = 1:length(nadd2s)
    for jdx = 1:length(as_to_add)
        AtomNumberDensity = nadd2s(idx) / add^2;                                                             % Areal density of atoms
        as                = (as_to_add(jdx) * add);                                                          % Scattering length
        eps_dd            = add/as;                                                                          % Relative interaction strength
        gs                = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as;                                 % Contact interaction strength
        gdd               = VacuumPermeability*DyMagneticMoment^2/3;
        MeanWidth         = var_widths(jdx, idx) * lz;                                                       % Mean width of Gaussian ansatz
    
        [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, alpha, phi);            % DDI potential in k-space
        
        % == Quantum Fluctuations term == %
        gammaQF           = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2));
        gamma5            = sqrt(2/5) / (sqrt(pi) * MeanWidth)^(3/2);
        gQF               = gamma5 * gammaQF;
        
        % == Dispersion relation == %
        DeltaK            = ((PlanckConstantReduced^2 .* k.^2) ./ (2 * Dy164Mass)) + ((2 * AtomNumberDensity) .* Ukk) + (3 * gQF * AtomNumberDensity^(3/2));
        EpsilonK          = sqrt(((PlanckConstantReduced^2 .* k.^2) ./ (2 * Dy164Mass)) .* DeltaK); 
        instability_boundary(jdx, idx)  = ~isreal(EpsilonK);
    end
end

nadd2s_from_figure     = [0.04974, 0.05383, 0.05655, 0.06609, 0.06916, 0.07291, 0.07836, 0.08517, 0.09063, 0.0978, 0.10459, 0.11345, 0.11822, 0.12231, 0.12674, 0.13117, 0.13560, 0.14003, 0.14548, 0.15127, 0.15775, 0.16660, 0.17546, 0.18364, 0.19557, 0.20579, 0.21839, 0.23850, 0.25144];
as_to_add_from_figure  = [0.76383, 0.76766, 0.76974, 0.77543, 0.77675, 0.77828, 0.78003, 0.78178, 0.78288, 0.7840, 0.78474, 0.78540, 0.78562, 0.78572, 0.78583, 0.78583, 0.78583, 0.78583, 0.78567, 0.78551, 0.78529, 0.78485, 0.78441, 0.78386, 0.78310, 0.78233, 0.78135, 0.77970, 0.77861];

figure(6)
clf
set(gcf,'Position',[50 50 950 750])
imagesc(nadd2s, as_to_add, instability_boundary);  % Specify x and y data for axes
hold on 
plot(nadd2s_from_figure, as_to_add_from_figure, 'r*-', 'LineWidth', 2);  % Plot the curve (red line)
set(gca, 'YDir', 'normal');              % Correct the y-axis direction
colorbar;                                % Add a colorbar
xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex');
ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex');
title('Roton instability boundary','fontsize',16,'interpreter','latex')
%%
function [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, alpha, phi)
   Go                        = sqrt(pi) * (k * MeanWidth/sqrt(2)) .* exp((k * MeanWidth/sqrt(2)).^2) .* erfc((k * MeanWidth/sqrt(2)));
   gamma4                    = 1/(sqrt(2*pi) * MeanWidth);
   Fka                       = (3 * cos(deg2rad(alpha))^2 - 1) + ((3 * Go) .* ((sin(deg2rad(alpha))^2 .* sin(deg2rad(phi))^2) - cos(deg2rad(alpha))^2));
   Ukk                       = (gs + (gdd * Fka)) * gamma4;
end

function ret = computeTotalEnergyPerParticle(x, as, AtomNumberDensity, wz, lz, gs, add, gdd, PlanckConstantReduced)
    eps_dd                     = add/as;                                                                     % Relative interaction strength
    MeanWidth                  = x * lz;
    gammaQF                    = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2));                   % Quantum Fluctuations term
    gamma4                     = 1/(sqrt(2*pi) * MeanWidth);
    gamma5                     = sqrt(2/5) / (sqrt(pi) * MeanWidth)^(3/2);
    gQF                        = gamma5 * gammaQF;
    Energy_AxialComponent      = (PlanckConstantReduced * wz) * ((lz^2/(4 * MeanWidth^2)) + (MeanWidth^2/(4 * lz^2)));
    Energy_TransverseComponent = (0.5 * (gs + (2*gdd)) * gamma4 * AtomNumberDensity) + ((2/5) * gQF * AtomNumberDensity^(3/2));
    ret                        = (Energy_AxialComponent + Energy_TransverseComponent) / (PlanckConstantReduced * wz);
end