% === Preliminaries ===
clear; clc;

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

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

% === System Parameters ===
AtomNumber         = 1E5;
wz                 = 2 * pi * 72.4;
lz                 = sqrt(PlanckConstantReduced/(Dy164Mass * wz));
Trapsize           = 7.5815 * lz;
theta              = 0;
phi                = 0;
MeanWidth          = 5.7304888515 * lz;
k                  = linspace(0, 3e6, 1000);
AtomNumberDensity  = AtomNumber / Trapsize^2;
add                = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2);
gdd                = VacuumPermeability*DyMagneticMoment^2/3;

% === eps_dd sweep parameters ===
eps_dd_values = linspace(1.15, 1.2755, 40);
max_kval = max(k .* add);
max_yval = 200;  % Hz - manually set this based on known safe upper bound
frames(length(eps_dd_values)) = struct('cdata', [], 'colormap', []);
pauseFrameIdx = 0;

% History for trail
xHist = {};
yHist = {};

% === Plotting ===
figure(1)
set(gcf,'Position',[100 100 950 750])
set(gcf, 'Color', 'w');       % Set figure background to white
set(gca, 'Color', 'w');       % Set axes background to white


for idx = 1:length(eps_dd_values)
    eps_dd = eps_dd_values(idx);
    as = add / eps_dd;
    gs = 4 * pi * PlanckConstantReduced^2 / Dy164Mass * as;

    % Compute DDI potential
    [~, ~, ~, Ukk] = computePotentialInMomentumSpace(k, gs, gdd, MeanWidth, theta, phi);

    % Quantum fluctuations
    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));
    
    % Check for instability
    if any(~isreal(DeltaK))
        disp(['Instability reached at eps_dd = ', num2str(eps_dd)]);
        break
    end

    EpsilonK = sqrt(((PlanckConstantReduced^2 .* k.^2) ./ (2 * Dy164Mass)) .* DeltaK);

    % Store current curve
    xvals = k .* add;
    yvals = EpsilonK ./ PlanckConstant;
    xHist{end+1} = xvals;
    yHist{end+1} = yvals;

    % Plot all past curves
    clf;
    hold on;
    Nfade = length(xHist);
    cmap = coolwarm(Nfade);  % Use custom coolwarm colormap
    for n = 1:5:Nfade-1
        plot(xHist{n}, yHist{n}, 'Color', [cmap(n,:), 0.2 + 0.6*(n/Nfade)], 'LineWidth', 1.5);
    end

    % Highlight current curve
    minDelta = min(DeltaK);
    [~, minIdx] = min(DeltaK);

    if (minDelta < 0.0001 * max(DeltaK)) && (minIdx > 10)
        plot(xvals, yvals, 'Color', [180, 4, 38]/255, 'LineWidth', 5);  % Orange

        % Arrow
        x_arrow = xvals(minIdx);
        y_arrow = yvals(minIdx);
        dx = 0;
        dy = -0.1 * max_yval;

        quiver(x_arrow, y_arrow - dy/2, dx, dy, ...
            'AutoScale', 'off', 'MaxHeadSize', 1.5, ...
            'Color', 'k', 'LineWidth', 2);

        text(x_arrow, y_arrow - dy/2 + 0.05 * max_yval, ...
            '$2\pi/\lambda_{\mathrm{rot}}$', ...
            'Interpreter', 'latex', 'FontSize', 16, 'Color', 'k', ...
            'HorizontalAlignment', 'center');

        pauseFrameIdx = idx;
    else
        plot(xvals, yvals, 'Color', [59, 76, 192]/255, 'LineWidth', 5);  % Blue
    end

    % Axes
    % title(['$a_s = ', num2str(round(1/eps_dd,3)), 'a_{dd}, \quad na_{dd}^2 = ', num2str(round(AtomNumberDensity * add^2,4)), '$'], ...
    %    'fontsize', 16, 'interpreter', 'latex')
    % xlabel('$k_{\rho}a_{dd}$','fontsize',26,'interpreter','latex')
    % ylabel('$\epsilon(k_{\rho})/h$ (Hz)','fontsize',26,'interpreter','latex')
    
    xlim([0, max_kval])
    ylim([0, max_yval])
    set(gca, 'XTick', [], 'YTick', [], 'Color', 'w', 'XColor', 'none', 'YColor', 'none');
    % Arrow positions (in normalized figure coordinates)
    ax = gca;
    ax.Units = 'normalized';
    pos = ax.Position;
    
    % X-axis arrow (horizontal)
    annotation('arrow', ...
        [pos(1), pos(1) + pos(3)*0.95], ...
        [pos(2), pos(2)], ...
        'Color', 'k', 'LineWidth', 2);
    
    % Y-axis arrow (vertical)
    annotation('arrow', ...
        [pos(1), pos(1)], ...
        [pos(2), pos(2) + pos(4)*0.95], ...
        'Color', 'k', 'LineWidth', 2);
    
    % x-label
    text(pos(1) + pos(3)*0.5, pos(2) - 0.15, 'Momentum', ...
        'Units', 'normalized', ...
        'HorizontalAlignment', 'center', ...
        'VerticalAlignment', 'top', ...
        'FontName', 'Bahnschrift' , ...
        'FontSize', 26);
    
    % y-label
    text(pos(1) - 0.2, pos(2) + pos(4)*0.5, 'Energy', ...
        'Units', 'normalized', ...
        'HorizontalAlignment', 'center', ...
        'VerticalAlignment', 'middle', ...
        'FontName', 'Bahnschrift' , ...
        'FontSize', 26, ...
        'Rotation', 90);

    
    drawnow;

    frames(idx) = getframe(gcf);
end

%% === Save as GIF ===
filename = 'roton_dispersion.gif';
delayTime = 0.15;
pauseRepeats = 10;

for i = 1:length(frames)
    [A, map] = rgb2ind(frame2im(frames(i)), 256);
    if i == 1
        imwrite(A, map, filename, 'gif', 'LoopCount', Inf, 'DelayTime', delayTime);
    elseif i == pauseFrameIdx
        for j = 1:pauseRepeats
            imwrite(A, map, filename, 'gif', 'WriteMode', 'append', 'DelayTime', delayTime);
        end
    else
        imwrite(A, map, filename, 'gif', 'WriteMode', 'append', 'DelayTime', delayTime);
    end
end

disp(['GIF saved as ', filename]);

%% === Helper Function ===
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 cmap = coolwarm(n)
    % Generates a coolwarm colormap with n colors
    c1 = [59, 76, 192]/255;
    c2 = [180, 4, 38]/255;
    r = linspace(c1(1), c2(1), n)';
    g = linspace(c1(2), c2(2), n)';
    b = linspace(c1(3), c2(3), n)';
    cmap = [r g b];
end