From 4583f1a6439f3e5279b8caf10d365c61f5ee2bbf Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Mon, 5 May 2025 18:49:45 +0200 Subject: [PATCH] Corrected script for calculating the dispersion relation in 2D - thanks to Tom Bland! --- .gitignore | 2 + .../RotonInstability/DipolarDispersion2D.m | 84 ++++++++++++------- 2 files changed, 58 insertions(+), 28 deletions(-) diff --git a/.gitignore b/.gitignore index ab48627..632af07 100644 --- a/.gitignore +++ b/.gitignore @@ -10,5 +10,7 @@ Time-Series-Analyzer/Time-Series-Data *.mp4 *.bat *.json +*.txt +*.asv .ipynb_checkpoints/ .vscode/ \ No newline at end of file diff --git a/Estimations/RotonInstability/DipolarDispersion2D.m b/Estimations/RotonInstability/DipolarDispersion2D.m index b487963..dc07a5a 100644 --- a/Estimations/RotonInstability/DipolarDispersion2D.m +++ b/Estimations/RotonInstability/DipolarDispersion2D.m @@ -31,7 +31,7 @@ DyMagneticMoment = 9.93*BohrMagneton; %% 2-D DDI Potential in k-space, with Gaussian ansatz width determined by constrained minimization % 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)); % % Defining a harmonic oscillator length - heredue to the choice of w0, l0 % is 1 micrometer @@ -40,8 +40,8 @@ wz = 2 * pi * 300; lz = sqrt(PlanckConstantReduced/(Dy164Mass * wz)); % Defining a harmonic oscillator length % Number of grid points in each direction -Params.Nx = 128; -Params.Ny = 128; +Params.Nx = 256; +Params.Ny = 256; Params.Lx = 150*1e-6; Params.Ly = 150*1e-6; [Transf] = setupSpace(Params); @@ -49,7 +49,7 @@ Params.Ly = 150*1e-6; nadd2s = 0.110; % Number density * add^2 as_to_add = 0.782; % 1/edd 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 gdd = VacuumPermeability*DyMagneticMoment^2/3; @@ -74,19 +74,19 @@ sigma = fmincon(TotalEnergyPerParticle, x0, Aineq, Bineq, Aeq, MeanWidth = sigma * lz; % == 2-D DDI Potential in k-space == % -VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth); -VDk_fftshifted = fftshift(VDk); +VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth); % Transf.K is in shifted space, so VDk is shifted +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 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 -set(gca, 'YDir', 'normal'); % Correct the y-axis direction +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 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(['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 == % 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)); gs_tilde = gs / (sqrt(2*pi) * MeanWidth); +gdd_tilde = gdd / (sqrt(2*pi) * MeanWidth); % you were missing this! % == 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)); + 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); 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 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 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(['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 -% 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) -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 % 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 theta and eta values +% Loop over theta and phi values for theta = theta_values - for eta = eta_values - % Update Params with current theta and eta + for phi = phi_values + % Update Params with current theta and phi 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 == % VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth); VDk_fftshifted = fftshift(VDk); @@ -147,28 +174,29 @@ for theta = theta_values EpsilonK = zeros(length(Transf.ky), length(Transf.kx)); gs_tilde = gs / (sqrt(2*pi) * MeanWidth); + gdd_tilde = gdd / (sqrt(2*pi) * MeanWidth); % you were missing this! % == 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)); + 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); 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 - figure(10) + figure(14) clf 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 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(['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 % frame = getframe(gcf); % Capture the current figure @@ -224,7 +252,7 @@ function VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth) Qsq = QX.^2 + QY.^2; 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 Fpar = -1 + 3*sqrt(pi)*QDsq.*erfcx(absQ)./absQ; % Scaled complementary error function erfcx(x) = e^(x^2) * erfc(x)