karthik/Calculations
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

Corrected script for calculating the dispersion relation in 2D - thanks to Tom Bland!

Browse Source
This commit is contained in:
Karthik 2025-05-05 18:49:45 +02:00
parent 2d6fb3f15c
commit 4583f1a643
2 changed files with 58 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.gitignore vendored
View File

@ -10,5 +10,7 @@ Time-Series-Analyzer/Time-Series-Data
*.mp4 *.mp4
*.bat *.bat
*.json *.json
*.txt
*.asv
.ipynb_checkpoints/ .ipynb_checkpoints/
.vscode/ .vscode/

84
Estimations/RotonInstability/DipolarDispersion2D.m
View File

@ -31,7 +31,7 @@ DyMagneticMoment = 9.93*BohrMagneton;
%% 2-D DDI Potential in k-space, with Gaussian ansatz width determined by constrained minimization %% 2-D DDI Potential in k-space, with Gaussian ansatz width determined by constrained minimization
% With user-defined values of interaction, density and tilt % With user-defined values of interaction, density and tilt
% w0 = 2*pi*61.6316; % Angular frequency unit [s^-1] % w0 = 2*pi*61.6316; % Angular frequency unit [s^-1]
% l0 = sqrt(PlanckConstantReduced/(Dy164Mass*w0)); % l0 = sqrt(PlanckConstantReduced/(Dy164Mass*w0));
% % Defining a harmonic oscillator length - heredue to the choice of w0, l0 % % Defining a harmonic oscillator length - heredue to the choice of w0, l0
% is 1 micrometer % is 1 micrometer
@ -40,8 +40,8 @@ wz = 2 * pi * 300;
lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length
% Number of grid points in each direction % Number of grid points in each direction
Params.Nx = 128; Params.Nx = 256;
Params.Ny = 128; Params.Ny = 256;
Params.Lx = 150*1e-6; Params.Lx = 150*1e-6;
Params.Ly = 150*1e-6; Params.Ly = 150*1e-6;
[Transf] = setupSpace(Params); [Transf] = setupSpace(Params);
@ -49,7 +49,7 @@ Params.Ly = 150*1e-6;
nadd2s = 0.110; % Number density * add^2 nadd2s = 0.110; % Number density * add^2
as_to_add = 0.782; % 1/edd as_to_add = 0.782; % 1/edd
Params.theta = 57; % Polar angle of dipole moment Params.theta = 57; % Polar angle of dipole moment
Params.eta = 0; % Azimuthal angle of dipole moment Params.phi = 0; % Azimuthal angle of dipole moment
add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length add = VacuumPermeability*DyMagneticMoment^2*Dy164Mass/(12*pi*PlanckConstantReduced^2); % Dipole length
gdd = VacuumPermeability*DyMagneticMoment^2/3; gdd = VacuumPermeability*DyMagneticMoment^2/3;
@ -74,19 +74,19 @@ sigma = fmincon(TotalEnergyPerParticle, x0, Aineq, Bineq, Aeq,
MeanWidth = sigma * lz; MeanWidth = sigma * lz;
% == 2-D DDI Potential in k-space == % % == 2-D DDI Potential in k-space == %
VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth); VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth); % Transf.K is in shifted space, so VDk is shifted
VDk_fftshifted = fftshift(VDk); VDk_fftshifted = fftshift(VDk); % This is in normal space, but dispersion is calculated in shifted space, so unshifted VDk needs to be used there
figure(8) figure(11)
clf clf
set(gcf,'Position',[50 50 950 750]) set(gcf,'Position',[50 50 950 750])
imagesc(fftshift(Transf.kx)*1e-6, fftshift(Transf.ky)*1e-6, VDk_fftshifted); % Specify x and y data for axes imagesc(fftshift(Transf.kx)*1e-6, fftshift(Transf.ky)*1e-6, VDk_fftshifted); % Specify x and y data for axes
set(gca, 'YDir', 'normal'); % Correct the y-axis direction set(gca, 'YDir', 'normal'); % Correct the y-axis direction
cbar1 = colorbar; cbar1 = colorbar;
cbar1.Label.Interpreter = 'latex'; cbar1.Label.Interpreter = 'latex';
xlabel('$k_x l_o$','fontsize',16,'interpreter','latex'); xlabel('$k_x l_o$','fontsize',16,'interpreter','latex');
ylabel('$k_y l_o$','fontsize',16,'interpreter','latex'); ylabel('$k_y l_o$','fontsize',16,'interpreter','latex');
title(['2-D DDI Potential: $\theta = ',num2str(Params.theta), '; \eta = ', num2str(Params.eta),'$'],'fontsize',16,'interpreter','latex') title(['2-D DDI Potential: $\theta = ',num2str(Params.alpha), '; \phi = ', num2str(Params.phi),'$'],'fontsize',16,'interpreter','latex')
% == Quantum Fluctuations term == % % == Quantum Fluctuations term == %
gammaQF = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2)); gammaQF = (32/3) * gs * (as^3/pi)^(1/2) * (1 + ((3/2) * eps_dd^2));
@ -95,47 +95,74 @@ gQF = gamma5 * gammaQF;
EpsilonK = zeros(length(Transf.ky), length(Transf.kx)); EpsilonK = zeros(length(Transf.ky), length(Transf.kx));
gs_tilde = gs / (sqrt(2*pi) * MeanWidth); gs_tilde = gs / (sqrt(2*pi) * MeanWidth);
gdd_tilde = gdd / (sqrt(2*pi) * MeanWidth); % you were missing this!
% == Dispersion relation == % % == Dispersion relation == %
for idx = 1:length(Transf.kx) for idx = 1:length(Transf.kx)
for jdx = 1:length(Transf.ky) 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)); DeltaK = ((PlanckConstantReduced^2 .* (Transf.kx(idx).^2 + Transf.ky(jdx).^2)) ./ (2 * Dy164Mass)) + (2 * AtomNumberDensity * gs_tilde) + ((2 * AtomNumberDensity) * gdd_tilde .* VDk(jdx, idx)) + (3 * gQF * AtomNumberDensity^(3/2));
EpsilonK(jdx, idx) = sqrt(((PlanckConstantReduced^2 .* (Transf.kx(idx).^2 + Transf.ky(jdx).^2)) ./ (2 * Dy164Mass)) .* DeltaK); EpsilonK(jdx, idx) = sqrt(((PlanckConstantReduced^2 .* (Transf.kx(idx).^2 + Transf.ky(jdx).^2)) ./ (2 * Dy164Mass)) .* DeltaK);
end end
end end
EpsilonK = double(imag(EpsilonK) ~= 0); % 'isreal' returns 0 for complex numbers and 1 for real numbers, so we negate it EpsilonK = real(EpsilonK);
figure(9) figure(12)
clf clf
set(gcf,'Position',[50 50 950 750]) set(gcf,'Position',[50 50 950 750])
imagesc(fftshift(Transf.kx)*1e-6, fftshift(Transf.ky)*1e-6, EpsilonK); % Specify x and y data for axes imagesc(fftshift(Transf.kx)*1e-6, fftshift(Transf.ky)*1e-6, fftshift(EpsilonK)); % Specify x and y data for axes
set(gca, 'YDir', 'normal'); % Correct the y-axis direction set(gca, 'YDir', 'normal'); % Correct the y-axis direction
cbar1 = colorbar; cbar1 = colorbar;
cbar1.Label.Interpreter = 'latex'; cbar1.Label.Interpreter = 'latex';
xlabel('$k_x l_o$','fontsize',16,'interpreter','latex'); xlabel('$k_x l_o$','fontsize',16,'interpreter','latex');
ylabel('$k_y l_o$','fontsize',16,'interpreter','latex'); ylabel('$k_y l_o$','fontsize',16,'interpreter','latex');
title(['2-D Dispersion: $\theta = ',num2str(Params.theta), '; \eta = ', num2str(Params.eta),'$'],'fontsize',16,'interpreter','latex') title(['2-D Dispersion: $\theta = ',num2str(Params.alpha), '; \phi = ', num2str(Params.phi),'$'],'fontsize',16,'interpreter','latex')
%% 3-D Plot
% Assuming Transf.kx, Transf.ky, and EpsilonK are already defined
% Apply fftshift and scaling
kx_shifted = fftshift(Transf.kx) * 1e-6; % Convert to desired units
ky_shifted = fftshift(Transf.ky) * 1e-6; % Convert to desired units
EpsilonK_shifted = fftshift(EpsilonK); % This is the energy or function you want to plot
% Create meshgrid for the kx and ky values
[Kx, Ky] = meshgrid(kx_shifted, ky_shifted);
% Create the 3D plot
figure(13)
clf
set(gcf,'Position',[50 50 950 750])
surf(Kx, Ky, EpsilonK_shifted);
shading interp; % Interpolates shading between grid points
% Label axes
xlabel('$k_x l_o$','fontsize',16,'interpreter','latex');
ylabel('$k_y l_o$','fontsize',16,'interpreter','latex');
zlabel('$\epsilon(k)$','fontsize',16,'interpreter','latex');
title(['2-D Dispersion: $\theta = ',num2str(Params.alpha), '; \phi = ', num2str(Params.phi),'$'],'fontsize',16,'interpreter','latex')
colorbar; % Optional, to show a color scale for the surface plot
% Apply a different colormap
colormap('jet'); % 'jet', 'parula', 'hot', etc.
colorbar; % Add colorbar to see the range of values
%% Cycle through angles %% Cycle through angles
% Define values for theta and eta % Define values for theta and phi
theta_values = 0:10:90; % Range of theta values (you can modify this) theta_values = 0:10:90; % Range of theta values (you can modify this)
eta_values = 0:10:90; % Range of eta values (you can modify this) phi_values = 0:10:90; % Range of phi values (you can modify this)
% Set up VideoWriter object to produce a movie % Set up VideoWriter object to produce a movie
% v = VideoWriter('potential_movie', 'MPEG-4'); % Create a video object % v = VideoWriter('potential_movie', 'MPEG-4'); % Create a video object
% v.FrameRate = 5; % Frame rate of the video % v.FrameRate = 5; % Frame rate of the video
% open(v); % Open the video file % open(v); % Open the video file
% Loop over theta and eta values % Loop over theta and phi values
for theta = theta_values for theta = theta_values
for eta = eta_values for phi = phi_values
% Update Params with current theta and eta % Update Params with current theta and phi
Params.theta = theta; Params.theta = theta;
Params.eta = eta; Params.phi = phi;
% Compute the potential for the current theta and eta % Compute the potential for the current theta and phi
% == 2-D DDI Potential in k-space == % % == 2-D DDI Potential in k-space == %
VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth); VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth);
VDk_fftshifted = fftshift(VDk); VDk_fftshifted = fftshift(VDk);
@ -147,28 +174,29 @@ for theta = theta_values
EpsilonK = zeros(length(Transf.ky), length(Transf.kx)); EpsilonK = zeros(length(Transf.ky), length(Transf.kx));
gs_tilde = gs / (sqrt(2*pi) * MeanWidth); gs_tilde = gs / (sqrt(2*pi) * MeanWidth);
gdd_tilde = gdd / (sqrt(2*pi) * MeanWidth); % you were missing this!
% == Dispersion relation == % % == Dispersion relation == %
for idx = 1:length(Transf.kx) for idx = 1:length(Transf.kx)
for jdx = 1:length(Transf.ky) 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)); DeltaK = ((PlanckConstantReduced^2 .* (Transf.kx(idx).^2 + Transf.ky(jdx).^2)) ./ (2 * Dy164Mass)) + (2 * AtomNumberDensity * gs_tilde) + ((2 * AtomNumberDensity) * gdd_tilde .* VDk(jdx, idx)) + (3 * gQF * AtomNumberDensity^(3/2));
EpsilonK(jdx, idx) = sqrt(((PlanckConstantReduced^2 .* (Transf.kx(idx).^2 + Transf.ky(jdx).^2)) ./ (2 * Dy164Mass)) .* DeltaK); EpsilonK(jdx, idx) = sqrt(((PlanckConstantReduced^2 .* (Transf.kx(idx).^2 + Transf.ky(jdx).^2)) ./ (2 * Dy164Mass)) .* DeltaK);
end end
end end
EpsilonK = double(imag(EpsilonK) ~= 0); % 'isreal' returns 0 for complex numbers and 1 for real numbers, so we negate it EpsilonK = real(EpsilonK);
% Plot the result % Plot the result
figure(10) figure(14)
clf clf
set(gcf,'Position',[50 50 950 750]) set(gcf,'Position',[50 50 950 750])
imagesc(fftshift(Transf.kx)*1e-6, fftshift(Transf.ky)*1e-6, EpsilonK); % Specify x and y data for axes imagesc(fftshift(Transf.kx)*1e-6, fftshift(Transf.ky)*1e-6, fftshift(EpsilonK)); % Specify x and y data for axes
set(gca, 'YDir', 'normal'); % Correct the y-axis direction set(gca, 'YDir', 'normal'); % Correct the y-axis direction
cbar1 = colorbar; cbar1 = colorbar;
cbar1.Label.Interpreter = 'latex'; cbar1.Label.Interpreter = 'latex';
xlabel('$k_x l_o$','fontsize',16,'interpreter','latex'); xlabel('$k_x l_o$','fontsize',16,'interpreter','latex');
ylabel('$k_y l_o$','fontsize',16,'interpreter','latex'); ylabel('$k_y l_o$','fontsize',16,'interpreter','latex');
title(['2-D Dispersion: $\theta = ',num2str(Params.theta), '; \eta = ', num2str(Params.eta),'$'],'fontsize',16,'interpreter','latex') title(['2-D Dispersion: $\theta = ',num2str(Params.theta), '; \phi = ', num2str(Params.phi),'$'],'fontsize',16,'interpreter','latex')
% Capture the frame and write to video % Capture the frame and write to video
% frame = getframe(gcf); % Capture the current figure % frame = getframe(gcf); % Capture the current figure
@ -224,7 +252,7 @@ function VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth)
Qsq = QX.^2 + QY.^2; Qsq = QX.^2 + QY.^2;
absQ = sqrt(Qsq); absQ = sqrt(Qsq);
QDsq = QX.^2*cos(Params.eta)^2 + QY.^2*sin(Params.eta)^2; % eta here is the azimuthal angle of the dipole moment (angle from the x-axis) QDsq = QX.^2*cos(Params.phi)^2 + QY.^2*sin(Params.phi)^2; % phi here is the azimuthal angle of the dipole moment (angle from the x-axis)
% Bare interaction % Bare interaction
Fpar = -1 + 3*sqrt(pi)*QDsq.*erfcx(absQ)./absQ; % Scaled complementary error function erfcx(x) = e^(x^2) * erfc(x) Fpar = -1 + 3*sqrt(pi)*QDsq.*erfcx(absQ)./absQ; % Scaled complementary error function erfcx(x) = e^(x^2) * erfc(x)