%% 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;

%% 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
add          = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
gdd          = VacuumPermeability*DyMagneticMoment^2/3;

nadd2s       = 0.005:0.0025:0.5;
as_to_add    = 0.3:0.025:0.95;
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
        gs                     = 4 * pi * PlanckConstantReduced^2/Dy164Mass * as;                                 % Contact interaction strength
        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(7)
clf
set(gcf,'Position',[50 50 1850 750])

theta                = 66;                                                                               % 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));
ScatteringLengths    = zeros(length(as_to_add), 1);
AtomNumber           = zeros(length(nadd2s), 1);
w0                   = 2 * pi * 61.6316;                                                           % Trap frequency in the tight confinement direction
l0                   = sqrt(PlanckConstantReduced/(Dy164Mass * w0));                               % Defining a harmonic oscillator length
tsize                = 10 * l0;

for idx = 1:length(nadd2s)
    for jdx = 1:length(as_to_add)
        AtomNumberDensity      = nadd2s(idx) / add^2;                                                             % Areal density of atoms
        AtomNumber(idx)        = AtomNumberDensity*tsize^2; 
        as                     = (as_to_add(jdx) * add);                                                          % Scattering length
        ScatteringLengths(jdx) = as/BohrRadius;
        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, 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));
        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

subplot(1, 2, 1);  % 1 row, 2 columns, first subplot
imagesc(nadd2s, as_to_add, instability_boundary);  % Specify x and y data for axes
set(gca, 'YDir', 'normal');              % Correct the y-axis direction
colorbar;                                % Add a colorbar
caxis([0 1])
xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex');
ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex');

title(['Along Y: $\theta = ',num2str(theta), '; \phi = ', num2str(phi),'$'],'fontsize',16,'interpreter','latex')

theta                = 66;                                                                               % Polar angle of dipole moment
phi                  = 90;                                                                               % 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));
ScatteringLengths    = zeros(length(as_to_add), 1);
AtomNumber           = zeros(length(nadd2s), 1);
w0                   = 2 * pi * 61.6316;                                                           % Trap frequency in the tight confinement direction
l0                   = sqrt(PlanckConstantReduced/(Dy164Mass * w0));                               % Defining a harmonic oscillator length
tsize                = 10 * l0;

for idx = 1:length(nadd2s)
    for jdx = 1:length(as_to_add)
        AtomNumberDensity      = nadd2s(idx) / add^2;                                                             % Areal density of atoms
        AtomNumber(idx)        = AtomNumberDensity*tsize^2; 
        as                     = (as_to_add(jdx) * add);                                                          % Scattering length
        ScatteringLengths(jdx) = as/BohrRadius;
        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, 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));
        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

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

title(['Along X: $\theta = ',num2str(theta), '; \phi = ', num2str(phi),'$'],'fontsize',16,'interpreter','latex')

%{
imagesc(AtomNumber*1E-5, ScatteringLengths, instability_boundary);   % Specify x and y data for axes
set(gca, 'YDir', 'normal');              % Correct the y-axis direction
cbar1 = colorbar;
cbar1.Label.Interpreter = 'latex';
caxis([0 1])
% ylabel(cbar1,'$(\times 10^{-31})$','FontSize',16,'Rotation',270)
xlabel(' Atom number for a trap area of 100$\mu m^2 ~ (\times 10^5)$','fontsize',16,'interpreter','latex');
ylabel('Scattering length ($\times a_0$)','fontsize',16,'interpreter','latex');
%}

sgtitle('Mean-field instability boundary','fontsize',16,'interpreter','latex')

%% Cycle through angles

% Define values for theta and phi
theta_values      = 0:2:90;  % Range of theta values (you can modify this)

% Set up VideoWriter object to produce a movie
v = VideoWriter('rib_movie', 'MPEG-4'); % Create a video object
v.FrameRate = 5;  % Frame rate of the video
open(v);           % Open the video file

for theta = theta_values
    figure(7)
    clf
    set(gcf,'Position',[50 50 1850 750])
    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));
    ScatteringLengths    = zeros(length(as_to_add), 1);
    AtomNumber           = zeros(length(nadd2s), 1);
    w0                   = 2 * pi * 61.6316;                                                           % Trap frequency in the tight confinement direction
    l0                   = sqrt(PlanckConstantReduced/(Dy164Mass * w0));                               % Defining a harmonic oscillator length
    tsize                = 10 * l0;
    
    for idx = 1:length(nadd2s)
        for jdx = 1:length(as_to_add)
            AtomNumberDensity      = nadd2s(idx) / add^2;                                                             % Areal density of atoms
            AtomNumber(idx)        = AtomNumberDensity*tsize^2; 
            as                     = (as_to_add(jdx) * add);                                                          % Scattering length
            ScatteringLengths(jdx) = as/BohrRadius;
            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, 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));
            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
    
    subplot(1, 2, 1);  % 1 row, 2 columns, first subplot
    imagesc(nadd2s, as_to_add, instability_boundary);  % Specify x and y data for axes
    set(gca, 'YDir', 'normal');              % Correct the y-axis direction
    colorbar;                                % Add a colorbar
    caxis([0 1])
    xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex');
    ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex');
    
    title(['Along Y: $\theta = ',num2str(theta), '; \phi = ', num2str(phi),'$'],'fontsize',16,'interpreter','latex')
    
    phi                  = 90;                                                                               % 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));
    ScatteringLengths    = zeros(length(as_to_add), 1);
    AtomNumber           = zeros(length(nadd2s), 1);
    w0                   = 2 * pi * 61.6316;                                                           % Trap frequency in the tight confinement direction
    l0                   = sqrt(PlanckConstantReduced/(Dy164Mass * w0));                               % Defining a harmonic oscillator length
    tsize                = 10 * l0;
    
    for idx = 1:length(nadd2s)
        for jdx = 1:length(as_to_add)
            AtomNumberDensity      = nadd2s(idx) / add^2;                                                             % Areal density of atoms
            AtomNumber(idx)        = AtomNumberDensity*tsize^2; 
            as                     = (as_to_add(jdx) * add);                                                          % Scattering length
            ScatteringLengths(jdx) = as/BohrRadius;
            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, 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));
            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
    
    % set(gcf,'Position',[50 50 950 750])
    subplot(1, 2, 2);  % 1 row, 2 columns, first subplot
    imagesc(nadd2s, as_to_add, instability_boundary);  % Specify x and y data for axes
    set(gca, 'YDir', 'normal');              % Correct the y-axis direction
    colorbar;                                % Add a colorbar
    caxis([0 1])
    xlabel('$na_{dd}^2$','fontsize',16,'interpreter','latex');
    ylabel('$a_s/a_{dd}$','fontsize',16,'interpreter','latex');
    
    title(['Along X: $\theta = ',num2str(theta), '; \phi = ', num2str(phi),'$'],'fontsize',16,'interpreter','latex')
    
    %{
    imagesc(AtomNumber*1E-5, ScatteringLengths, instability_boundary);   % Specify x and y data for axes
    set(gca, 'YDir', 'normal');              % Correct the y-axis direction
    cbar1 = colorbar;
    cbar1.Label.Interpreter = 'latex';
    caxis([0 1])
    % ylabel(cbar1,'$(\times 10^{-31})$','FontSize',16,'Rotation',270)
    xlabel(' Atom number for a trap area of 100$\mu m^2 ~ (\times 10^5)$','fontsize',16,'interpreter','latex');
    ylabel('Scattering length ($\times a_0$)','fontsize',16,'interpreter','latex');
    %}
    
    % Capture the frame and write to video
    frame = getframe(gcf); % Capture the current figure
    writeVideo(v, frame);  % Write the frame to the video

    % sgtitle('Mean-field instability boundary','fontsize',16,'interpreter','latex')
end

% Close the video file
close(v);

%%
function [Go,gamma4,Fka,Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, theta, 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(theta))^2 - 1) + ((3 * Go) .* ((sin(deg2rad(theta))^2 .* sin(deg2rad(phi))^2) - cos(deg2rad(theta))^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