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

Corrected pre-factors, added lines to get us and vs with correct normalization.

Browse Source
This commit is contained in:
Karthik 2024-11-12 20:16:09 +01:00
parent e08cd00cac
commit 567193129b
1 changed files with 33 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

37
Bogoliubov-deGennes-Solver/solveBogoliubovdeGennesIn2D.m
View File

@ -9,8 +9,8 @@ KEop = 0.5*kx.^2 + 0.5*ky.^2;
Ez = (0.25*VParams.sigma^2) + (0.25*wz_tilde^2*VParams.sigma^2); Ez = (0.25*VParams.sigma^2) + (0.25*wz_tilde^2*VParams.sigma^2);
muchem_tilde = muchem - Ez; muchem_tilde = muchem - Ez;
gs_pf_2D = 1/(sqrt(2*pi)*VParams.sigma); g_pf_2D = 1/(sqrt(2*pi)*VParams.sigma);
gdd_pf_2D = sqrt(2/5)/(pi^(3/4)*VParams.sigma^(3/2)); gQF_pf_2D = sqrt(2/5)/(pi^(3/4)*VParams.sigma^(3/2));
% eigs only works with column vectors % eigs only works with column vectors
psi = psi.'; psi = psi.';
@ -21,11 +21,13 @@ VDk = VDk.';
frho = fftn(abs(psi).^2); frho = fftn(abs(psi).^2);
Phi = real(ifftn(frho.*VDk)); Phi = real(ifftn(frho.*VDk));
% Operators
H = @(w) real(ifft(KEop.*fft(w))); H = @(w) real(ifft(KEop.*fft(w)));
C = @(w) (((gs_pf_2D*gs*abs(psi).^2) + (gdd_pf_2D*gdd*Phi)).*w) + (gammaQF.* abs(psi).^3 .*w); C = @(w) (((g_pf_2D*gs*abs(psi).^2) + (g_pf_2D*gdd*Phi)).*w) + (gQF_pf_2D*gammaQF.*abs(psi).^3.*w);
muHC = @(w) (-muchem_tilde .* w) + H(w) + C(w); muHC = @(w) (-muchem_tilde .* w) + H(w) + C(w);
X = @(w,psi) (psi.*real(ifft(VDk.*fft(psi.*w)))) + (3/2)*(gammaQF.* abs(psi).^3).*w; X = @(w,psi) (psi.*real(ifft(VDk.*fft(psi.*w)))) + (3/2)*(gQF_pf_2D*gammaQF.* abs(psi).^3).*w;
% Operate on g
BdG = @(g) muHC(muHC(g) + (2.*X(g))); BdG = @(g) muHC(muHC(g) + (2.*X(g)));
syssize = size(psi); syssize = size(psi);
@ -44,5 +46,32 @@ clear D;
% Eigenvalues % Eigenvalues
evals = sqrt(evals); evals = sqrt(evals);
BogModeCutoff = 1;
if BogModeCutoff == 1
g(:,real(evals)<1e-6) = [];
evals(real(evals)<1e-6) = [];
Neigs = length(evals);
end
% Obtain f from g
f = zeros(size(g));
for ii = 1:Neigs
f(:,ii) = muHC(g(:,ii)) + (2.*X(g(:,ii)));
end
% Obtain u and v from f and g
u = (f + g)/2;
v = (f - g)/2;
% Renormalize to \int |u|^2 - |v|^2 = 1
for ii=1:Neigs
normalization = sum(abs(u(:,ii)).^2 - abs(v(:,ii)).^2)*Transf.dz;
u(:,ii) = u(:,ii)/sqrt(normalization);
v(:,ii) = v(:,ii)/sqrt(normalization);
end
modes.u = u'; modes.v = v';
modes.g = g'; modes.f = f';
end end