diff --git a/Quasi2DBogoliubovSpectrum.m b/Quasi2DBogoliubovSpectrum.m
index 2d96620..e21d5a3 100644
--- a/Quasi2DBogoliubovSpectrum.m
+++ b/Quasi2DBogoliubovSpectrum.m
@@ -28,48 +28,49 @@ Dy164Mass                    = 163.929174751*1.660539066E-27;
 Dy164IsotopicAbundance       = 0.2826;
 DyMagneticMoment             = 9.93*9.274009994E-24;
 
-%% Dispersion relation of the quasiparticle excitations
-AtomNumber        = 1E5;
-wz                = 2*pi*72.4;
+%% 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.4*BohrRadius;                                                                % Scattering length
-Trapsize          = 7.6;
-alpha             = 0;
-phi               = 0;
-MeanWidth         = 2.8215042184E3*lz;
-k                 = linspace(0, 1e7, 1000);
+as                = 102.515*BohrRadius;                                                              % Scattering length
+Trapsize          = 7.5815;                                                                          % 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, 3e6, 1000);                                                          % Vector of magnitudes of k vector
 
-
-AtomNumberDensity = AtomNumber / (Trapsize * lz)^2;
+AtomNumberDensity = AtomNumber / (Trapsize * lz)^2;                                                  % Areal density of atoms
 add               = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
-eps_dd            = add/as;
+eps_dd            = add/as;                                                                          % Relative interaction strength
 gs                = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as;                                 % Contact interaction strength
+gdd               = VacuumPermeability*DyMagneticMoment^2/3;
 
-[fk,Fka,Ukk] = computePotentialInMomentumSpace(k, lz, alpha, phi, gs, eps_dd);
+[Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, alpha, phi);            % DDI potential in k-space
 
 % == Quantum Fluctuations term == %
-gQF               = ((256 * PlanckConstantReduced^2) / (15*Dy164Mass*MeanWidth^3)) * as^(5/2) * (1 + ((3/2) * eps_dd^2));
-
+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',[100 100 950 750])
-% xvals = (k .* lz/sqrt(2));
+set(gcf,'Position',[50 50 950 750])
 xvals = (k .* add);
-yvals = EpsilonK ./ (PlanckConstantReduced * wz);
+yvals = EpsilonK ./ PlanckConstant;
 plot(xvals, yvals,LineWidth=2.0)
-% xlim([3.45, 3.65])
-% ylim([0, 0.001])
-title(horzcat(['$a_s = ',num2str(1/eps_dd),'a_{dd}, '], ['na_{dd}^2 = ',num2str(AtomNumberDensity * add^2),'$']),'fontsize',16,'interpreter','latex')
-xlabel('$ka_{dd}$','fontsize',16,'interpreter','latex')
-ylabel('$\epsilon(k)/\hbar \omega_z$','fontsize',16,'interpreter','latex')
+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
 
 %%
-function [fk,Fka,Ukk] = computePotentialInMomentumSpace(k, lz, alpha, phi, gs, eps_dd)
-    fk                = (3 * sqrt(pi)) * (k .* lz/sqrt(2)) .* exp((k .* lz/sqrt(2)).^2) .* erfc((k .* lz/sqrt(2))) ;
-    Fka               = (fk .* sin(deg2rad(phi))^2 - 1) + (cos(deg2rad(alpha))^2 .* (3 - (fk .* (sin(deg2rad(phi))^2 + 1))));
-    Ukk               = (gs/ (sqrt(2 * pi) * lz)) .* (1 + (eps_dd .* Fka));
+
+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
\ No newline at end of file