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

More corrections.

Browse Source
This commit is contained in:
Karthik 2024-11-13 14:06:07 +01:00
parent 4f8d7614a9
commit f8b6d72f52
1 changed files with 10 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -5,7 +5,7 @@ gs = Params.gs;
gdd = Params.gdd;
gammaQF = Params.gammaQF;
KEop = 0.5*kx.^2 + 0.5*ky.^2;
KEop = 0.5*KX.^2 + 0.5*KY.^2;
Ez = (0.25*VParams.sigma^2) + (0.25*wz_tilde^2*VParams.sigma^2);
muchem_tilde = muchem - Ez;
@ -23,15 +23,15 @@ Phi = real(ifftn(frho.*VDk));
% Operators
H = @(w) real(ifft(KEop.*fft(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);
X = @(w,psi) (psi.*real(ifft(VDk.*fft(psi.*w)))) + (3/2)*(gQF_pf_2D*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);
X = @(w,psi) (psi.*real(ifft(VDk.*fft(psi.*w)))) + (3/2)*(gQF_pf_2D*gammaQF*abs(psi).^3).*w;
% Operate on g
% Operate in order on g
BdG = @(g) muHC(muHC(g) + (2.*X(g)));
syssize = size(psi);
opts.v0 = psi(:);
opts.v0 = psi(:,:);
opts.tol = 1e-16;
opts.disp = 1;
opts.issym = 0;
@ -46,17 +46,10 @@ clear D;
% Eigenvalues
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) = (1/evals(ii)) * (muHC(g(:,ii)) + (2.*X(g(:,ii))));
f(:,:,ii) = (1/evals(ii)) * (muHC(g(:,:,ii)) + (2.*X(g(:,:,ii))));
end
% Obtain u and v from f and g
@ -65,9 +58,9 @@ 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);
normalization = sum(abs(u(:,:,ii)).^2 - abs(v(:,:,ii)).^2);
u(:,:,ii) = u(:,:,ii)/sqrt(normalization);
v(:,:,ii) = v(:,:,ii)/sqrt(normalization);
end
modes.u = u'; modes.v = v';