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

%% 2-D DDI Potential in k-space, with Gaussian ansatz width determined by constrained minimization

wz           = 2 * pi * 300;                                                                    % Trap frequency in the tight confinement direction
lz           = sqrt(PlanckConstantReduced/(Dy164Mass * wz));                                    % Defining a harmonic oscillator length

% Number of grid points in each direction
Params.Nx    = 128;
Params.Ny    = 128;

% Dimensions (in units of l0)
w0           = 2*pi*61.6316;                               % Angular frequency unit [s^-1]
l0           = sqrt(PlanckConstantReduced/(Dy164Mass*w0)); % Defining a harmonic oscillator length
Params.Lx    = 150*l0;
Params.Ly    = 150*l0;
[Transf]     = setupSpace(Params);

nadd2s       = 0.05:0.001:0.25;
as_to_add    = 0.74:0.001:0.79;

add          = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
gdd          = VacuumPermeability*DyMagneticMoment^2/3;

%%
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(10)
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');

%% Chosen values of interaction, density and tilt

nadd              = 0.110;
asadd             = 0.782;
Params.alpha      = 0;
Params.phi        = 0;

[~, nadd2sidx]    = min(abs(nadd2s - nadd));
[~, asaddidx]     = min(abs(as_to_add - asadd));

AtomNumberDensity = nadd2s(nadd2sidx) / add^2;                                                       % Areal density of atoms
as                = (as_to_add(asaddidx) * 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(asaddidx, nadd2sidx)*lz;

% == 2-D DDI Potential in k-space == %
VDk               = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth);
VDk_fftshifted    = fftshift(VDk);

figure(11)
clf
set(gcf,'Position',[50 50 950 750])
imagesc(fftshift(Transf.kx)*l0, fftshift(Transf.ky)*l0, VDk_fftshifted);   % Specify x and y data for axes
set(gca, 'YDir', 'normal');                                     % Correct the y-axis direction
cbar1 = colorbar;
cbar1.Label.Interpreter = 'latex';
xlabel('$k_x l_o$','fontsize',16,'interpreter','latex');
ylabel('$k_y l_o$','fontsize',16,'interpreter','latex');
title(['$\theta = ',num2str(Params.alpha), '; \phi = ', num2str(Params.phi),'$'],'fontsize',16,'interpreter','latex')

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

EpsilonK          = zeros(length(Transf.ky), length(Transf.kx));
gs_tilde          = gs / (sqrt(2*pi) * MeanWidth);

% == Dispersion relation == %
for idx = 1:length(Transf.kx)
    for jdx = 1:length(Transf.ky)
        DeltaK              = ((PlanckConstantReduced^2 .* (Transf.kx(idx).^2 + Transf.ky(jdx).^2)) ./ (2 * Dy164Mass)) + (2 * AtomNumberDensity * gs_tilde) + ((2 * AtomNumberDensity) .* VDk_fftshifted(jdx, idx)) + (3 * gQF * AtomNumberDensity^(3/2));
        EpsilonK(jdx, idx)  = sqrt(((PlanckConstantReduced^2 .* (Transf.kx(idx).^2 + Transf.ky(jdx).^2)) ./ (2 * Dy164Mass)) .* DeltaK); 
    end
end

EpsilonK = double(imag(EpsilonK) ~= 0);  % 'isreal' returns 0 for complex numbers and 1 for real numbers, so we negate it

figure(12)
clf
set(gcf,'Position',[50 50 950 750])
imagesc(fftshift(Transf.kx)*l0, fftshift(Transf.ky)*l0, EpsilonK);   % Specify x and y data for axes
set(gca, 'YDir', 'normal');              % Correct the y-axis direction
cbar1 = colorbar;
cbar1.Label.Interpreter = 'latex';
xlabel('$k_x l_o$','fontsize',16,'interpreter','latex');
ylabel('$k_y l_o$','fontsize',16,'interpreter','latex');
title(['$\theta = ',num2str(Params.alpha), '; \phi = ', num2str(Params.phi),'$'],'fontsize',16,'interpreter','latex')

%%

nadd              = 0.110;
asadd             = 0.782;
% Define values for alpha and phi
alpha_values      = 0:5:90;  % Range of alpha values (you can modify this)
phi_values        = 0:2:90;   % Range of phi values (you can modify this)

[~, nadd2sidx]    = min(abs(nadd2s - nadd));
[~, asaddidx]     = min(abs(as_to_add - asadd));

AtomNumberDensity = nadd2s(nadd2sidx) / add^2;                                                       % Areal density of atoms
as                = (as_to_add(asaddidx) * 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(asaddidx, nadd2sidx)*lz;

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

% Loop over alpha and phi values
for alpha = alpha_values
    for phi = phi_values
        % Update Params with current alpha and phi
        Params.alpha = alpha;
        Params.phi = phi;

        % Compute the potential for the current alpha and phi
        % == 2-D DDI Potential in k-space == %
        VDk               = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth);
        VDk_fftshifted    = fftshift(VDk);
        
        % == 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;
        
        EpsilonK          = zeros(length(Transf.ky), length(Transf.kx));
        gs_tilde          = gs / (sqrt(2*pi) * MeanWidth);
        
        % == Dispersion relation == %
        for idx = 1:length(Transf.kx)
            for jdx = 1:length(Transf.ky)
                DeltaK              = ((PlanckConstantReduced^2 .* (Transf.kx(idx).^2 + Transf.ky(jdx).^2)) ./ (2 * Dy164Mass)) + (2 * AtomNumberDensity * gs_tilde) + ((2 * AtomNumberDensity) .* VDk_fftshifted(jdx, idx)) + (3 * gQF * AtomNumberDensity^(3/2));
                EpsilonK(jdx, idx)  = sqrt(((PlanckConstantReduced^2 .* (Transf.kx(idx).^2 + Transf.ky(jdx).^2)) ./ (2 * Dy164Mass)) .* DeltaK); 
            end
        end

        EpsilonK = double(imag(EpsilonK) ~= 0);  % 'isreal' returns 0 for complex numbers and 1 for real numbers, so we negate it

        % Plot the result
        figure(12)
        clf
        set(gcf,'Position',[50 50 950 750])
        imagesc(fftshift(Transf.kx)*l0, fftshift(Transf.ky)*l0, EpsilonK);   % Specify x and y data for axes
        set(gca, 'YDir', 'normal');              % Correct the y-axis direction
        cbar1 = colorbar;
        cbar1.Label.Interpreter = 'latex';
        xlabel('$k_x$','fontsize',16,'interpreter','latex');
        ylabel('$k_y$','fontsize',16,'interpreter','latex');
        title(['$\theta = ',num2str(Params.alpha), '; \phi = ', num2str(Params.phi),'$'],'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
    end
end

% Close the video file
close(v);

%%
function [Transf] = setupSpace(Params)
    
    Transf.Xmax            = 0.5*Params.Lx;
    Transf.Ymax            = 0.5*Params.Ly;
        
    Nx                     = Params.Nx; 
    Ny                     = Params.Ny;
    
    % Fourier grids
    x                      = linspace(-0.5*Params.Lx,0.5*Params.Lx-Params.Lx/Params.Nx,Params.Nx);
    Kmax                   = pi*Params.Nx/Params.Lx;
    kx                     = linspace(-Kmax,Kmax,Nx+1);
    kx                     = kx(1:end-1); 
    dkx                    = kx(2)-kx(1);
    kx                     = fftshift(kx);
    
    y                      = linspace(-0.5*Params.Ly,0.5*Params.Ly-Params.Ly/Params.Ny,Params.Ny);
    Kmax                   = pi*Params.Ny/Params.Ly;
    ky                     = linspace(-Kmax,Kmax,Ny+1);
    ky                     = ky(1:end-1); 
    dky                    = ky(2)-ky(1);
    ky                     = fftshift(ky);
    
    [Transf.X,Transf.Y]    = ndgrid(x,y);
    [Transf.KX,Transf.KY]  = ndgrid(kx,ky);
    Transf.x               = x; 
    Transf.y               = y;
    Transf.kx              = kx;
    Transf.ky              = ky;
    Transf.dx              = x(2)-x(1); 
    Transf.dy              = y(2)-y(1);
    Transf.dkx             = dkx; 
    Transf.dky             = dky;
end

function VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth)
% == Calculating the DDI potential in Fourier space with appropriate cutoff == %
    
    % Interaction in K space
    QX              = Transf.KX*MeanWidth/sqrt(2);
    QY              = Transf.KY*MeanWidth/sqrt(2);

    Qsq             = QX.^2 + QY.^2;
    absQ            = sqrt(Qsq);
    QDsq            = QX.^2*cos(Params.phi)^2 + QY.^2*sin(Params.phi)^2;
    
    % Bare interaction
    Fpar            = -1 + 3*sqrt(pi)*QDsq.*erfcx(absQ)./absQ; % Scaled complementary error function erfcx(x) = e^(x^2) * erfc(x)
    Fperp           = 2 - 3*sqrt(pi).*absQ.*erfcx(absQ);
    Fpar(absQ == 0) = -1; 
    
    % Full DDI
    VDk             = (Fpar*sin(Params.alpha)^2 + Fperp*cos(Params.alpha)^2);
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