From 6093c700aba6d7c17039695663f799e85ed8bc21 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Fri, 10 Jan 2025 15:51:44 +0100 Subject: [PATCH] Computes for specified edd and/or density. --- .../DipolarDispersion2D.m | 76 +++++-------------- 1 file changed, 18 insertions(+), 58 deletions(-) diff --git a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersion2D.m b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersion2D.m index 7113df3..95a0a79 100644 --- a/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersion2D.m +++ b/Estimations/DipolarDispersionAndRotonInstabilityBoundary/DipolarDispersion2D.m @@ -30,29 +30,27 @@ DyMagneticMoment = 9.93*BohrMagneton; %% 2-D DDI Potential in k-space, with Gaussian ansatz width determined by constrained minimization +w0 = 2*pi*61.6316; % Angular frequency unit [s^-1] +l0 = sqrt(PlanckConstantReduced/(Dy164Mass*w0)); % Defining a harmonic oscillator length + 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; +nadd2s = 0.110; +as_to_add = 0.782; +Params.alpha = 10; +Params.phi = 0; 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 = []; @@ -63,42 +61,16 @@ 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'); +AtomNumberDensity = nadd2s / add^2; % Areal density of atoms +as = as_to_add * add; % Scattering length +eps_dd = add/as; % Relative interaction strength +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); %% 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; +MeanWidth = sigma * lz; % == 2-D DDI Potential in k-space == % VDk = compute2DPotentialInMomentumSpace(Transf, Params, MeanWidth); @@ -108,7 +80,7 @@ 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 +set(gca, 'YDir', 'normal'); % Correct the y-axis direction cbar1 = colorbar; cbar1.Label.Interpreter = 'latex'; xlabel('$k_x l_o$','fontsize',16,'interpreter','latex'); @@ -131,13 +103,13 @@ for idx = 1:length(Transf.kx) end end -EpsilonK = double(imag(EpsilonK) ~= 0); % 'isreal' returns 0 for complex numbers and 1 for real numbers, so we negate it +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 +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'); @@ -146,22 +118,10 @@ title(['$\theta = ',num2str(Params.alpha), '; \phi = ', num2str(Params.phi),'$'] %% -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